Analyzing single-cell RNA-seq data with Bioconductor (miscellaneous)

2.1 Using fixed thresholds

One alternative strategy is to set pre-defined thresholds on each QC metric. For example, we might remove all cells with library sizes below 100000 and numbers of expressed genes below 4000. This generally requires substantial experience to determine appropriate thresholds for each experimental protocol and biological system. Indeed, even with the same protocol and system, the appropriate threshold can vary from run to run due to the vagaries of RNA capture and sequencing.

≣ Table of Contents
◂ Using fixed thresholds
▴ Alternative approaches to quality control

2.2 Using PCA-based outliers

Another strategy is to perform a principal components analysis (PCA) based on the quality metrics for each cell, e.g., the total number of reads, the total number of features and the proportion of mitochondrial or spike-in reads. Outliers on a PCA plot may be indicative of low-quality cells that have aberrant technical properties compared to the (presumed) majority of high-quality cells. This is demonstrated below on a brain cell dataset from Tasic et al. (2016), using functions from the scater package (???).

# Obtaining the dataset.
library(scRNAseq)
data(allen)

# Setting up the data.
library(SingleCellExperiment)
sce.allen <- as(allen, "SingleCellExperiment")
assayNames(sce.allen) <- "counts"
isSpike(sce.allen, "ERCC") <- grep("ERCC", rownames(sce.allen))

# Computing the QC metrics and running PCA.
library(scater)
sce.allen <- calculateQCMetrics(sce.allen, feature_controls=list(ERCC=isSpike(sce.allen)))
sce.allen <- runPCA(sce.allen, pca_data_input="coldata", detect_outliers=TRUE)
table(sce.allen$outlier)

## 
## FALSE  TRUE 
##   374     5

Methods like PCA-based outlier detection and support vector machines can provide more power to distinguish low-quality cells from high-quality counterparts (Ilicic et al. 2016). This is because they are able to detect subtle patterns across many quality metrics simultaneously. However, this comes at some cost to interpretability, as the reason for removing a given cell may not always be obvious. Users interested in the more sophisticated approaches are referred to the scater and cellity packages.

For completeness, we note that outliers can also be identified from PCA on the gene expression profiles, rather than QC metrics. We consider this to be a risky strategy as it can remove high-quality cells in rare populations.

3.1 Motivation

Scaling normalization strategies for scRNA-seq data can be broadly divided into two classes. The first class assumes that there exists a subset of genes that are not DE between samples, as previously described. The second class uses the fact that the same amount of spike-in RNA was added to each cell (A. T. L. Lun et al. 2017). Differences in the coverage of the spike-in transcripts can only be due to cell-specific biases, e.g., in capture efficiency or sequencing depth. Scaling normalization is then applied to equalize spike-in coverage across cells.

The choice between these two normalization strategies depends on the biology of the cells and the features of interest. If the majority of genes are expected to be DE and there is no reliable house-keeping set, spike-in normalization may be the only option for removing cell-specific biases. Spike-in normalization should also be used if differences in the total RNA content of individual cells are of interest. In any particular cell, an increase in the amount of endogenous RNA will not increase spike-in coverage (with or without library quantification). Thus, the former will not be represented as part of the bias in the latter, which means that the effects of total RNA content on expression will not be removed upon scaling. With non-DE normalization, an increase in RNA content will systematically increase the expression of all genes in the non-DE subset, such that it will be treated as bias and removed.

≣ Table of Contents
▸ Calculating spike-in size factors
◂ Motivation
▴ Normalizing based on spike-in coverage

3.2 Setting up the data

We demonstrate the use of spike-in normalization on a dataset involving different cell types – namely, mouse embryonic stem cells (mESCs) and mouse embryonic fibroblasts (MEFs) (Islam et al. 2011). The count table was obtained from the NCBI Gene Expression Omnibus (GEO) as a supplementary file using the accession number GSE29087. We load the counts into R, using colClasses to speed up read.table by pre-defining the type of each column. We also specify the rows corresponding to spike-in transcripts.

library(SingleCellExperiment)
counts <- read.table("GSE29087_L139_expression_tab.txt.gz", 
    colClasses=c(list("character", NULL, NULL, NULL, NULL, NULL, NULL), 
        rep("integer", 96)), skip=6, sep='\t', row.names=1)
is.spike <- grep("SPIKE", rownames(counts)) 
sce.islam <- SingleCellExperiment(list(counts=as.matrix(counts)))
isSpike(sce.islam, "spike") <- is.spike
dim(sce.islam)

## [1] 22936    96

We perform some quality control to remove low-quality cells using the calculateQCMetrics function. Outliers are identified within each cell type to avoid issues with systematic differences in the metrics between cell types. The negative control wells do not contain any cells and are useful for quality control (as they should manifest as outliers for the various metrics), but need to be removed prior to downstream analysis.

library(scater)
sce.islam <- calculateQCMetrics(sce.islam, feature_controls=list(spike=is.spike))
sce.islam$grouping <- rep(c("mESC", "MEF", "Neg"), c(48, 44, 4))

libsize.drop <- isOutlier(sce.islam$total_counts, nmads=3, type="lower", 
    log=TRUE, batch=sce.islam$grouping)
feature.drop <- isOutlier(sce.islam$total_features, nmads=3, type="lower", 
    log=TRUE, batch=sce.islam$grouping)
spike.drop <- isOutlier(sce.islam$pct_counts_spike, nmads=3, type="higher", 
    batch=sce.islam$grouping)
    
sce.islam <- sce.islam[,!(libsize.drop | feature.drop | 
    spike.drop | sce.islam$grouping=="Neg")]
data.frame(ByLibSize=sum(libsize.drop), ByFeature=sum(feature.drop),
    BySpike=sum(spike.drop), Remaining=ncol(sce.islam))

##   ByLibSize ByFeature BySpike Remaining
## 1         4         6      12        77

≣ Table of Contents
◂ Setting up the data
▴ Normalizing based on spike-in coverage

3.3 Calculating spike-in size factors

We apply the computeSpikeFactors method to estimate size factors for all cells. This method computes the total count over all spike-in transcripts in each cell, and calculates size factors to equalize the total spike-in count across cells. Here, we set general.use=TRUE as we intend to apply the spike-in factors to all counts.

library(scran)
sce.islam <- computeSpikeFactors(sce.islam, general.use=TRUE)

Running normalize will use the spike-in-based size factors to compute normalized log-expression values. Unlike the previous analyses, we do not have to define separate size factors for the spike-in transcripts. This is because the relevant factors are already being used for all genes and spike-in transcripts when general.use=TRUE. (The exception is if the experiment uses multiple spike-in sets that behave differently and need to be normalized separately.)

sce.islam <- normalize(sce.islam)

For comparison, we also compute the deconvolution size factors (A. T. Lun, Bach, and Marioni 2016) and plot them against the spike-in factors. We observe a negative correlation between the two sets of values (Figure 1). This is because MEFs contain more endogenous RNA, which reduces the relative spike-in coverage in each library (thereby decreasing the spike-in size factors) but increases the coverage of endogenous genes (thus increasing the deconvolution size factors). If the spike-in size factors were applied to the counts, the expression values in MEFs would be scaled up while expression in mESCs would be scaled down. However, the opposite would occur if deconvolution size factors were used.

colours <- c(mESC="red", MEF="grey")
deconv.sf <- computeSumFactors(sce.islam, sf.out=TRUE, cluster=sce.islam$grouping)
plot(sizeFactors(sce.islam), deconv.sf, col=colours[sce.islam$grouping], pch=16, 
    log="xy", xlab="Size factor (spike-in)", ylab="Size factor (deconvolution)")
legend("bottomleft", col=colours, legend=names(colours), pch=16)

Figure 1: Size factors from spike-in normalization, plotted against the size factors from deconvolution for all cells in the mESC/MEF dataset. Axes are shown on a log-scale, and cells are coloured according to their identity. Deconvolution size factors were computed with small pool sizes owing to the low number of cells of each type.

Whether or not total RNA content is relevant – and thus, the choice of normalization strategy – depends on the biological hypothesis. In the HSC and brain analyses, variability in total RNA across the population was treated as noise and removed by non-DE normalization. This may not always be appropriate if total RNA is associated with a biological difference of interest. For example, Islam et al. (2011) observe a 5-fold difference in total RNA between mESCs and MEFs. Similarly, the total RNA in a cell changes across phases of the cell cycle (Buettner et al. 2015). Spike-in normalization will preserve these differences in total RNA content such that the corresponding biological groups can be easily resolved in downstream analyses.

Comments from Aaron:

• We only use genes with average counts greater than 1 (as specified in min.mean) to compute the deconvolution size factors. This avoids problems with discreteness as mentioned in our previous uses of computeSumFactors.
• Setting sf.out=TRUE will directly return the size factors, rather than a SingleCellExperiment object containing those factors. This is more convenient when only the size factors are required for further analysis.

4.1 Setting up the data

Highly variable genes (HVGs) are defined as genes with biological components that are significantly greater than zero. These genes are interesting as they drive differences in the expression profiles between cells, and should be prioritized for further investigation. Formal detection of HVGs allows us to avoid genes that are highly variable due to technical factors such as sampling noise during RNA capture and library preparation. This adds another level of statistical rigour to our previous analyses, in which we only modelled the technical component.

To demonstrate, we use data from haematopoietic stem cells (HSCs) (Wilson et al. 2015), generated using the Smart-seq2 protocol (Picelli et al. 2014) with ERCC spike-ins. Counts were obtained from NCBI GEO as a supplementary file using the accession number GSE61533. Our first task is to load the count matrix into memory. In this case, some work is required to retrieve the data from the Gzip-compressed Excel format.

library(R.utils)
gunzip("GSE61533_HTSEQ_count_results.xls.gz", remove=FALSE, overwrite=TRUE)
library(readxl)
all.counts <- as.data.frame(read_excel('GSE61533_HTSEQ_count_results.xls', sheet=1))
rownames(all.counts) <- all.counts$ID
all.counts <- as.matrix(all.counts[,-1])

We store the results in a SingleCellExperiment object and identify the rows corresponding to the spike-ins based on the row names.

sce.hsc <- SingleCellExperiment(list(counts=all.counts))
dim(sce.hsc)

## [1] 38498    96

is.spike <- grepl("^ERCC", rownames(sce.hsc))
isSpike(sce.hsc, "ERCC") <- is.spike
summary(is.spike)

##    Mode   FALSE    TRUE 
## logical   38406      92

For each cell, we calculate quality control metrics using the calculateQCMetrics function as previously described. We filter out HSCs that are outliers for any metric, under the assumption that these represent low-quality libraries.

sce.hsc <- calculateQCMetrics(sce.hsc, feature_controls=list(ERCC=is.spike))
libsize.drop <- isOutlier(sce.hsc$total_counts, nmads=3, type="lower", log=TRUE)
feature.drop <- isOutlier(sce.hsc$total_features, nmads=3, type="lower", log=TRUE)
spike.drop <- isOutlier(sce.hsc$pct_counts_ERCC, nmads=3, type="higher")
sce.hsc <- sce.hsc[,!(libsize.drop | feature.drop | spike.drop)]
data.frame(ByLibSize=sum(libsize.drop), ByFeature=sum(feature.drop),
    BySpike=sum(spike.drop), Remaining=ncol(sce.hsc))

##   ByLibSize ByFeature BySpike Remaining
## 1         2         2       3        92

We remove genes that are not expressed in any cell to reduce computational work in downstream steps.

to.keep <- nexprs(sce.hsc, byrow=TRUE) > 0
sce.hsc <- sce.hsc[to.keep,]
summary(to.keep)

##    Mode   FALSE    TRUE 
## logical   17229   21269

We apply the deconvolution method to compute size factors for the endogenous genes (A. T. Lun, Bach, and Marioni 2016). Separate size factors for the spike-in transcripts are also calculated, as previously discussed. We then calculate log-transformed normalized expression values for further use.

sce.hsc <- computeSumFactors(sce.hsc)
summary(sizeFactors(sce.hsc))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4076  0.8101  0.9546  1.0000  1.1585  2.0264

sce.hsc <- computeSpikeFactors(sce.hsc, type="ERCC", general.use=FALSE)
summary(sizeFactors(sce.hsc, "ERCC"))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2562  0.6198  0.8623  1.0000  1.2122  3.0289

sce.hsc <- normalize(sce.hsc)

≣ Table of Contents
◂ Setting up the data
▴ Detecting highly variable genes

4.2 Testing for significantly positive biological components

We fit a mean-variance trend to the spike-in transcripts to quantify the technical component of the variance, as previously described. The biological component for each gene is defined as the difference between its total variance and the fitted value of the trend (Figure 2).

var.fit <- trendVar(sce.hsc, parametric=TRUE, span=0.3)
var.out <- decomposeVar(sce.hsc, var.fit)
plot(var.out$mean, var.out$total, pch=16, cex=0.6, xlab="Mean log-expression", 
    ylab="Variance of log-expression")
curve(var.fit$trend(x), col="dodgerblue", lwd=2, add=TRUE)
cur.spike <- isSpike(sce.hsc)
points(var.out$mean[cur.spike], var.out$total[cur.spike], col="red", pch=16)

Figure 2: Variance of normalized log-expression values for each gene in the HSC dataset, plotted against the mean log-expression. The blue line represents the mean-dependent trend fitted to the variances of the spike-in transcripts (red).

We define HVGs as those genes that have a biological component that is significantly greater than zero. We use a false discovery rate (FDR) of 5% after correcting for multiple testing with the Benjamini-Hochberg method.

hvg.out <- var.out[which(var.out$FDR <= 0.05),]
nrow(hvg.out)

## [1] 507

We rank the results to focus on genes with larger biological components. This highlights an interesting aspect of the underlying hypothesis test, which is based on the ratio of the total variance to the expected technical variance. Ranking based on p-value tends to prioritize HVGs that are more likely to be true positives but, at the same time, less likely to be interesting. This is because the ratio can be very large for HVGs that have very low total variance and do not contribute much to the cell-cell heterogeneity.

hvg.out <- hvg.out[order(hvg.out$bio, decreasing=TRUE),] 
write.table(file="hsc_hvg.tsv", hvg.out, sep="\t", quote=FALSE, col.names=NA)
head(hvg.out)

##              mean    total       bio      tech      p.value          FDR
## Fos      6.461286 19.56476 12.805943  6.758821 1.119519e-18 7.660040e-16
## Dusp1    6.821227 15.62737 10.101224  5.526149 8.374847e-18 4.801051e-15
## Rgs1     5.311914 20.30772 10.004682 10.303040 9.682171e-08 1.134633e-05
## Ppp1r15a 6.666806 14.52844  8.481203  6.047234 1.657304e-12 4.687077e-10
## Ly6a     8.404743 10.07457  8.075840  1.998731 1.917314e-50 5.809736e-47
## Egr1     6.714811 13.84811  7.964822  5.883286 6.190503e-12 1.601302e-09

We check the distribution of expression values for the genes with the largest biological components. This ensures that the variance estimate is not driven by one or two outlier cells (Figure 3).

fontsize <- theme(axis.text=element_text(size=12), axis.title=element_text(size=16))
plotExpression(sce.hsc, features=rownames(hvg.out)[1:10]) + fontsize

Figure 3: Violin plots of normalized log-expression values for the top 10 genes with the largest biological components in the HSC dataset. Each point represents the log-expression value in a single cell.

There are many other strategies for defining HVGs, based on a variety of metrics:

• the coefficient of variation (Brennecke et al. 2013; Kołodziejczyk et al. 2015; Kim et al. 2015)
• the dispersion parameter in the negative binomial distribution (McCarthy, Chen, and Smyth 2012)
• a proportion of total variability (Vallejos, Marioni, and Richardson 2015)

Some of these methods are available in scran – for example, see DM or technicalCV2 for calculations based on the coefficient of variation. Here, we use the variance of the log-expression values because the log-transformation protects against genes with strong expression in only one or two cells. This ensures that the set of top HVGs is not dominated by genes with (mostly uninteresting) outlier expression patterns.

5.1 Trend fitting when spike-ins are unavailable

If spike-in RNA has not been added in appropriate quantities (or at all), an alternative approach is to fit the trend to the variance estimates of the endogenous genes. This is done using the use.spikes=FALSE setting in trendVar, as shown below for the HSC dataset.

var.fit.nospike <- trendVar(sce.hsc, parametric=TRUE, use.spikes=FALSE, span=0.2)
var.out.nospike <- decomposeVar(sce.hsc, var.fit.nospike)

The simplest interpretation of the results assumes that the majority of genes are not variably expressed. This means that the technical component dominates the total variance for most genes, such that the fitted trend can be treated as an estimate of the technical component. In Figure 11, the trend passes through or close to most of the spike-in variances, indicating that our assumption is valid.

plot(var.out.nospike$mean, var.out.nospike$total, pch=16, cex=0.6, 
    xlab="Mean log-expression", ylab="Variance of log-expression")
curve(var.fit.nospike$trend(x), col="dodgerblue", lwd=2, add=TRUE)
points(var.out.nospike$mean[cur.spike], var.out.nospike$total[cur.spike], col="red", pch=16)

Figure 4: Variance of normalized log-expression values for each gene in the 416B dataset, plotted against the mean log-expression. The blue line represents the mean-dependent trend fitted to the variances of the endogenous genes (black), with spike-in transcripts shown in red.

If our assumption does not hold, the output of decomposeVar is more difficult to interpret. The fitted value of the trend can no longer be generally interpreted as the technical component, as it contains some biological variation as well. Instead, recall that the biological component reported by decomposeVar represents the residual for each gene over the majority of genes with the same abundance. One could assume that the variabilities of most genes are driven by constitutive “house-keeping” processes, which are biological in origin but generally uninteresting. Any gene with an increase in its variance is relatively highly variable and can be prioritized for further study.

≣ Table of Contents
◂ Trend fitting when spike-ins are unavailable
▴ Advanced modelling of the technical noise

5.2 Blocking on uninteresting factors of variation

Our previous analysis of the 416B dataset specified design in trendVar to ensure that systematic differences between plates do not inflate the variance. However, this implicitly assumes that the trend is the same between plates, given that a single trend is fitted to the spike-in transcripts. This may not always be the case, e.g., when different amounts of spike-in RNA are added between batches. The use of a single trend would subsequently be inappropriate, resulting in inaccurate estimates of the technical component for each gene.

For datasets containing multiple batches, an alternative strategy is to perform trend fitting and variance decomposition separately for each batch. This accommodates differences in the mean-variance trends between batches, especially if a different amount of spike-in RNA was added to the cells in each batch. We demonstrate this approach by treating each plate in the 416B dataset as a different batch (Figure 5). This yields plate-specific estimates of the biological and technical components for each gene.

sce.416B <- readRDS("416B_data.rds") # Loading the saved object.
collected <- list()
par(mfrow=c(1,2))
for (plate in levels(sce.416B$Plate)) {
    cur.sce <- sce.416B[,sce.416B$Plate==plate]
    cur.sce <- normalize(cur.sce) 

    # Estimating the technical/biological components.
    cur.fit <- trendVar(cur.sce, parametric=TRUE, span=0.4)
    cur.out <- decomposeVar(cur.sce, cur.fit)
    collected[[plate]] <- cur.out 

    # Making a plot for diagnostic purposes.
    plot(cur.out$mean, cur.out$total, pch=16, cex=0.6, xlab="Mean log-expression", 
        ylab="Variance of log-expression", main=plate)
    curve(cur.fit$trend(x), col="dodgerblue", lwd=2, add=TRUE)
    cur.spike <- isSpike(cur.sce)
    points(cur.out$mean[cur.spike], cur.out$total[cur.spike], col="red", pch=16)
}

Figure 5: Variance of normalized log-expression values for each gene in each plate of the 416B dataset, plotted against the mean log-expression. The blue line represents the mean-dependent trend fitted to the variances of the spike-in transcripts (red).

Statistics are combined across multiple batches using the combineVar function. This function uses a weighted average across batches for the means and variances, and Fisher’s method for combining the p-values. These results can be used in downstream functions such as denoisePCA, or for detecting highly variable genes (see below).

comb.out <- do.call(combineVar, collected)
head(comb.out)

##                           mean       total         bio       tech   p.value
## ENSMUSG00000103377 0.008222727 0.012305609 -0.02679672 0.03910233       NaN
## ENSMUSG00000103147 0.033848916 0.069429347 -0.07460563 0.14403498       NaN
## ENSMUSG00000103161 0.005070313 0.004730286 -0.01684507 0.02157536       NaN
## ENSMUSG00000102331 0.018625989 0.032999668 -0.05186855 0.08486822 1.0000000
## ENSMUSG00000102948 0.058106895 0.085618238 -0.16558431 0.25120255 1.0000000
## Rp1                0.097109613 0.457054312  0.02088991 0.43616440 0.0934758
##                          FDR
## ENSMUSG00000103377       NaN
## ENSMUSG00000103147       NaN
## ENSMUSG00000103161       NaN
## ENSMUSG00000102331 1.0000000
## ENSMUSG00000102948 1.0000000
## Rp1                0.2528314

Comments from Aaron:

• We run normalize within the loop to ensure that the average abundances of the spike-in transcripts are comparable to the endogenous genes. This adjusts the size factors across cells in each batch so that the mean is equal to 1, for both the spike-in and gene-based sets of size factors. Log-normalized expression values are then recalculated using these centred size factors. This procedure avoids problems due to differences in the quantity of spike-in RNA between batches. In such cases, if the globally-centred size factors were used, there would be a systematic difference in the scaling of spike-in transcripts compared to endogenous genes. The fitted trend would then be shifted along the x-axis and fail to accurately capture the technical component for each gene.