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 described in the previous workflows. The second class uses the fact that the same amount of spike-in RNA was added to each cell (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.

2 Setting up the data

2.1 Obtaining the dataset

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.

library(BiocFileCache)
bfc <- BiocFileCache("raw_data", ask=FALSE)
islam.fname <- bfcrpath(bfc, file.path("ftp://ftp.ncbi.nlm.nih.gov/geo/series",
    "GSE29nnn/GSE29087/suppl/GSE29087_L139_expression_tab.txt.gz"))

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(islam.fname,
    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

2.2 Applying quality control

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)
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_by_counts, 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

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)
head(sizeFactors(sce.islam))    
## [1] 1.1486524 1.1274936 0.4285218 1.1011014 0.6646450 1.5932121
head(sizeFactors(sce.islam, "spike")) # same as general size factors.
## [1] 1.1486524 1.1274936 0.4285218 1.1011014 0.6646450 1.5932121

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 for this data set (Lun, Bach, and Marioni 2016).

deconv.sf <- computeSumFactors(sce.islam, sf.out=TRUE, 
    cluster=sce.islam$grouping)
head(deconv.sf)
## [1] 0.06424167 0.12172399 0.15547625 0.12876863 0.05680651 0.07888560

We observe a negative correlation between the two sets of size factors (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")
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:

  • 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 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 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.2 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.10-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.10-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] org.Mm.eg.db_3.8.2                    
##  [2] readxl_1.3.1                          
##  [3] scRNAseq_1.11.2                       
##  [4] TxDb.Mmusculus.UCSC.mm10.ensGene_3.4.0
##  [5] GenomicFeatures_1.37.1                
##  [6] Matrix_1.2-17                         
##  [7] edgeR_3.27.4                          
##  [8] limma_3.41.4                          
##  [9] BiocNeighbors_1.3.2                   
## [10] TENxBrainData_1.5.0                   
## [11] HDF5Array_1.13.2                      
## [12] rhdf5_2.29.0                          
## [13] batchelor_1.1.4                       
## [14] BiocSingular_1.1.3                    
## [15] scran_1.13.7                          
## [16] scater_1.13.7                         
## [17] ggplot2_3.2.0                         
## [18] SingleCellExperiment_1.7.0            
## [19] SummarizedExperiment_1.15.3           
## [20] DelayedArray_0.11.1                   
## [21] BiocParallel_1.19.0                   
## [22] matrixStats_0.54.0                    
## [23] GenomicRanges_1.37.11                 
## [24] GenomeInfoDb_1.21.1                   
## [25] org.Hs.eg.db_3.8.2                    
## [26] AnnotationDbi_1.47.0                  
## [27] IRanges_2.19.10                       
## [28] S4Vectors_0.23.12                     
## [29] Biobase_2.45.0                        
## [30] BiocGenerics_0.31.4                   
## [31] BiocFileCache_1.9.1                   
## [32] dbplyr_1.4.1                          
## [33] knitr_1.23                            
## [34] BiocStyle_2.13.2                      
## 
## loaded via a namespace (and not attached):
##   [1] AnnotationHub_2.17.3          RcppEigen_0.3.3.5.0          
##   [3] igraph_1.2.4.1                lazyeval_0.2.2               
##   [5] sp_1.3-1                      splines_3.6.0                
##   [7] digest_0.6.19                 htmltools_0.3.6              
##   [9] viridis_0.5.1                 magrittr_1.5                 
##  [11] memoise_1.1.0                 openxlsx_4.1.0.1             
##  [13] Biostrings_2.53.0             xts_0.11-2                   
##  [15] prettyunits_1.0.2             colorspace_1.4-1             
##  [17] blob_1.1.1                    rappdirs_0.3.1               
##  [19] haven_2.1.0                   xfun_0.7                     
##  [21] dplyr_0.8.1                   callr_3.2.0                  
##  [23] crayon_1.3.4                  RCurl_1.95-4.12              
##  [25] zoo_1.8-6                     glue_1.3.1                   
##  [27] gtable_0.3.0                  zlibbioc_1.31.0              
##  [29] XVector_0.25.0                car_3.0-3                    
##  [31] Rhdf5lib_1.7.2                DEoptimR_1.0-8               
##  [33] abind_1.4-5                   VIM_4.8.0                    
##  [35] scales_1.0.0                  pheatmap_1.0.12              
##  [37] DBI_1.0.0                     ggthemes_4.2.0               
##  [39] Rcpp_1.0.1                    laeken_0.5.0                 
##  [41] viridisLite_0.3.0             xtable_1.8-4                 
##  [43] progress_1.2.2                dqrng_0.2.1                  
##  [45] proxy_0.4-23                  foreign_0.8-71               
##  [47] bit_1.1-14                    rsvd_1.0.1                   
##  [49] vcd_1.4-4                     httr_1.4.0                   
##  [51] RColorBrewer_1.1-2            pkgconfig_2.0.2              
##  [53] XML_3.98-1.20                 nnet_7.3-12                  
##  [55] locfit_1.5-9.1                dynamicTreeCut_1.63-1        
##  [57] tidyselect_0.2.5              labeling_0.3                 
##  [59] rlang_0.3.4                   later_0.8.0                  
##  [61] munsell_0.5.0                 cellranger_1.1.0             
##  [63] tools_3.6.0                   RSQLite_2.1.1                
##  [65] ExperimentHub_1.11.1          ranger_0.11.2                
##  [67] evaluate_0.14                 stringr_1.4.0                
##  [69] yaml_2.2.0                    processx_3.3.1               
##  [71] bit64_0.9-7                   zip_2.0.2                    
##  [73] robustbase_0.93-5             purrr_0.3.2                  
##  [75] mime_0.7                      biomaRt_2.41.3               
##  [77] compiler_3.6.0                beeswarm_0.2.3               
##  [79] curl_3.3                      interactiveDisplayBase_1.23.0
##  [81] e1071_1.7-2                   smoother_1.1                 
##  [83] tibble_2.1.3                  statmod_1.4.32               
##  [85] stringi_1.4.3                 highr_0.8                    
##  [87] ps_1.3.0                      forcats_0.4.0                
##  [89] lattice_0.20-38               pillar_1.4.1                 
##  [91] simpleSingleCell_1.9.8        BiocManager_1.30.4           
##  [93] lmtest_0.9-37                 data.table_1.12.2            
##  [95] cowplot_0.9.4                 bitops_1.0-6                 
##  [97] irlba_2.3.3                   httpuv_1.5.1                 
##  [99] rtracklayer_1.45.1            R6_2.4.0                     
## [101] bookdown_0.11                 promises_1.0.1               
## [103] rio_0.5.16                    gridExtra_2.3                
## [105] vipor_0.4.5                   codetools_0.2-16             
## [107] boot_1.3-22                   MASS_7.3-51.4                
## [109] assertthat_0.2.1              destiny_2.15.0               
## [111] withr_2.1.2                   GenomicAlignments_1.21.2     
## [113] Rsamtools_2.1.2               GenomeInfoDbData_1.2.1       
## [115] hms_0.4.2                     grid_3.6.0                   
## [117] beachmat_2.1.1                class_7.3-15                 
## [119] rmarkdown_1.13                DelayedMatrixStats_1.7.1     
## [121] carData_3.0-2                 Rtsne_0.15                   
## [123] TTR_0.23-4                    scatterplot3d_0.3-41         
## [125] shiny_1.3.2                   ggbeeswarm_0.6.0

References

Buettner, F., K. N. Natarajan, F. P. Casale, V. Proserpio, A. Scialdone, F. J. Theis, S. A. Teichmann, J. C. Marioni, and O. Stegle. 2015. “Computational analysis of cell-to-cell heterogeneity in single-cell RNA-sequencing data reveals hidden subpopulations of cells.” Nat. Biotechnol. 33 (2):155–60.

Islam, S., U. Kjallquist, A. Moliner, P. Zajac, J. B. Fan, P. Lonnerberg, and S. Linnarsson. 2011. “Characterization of the single-cell transcriptional landscape by highly multiplex RNA-seq.” Genome Res. 21 (7):1160–7.

Lun, A. T. L., F. J. Calero-Nieto, L. Haim-Vilmovsky, B. Gottgens, and J. C. Marioni. 2017. “Assessing the reliability of spike-in normalization for analyses of single-cell RNA sequencing data.” Genome Res. 27 (11):1795–1806.

Lun, A. T., K. Bach, and J. C. Marioni. 2016. “Pooling across cells to normalize single-cell RNA sequencing data with many zero counts.” Genome Biol. 17 (April):75.