1 Overview

Here, we describe a few additional analyses that can be performed with single-cell RNA sequencing data. This includes detection of significant correlations between genes and regressing out the effect of cell cycle from the gene expression matrix.

2 Identifying correlated gene pairs with Spearman’s rho

scRNA-seq data is commonly used 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 haematopoietic stem cell (HSC) Smart-seq2 dataset (Wilson et al. 2015). Counts were obtained from NCBI GEO as a supplementary file using the accession number GSE61533, and are used to generate a SingleCellExperiment as shown below.

library(BiocFileCache)
bfc <- BiocFileCache("raw_data", ask=FALSE)
wilson.fname <- bfcrpath(bfc, file.path("ftp://ftp.ncbi.nlm.nih.gov/geo/series",
    "GSE61nnn/GSE61533/suppl/GSE61533_HTSEQ_count_results.xls.gz"))
library(R.utils)
wilson.name2 <- "GSE61533_HTSEQ_count_results.xls"
gunzip(wilson.fname, destname=wilson.name2, remove=FALSE, overwrite=TRUE)

library(readxl)
all.counts <- read_excel(wilson.name2)
gene.names <- all.counts$ID
all.counts <- as.matrix(all.counts[,-1])
rownames(all.counts) <- gene.names

library(SingleCellExperiment)
sce.hsc <- SingleCellExperiment(list(counts=all.counts))
is.spike <- grepl("^ERCC", rownames(sce.hsc))
sce.hsc <- splitAltExps(sce.hsc, ifelse(is.spike, "ERCC", "gene"))

library(scater)
sce.hsc <- addPerCellQC(sce.hsc)
spike.drop <- quickPerCellQC(colData(sce.hsc))
sce.hsc <- sce.hsc[,!spike.drop$discard]

library(scran)
sce.hsc <- computeSumFactors(sce.hsc)
sce.hsc <- logNormCounts(sce.hsc)

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)
## DataFrame with 6 rows and 6 columns
##         gene1       gene2       rho     p.value        FDR   limited
##   <character> <character> <numeric>   <numeric>  <numeric> <logical>
## 1       H2-Aa      H2-Ab1  0.529774 2.00000e-06 0.00189200      TRUE
## 2       H2-D1       H2-K1  0.436215 7.99999e-06 0.00378400     FALSE
## 3      H2-Ab1      H2-Eb1  0.422085 2.20000e-05 0.00693733     FALSE
## 4       H2-Aa      H2-Eb1  0.405762 4.40000e-05 0.01040599     FALSE
## 5       H2-Q6       H2-Q7  0.337576 9.43999e-04 0.17860462     FALSE
## 6       H2-K1       H2-K2  0.316598 2.05000e-03 0.32321634     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     942       4

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"))
## DataFrame with 1 row and 6 columns
##         gene1       gene2       rho     p.value         FDR   limited
##   <character> <character> <numeric>   <numeric>   <numeric> <logical>
## 1         Fos         Jun  0.449599 7.99999e-06 7.99999e-06     FALSE

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

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

Figure 1: Expression of Fos plotted against the expression of Jun for all cells in the HSC dataset.

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:

3 Comments on filtering by abundance

Low-abundance genes are problematic as zero or near-zero counts do not contain much information for reliable statistical inference. In applications involving hypothesis testing, these genes typically do not provide enough evidence to reject the null hypothesis yet they still increase the severity of the multiple testing correction. The discreteness of the counts may also interfere with statistical procedures, e.g., by compromising the accuracy of continuous approximations. Thus, low-abundance genes are often removed in many RNA-seq analysis pipelines before the application of downstream methods.

The “optimal” choice of filtering strategy depends on the downstream application. A more aggressive filter is usually required to remove discreteness and to avoid zeroes, e.g., for normalization purposes. By comparison, the filter statistic for hypothesis testing is mainly required to be independent of the test statistic under the null hypothesis (Bourgon, Gentleman, and Huber 2010). Given these differences in priorities, we (or the relevant function) will filter at each step as appropriate, rather than applying a single filter for the entire analysis. For example, computeSumFactors() will apply a somewhat stringent filter based on the average count, while fitTrendVar() will apply a relatively relaxed filter based on the average log-expression. Other applications will not do any abundance-based filtering at all (e.g., denoisePCA()) to preserve biological signal from lowly expressed genes.

Nonetheless, if global filtering is desired, it is simple to achieve by simply subsetting the SingleCellExperiment object. The example below demonstrates how we could remove genes with average counts less than 1 in the HSC dataset. The number of TRUE values in demo.keep corresponds to the number of retained rows/genes after filtering.

ave.counts <- calculateAverage(sce.hsc)
demo.keep <- ave.counts >= 1
filtered.sce.hsc <- sce.hsc[demo.keep,]
summary(demo.keep)
##    Mode   FALSE    TRUE 
## logical   24377   14029

4 Concluding remarks

All software packages used in this workflow are publicly available from the Comprehensive R Archive Network (https://cran.r-project.org) or the Bioconductor project (http://bioconductor.org). The specific version numbers of the packages used are shown below, along with the version of the R installation.

sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.12-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.12-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] scran_1.18.0                scater_1.18.0              
##  [3] ggplot2_3.3.2               SingleCellExperiment_1.12.0
##  [5] SummarizedExperiment_1.20.0 Biobase_2.50.0             
##  [7] GenomicRanges_1.42.0        GenomeInfoDb_1.26.0        
##  [9] IRanges_2.24.0              S4Vectors_0.28.0           
## [11] BiocGenerics_0.36.0         MatrixGenerics_1.2.0       
## [13] matrixStats_0.57.0          readxl_1.3.1               
## [15] R.utils_2.10.1              R.oo_1.24.0                
## [17] R.methodsS3_1.8.1           BiocFileCache_1.14.0       
## [19] dbplyr_1.4.4                knitr_1.30                 
## [21] BiocStyle_2.18.0           
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-6              bit64_4.0.5              
##  [3] httr_1.4.2                tools_4.0.3              
##  [5] R6_2.5.0                  irlba_2.3.3              
##  [7] vipor_0.4.5               DBI_1.1.0                
##  [9] colorspace_1.4-1          withr_2.3.0              
## [11] tidyselect_1.1.0          gridExtra_2.3            
## [13] bit_4.0.4                 curl_4.3                 
## [15] compiler_4.0.3            BiocNeighbors_1.8.0      
## [17] DelayedArray_0.16.0       labeling_0.4.2           
## [19] bookdown_0.21             scales_1.1.1             
## [21] rappdirs_0.3.1            stringr_1.4.0            
## [23] digest_0.6.27             rmarkdown_2.5            
## [25] XVector_0.30.0            pkgconfig_2.0.3          
## [27] htmltools_0.5.0           sparseMatrixStats_1.2.0  
## [29] highr_0.8                 limma_3.46.0             
## [31] rlang_0.4.8               RSQLite_2.2.1            
## [33] DelayedMatrixStats_1.12.0 farver_2.0.3             
## [35] generics_0.0.2            BiocParallel_1.24.0      
## [37] dplyr_1.0.2               RCurl_1.98-1.2           
## [39] magrittr_1.5              BiocSingular_1.6.0       
## [41] GenomeInfoDbData_1.2.4    scuttle_1.0.0            
## [43] Matrix_1.2-18             Rcpp_1.0.5               
## [45] ggbeeswarm_0.6.0          munsell_0.5.0            
## [47] viridis_0.5.1             lifecycle_0.2.0          
## [49] stringi_1.5.3             yaml_2.2.1               
## [51] edgeR_3.32.0              zlibbioc_1.36.0          
## [53] grid_4.0.3                blob_1.2.1               
## [55] dqrng_0.2.1               crayon_1.3.4             
## [57] lattice_0.20-41           cowplot_1.1.0            
## [59] beachmat_2.6.0            magick_2.5.0             
## [61] locfit_1.5-9.4            pillar_1.4.6             
## [63] igraph_1.2.6              codetools_0.2-16         
## [65] glue_1.4.2                evaluate_0.14            
## [67] BiocManager_1.30.10       vctrs_0.3.4              
## [69] cellranger_1.1.0          gtable_0.3.0             
## [71] purrr_0.3.4               assertthat_0.2.1         
## [73] xfun_0.18                 rsvd_1.0.3               
## [75] viridisLite_0.3.0         tibble_3.0.4             
## [77] beeswarm_0.2.3            memoise_1.1.0            
## [79] bluster_1.0.0             statmod_1.4.35           
## [81] ellipsis_0.3.1

References

Angel, P., and M. Karin. 1991. “The role of Jun, Fos and the AP-1 complex in cell-proliferation and transformation.” Biochim. Biophys. Acta 1072 (2-3): 129–57.

Bourgon, R., R. Gentleman, and W. Huber. 2010. “Independent filtering increases detection power for high-throughput experiments.” Proc. Natl. Acad. Sci. U.S.A. 107 (21): 9546–51.

Phipson, B., and G. K. Smyth. 2010. “Permutation P-values should never be zero: calculating exact P-values when permutations are randomly drawn.” Stat. Appl. Genet. Mol. Biol. 9: Article 39.

Simes, R. J. 1986. “An Improved Bonferroni Procedure for Multiple Tests of Significance.” Biometrika 73 (3): 751–54.

Wilson, N. K., D. G. Kent, F. Buettner, M. Shehata, I. C. Macaulay, F. J. Calero-Nieto, M. Sanchez Castillo, et al. 2015. “Combined single-cell functional and gene expression analysis resolves heterogeneity within stem cell populations.” Cell Stem Cell 16 (6): 712–24.