About This Document »

Package name simpleSingleCell
Built with Bioconductor (R) 3.6 (3.4.2)
Last Built Fri, 15 Dec 2017 15:26:00 -0800
Last Modified Fri, 15 Dec 2017 14:47:50 -0800 (r132255)
Source Package ../simpleSingleCell_1.0.5.tar.gz
Windows Binary ../NOT_SUPPORTED
R Script part3.R

To install this workflow under Bioconductor 3.6, start R and enter:

source("http://bioconductor.org/workflows.R")
workflowInstall("simpleSingleCell")

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

Aaron T. L. Lun1, Davis J. McCarthy2,3 and John C. Marioni1,2,4

1Cancer Research UK Cambridge Institute, Li Ka Shing Centre, Robinson Way, Cambridge CB2 0RE, United Kingdom
2EMBL European Bioinformatics Institute, Wellcome Genome Campus, Hinxton, Cambridge CB10 1SD, United Kingdom
3St Vincent's Institute of Medical Research, 41 Victoria Parade, Fitzroy, Victoria 3065, Australia
4Wellcome Trust Sanger Institute, Wellcome Genome Campus, Hinxton, Cambridge CB10 1SA, United Kingdom

Table of Contents
Alternative approaches to quality control
Top

1 Overview

The previous workflows focused on analyzing single-cell RNA-seq data with “standard” procedures. However, a number of alternative parameter settings and strategies can be used at some steps of the workflow. This workflow describes a few of these alternative settings as well as the rationale behind choosing them instead of the defaults.

Table of Contents
Normalizing based on spike-in coverage
Overview
Top
Using PCA-based outliers

2 Alternative approaches to quality control

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.

Table of Contents
Detecting highly variable genes
Alternative approaches to quality control
Top
Setting up the data

3 Normalizing based on spike-in coverage

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)
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.

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.

Table of Contents
Advanced modelling of the technical noise
Normalizing based on spike-in coverage
Top
Testing for significantly positive biological components

4 Detecting highly variable genes

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)
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).

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
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.

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.

Table of Contents
Identifying correlated gene pairs with Spearman’s rho
Detecting highly variable genes
Top
Blocking on uninteresting factors of variation

5 Advanced modelling of the technical noise

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)
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.

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)
}
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).

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.

Table of Contents
Blocking on the cell cycle phase
Advanced modelling of the technical noise
Top

6 Identifying correlated gene pairs with Spearman’s rho

Another use for scRNA-seq data is to identify correlations between the expression profiles of different genes. This is quantified by computing Spearman’s rho, which accommodates non-linear relationships in the expression values. Non-zero correlations between pairs of genes provide evidence for their co-regulation. However, the noise in the data requires some statistical analysis to determine whether a correlation is significantly non-zero.

To demonstrate, we use the correlatePairs function to identify significant correlations between the various histocompatability antigens in the HSC data set. The significance of each correlation is determined using a permutation test. For each pair of genes, the null hypothesis is that the expression profiles of two genes are independent. Shuffling the profiles and recalculating the correlation yields a null distribution that is used to obtain a p-value for each observed correlation value (Phipson and Smyth 2010).

set.seed(100)
var.cor <- correlatePairs(sce.hsc, subset.row=grep("^H2-", rownames(sce.hsc)))
head(var.cor)
##    gene1   gene2       rho      p.value          FDR limited
## 1 H2-Ab1  H2-Eb1 0.4968944 1.999998e-06 0.0004349996    TRUE
## 2  H2-Aa  H2-Ab1 0.4886488 1.999998e-06 0.0004349996    TRUE
## 3  H2-D1   H2-K1 0.4150548 2.999997e-05 0.0043499957   FALSE
## 4  H2-Aa  H2-Eb1 0.4093985 4.199996e-05 0.0045674954   FALSE
## 5 H2-Ab1 H2-DMb1 0.3581986 4.859995e-04 0.0422819577   FALSE
## 6  H2-Q6   H2-Q7 0.3432024 8.659991e-04 0.0627849372   FALSE

Correction for multiple testing across many gene pairs is performed by controlling the FDR at 5%.

sig.cor <- var.cor$FDR <= 0.05
summary(sig.cor)
##    Mode   FALSE    TRUE 
## logical     430       5

We can also compute correlations between specific pairs of genes, or between all pairs between two distinct sets of genes. The example below computes the correlation between Fos and Jun, which dimerize to form the AP-1 transcription factor (Angel and Karin 1991).

correlatePairs(sce.hsc, subset.row=cbind("Fos", "Jun"))
##   gene1 gene2       rho      p.value          FDR limited
## 1   Fos   Jun 0.4634958 5.999994e-06 5.999994e-06   FALSE

Examination of the expression profiles in Figure 6 confirms the presence of a modest correlation between these two genes.

plotExpression(sce.hsc, features="Fos", x="Jun")
Expression of _Fos_ plotted against the expression of _Jun_ for all cells in the HSC data set.

Figure 6: Expression of Fos plotted against the expression of Jun for all cells in the HSC data set.

The use of correlatePairs is primarily intended to identify correlated gene pairs for validation studies. Obviously, non-zero correlations do not provide evidence for a direct regulatory interaction, let alone specify causality. To construct regulatory networks involving many genes, we suggest using dedicated packages such as WCGNA.

Comments from Aaron:

  • We suggest only computing correlations between a subset of genes of interest, known either a priori or empirically defined, e.g., as HVGs. Computing correlations across all genes will take too long; unnecessarily increase the severity of the multiple testing correction; and may prioritize strong but uninteresting correlations, e.g., between tightly co-regulated house-keeping genes.
  • The correlatePairs function can also return gene-centric output by setting per.gene=TRUE. This calculates a combined p-value (Simes 1986) for each gene that indicates whether it is significantly correlated to any other gene. From a statistical perspective, this is a more natural approach to correcting for multiple testing when genes, rather than pairs of genes, are of interest.
  • The Limited field indicates whether the p-value was lower-bounded by the number of permutations. If this is TRUE for any non-significant gene at the chosen FDR threshold, consider increasing the number of permutations to improve power.

Table of Contents
Concluding remarks
Identifying correlated gene pairs with Spearman’s rho
Top

7 Blocking on the cell cycle phase

Cell cycle phase is usually uninteresting in studies focusing on other aspects of biology. However, the effects of cell cycle on the expression profile can mask other effects and interfere with the interpretation of the results. This cannot be avoided by simply removing cell cycle marker genes, as the cell cycle can affect a substantial number of other transcripts (Buettner et al. 2015). Rather, more sophisticated strategies are required, one of which is demonstrated below using data from a study of T Helper 2 (TH2) cells (Mahata et al. 2014). Buettner et al. (2015) have already applied quality control and normalized the data, so we can use them directly as log-expression values (accessible as Supplementary Data 1 of https://dx.doi.org/10.1038/nbt.3102).

incoming <- as.data.frame(read_excel("nbt.3102-S7.xlsx", sheet=1))
rownames(incoming) <- incoming[,1]
incoming <- incoming[,-1]
incoming <- incoming[,!duplicated(colnames(incoming))] # Remove duplicated genes.
sce.th2 <- SingleCellExperiment(list(logcounts=t(incoming)))

We empirically identify the cell cycle phase using the pair-based classifier in cyclone. The majority of cells in Figure 7 seem to lie in G1 phase, with small numbers of cells in the other phases.

library(org.Mm.eg.db)
ensembl <- mapIds(org.Mm.eg.db, keys=rownames(sce.th2), keytype="SYMBOL", column="ENSEMBL")

set.seed(100)
mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds", 
    package="scran"))
assignments <- cyclone(sce.th2, mm.pairs, gene.names=ensembl, assay.type="logcounts")

plot(assignments$score$G1, assignments$score$G2M, 
    xlab="G1 score", ylab="G2/M score", pch=16)
Cell cycle phase scores from applying the pair-based classifier on the T~H~2 dataset, where each point represents a cell.

Figure 7: Cell cycle phase scores from applying the pair-based classifier on the TH2 dataset, where each point represents a cell.

We can block directly on the phase scores in downstream analyses. This is more graduated than using a strict assignment of each cell to a specific phase, as the magnitude of the score considers the uncertainty of the assignment. The phase covariates in the design matrix will absorb any phase-related effects on expression such that they will not affect estimation of the effects of other experimental factors. Users should also ensure that the phase score is not confounded with other factors of interest. For example, model fitting is not possible if all cells in one experimental condition are in one phase, and all cells in another condition are in a different phase.

design <- model.matrix(~ G1 + G2M, assignments$score)
fit.block <- trendVar(sce.th2, design=design, parametric=TRUE, use.spikes=NA)
sce.th2.block <- denoisePCA(sce.th2, technical=fit.block$trend, design=design) 

The result of blocking on design is visualized with some PCA plots in Figure 8. Before removal, the distribution of cells along the first two principal components is strongly associated with their G1 and G2/M scores. This is no longer the case after removal, which suggests that the cell cycle effect has been mitigated.

sce.th2$G1score <- sce.th2.block$G1score <- assignments$score$G1
sce.th2$G2Mscore <- sce.th2.block$G2Mscore <- assignments$score$G2M

# Without blocking on phase score.
fit <- trendVar(sce.th2, parametric=TRUE, use.spikes=NA) 
sce.th2 <- denoisePCA(sce.th2, technical=fit$trend)
out <- plotReducedDim(sce.th2, use_dimred="PCA", ncomponents=2, colour_by="G1score", 
    size_by="G2Mscore") + fontsize + ggtitle("Before removal")

# After blocking on the phase score.
out2 <- plotReducedDim(sce.th2.block, use_dimred="PCA", ncomponents=2, 
    colour_by="G1score", size_by="G2Mscore") + fontsize + 
    ggtitle("After removal")
multiplot(out, out2, cols=2)