Chapter 5 Overview

5.1 Introduction

This chapter provides an overview of the framework of a typical scRNA-seq analysis workflow (Figure 5.1). Subsequent chapters will describe each analysis step in more detail.

Schematic of a typical scRNA-seq analysis workflow. Each stage (separated by dashed lines) consists of a number of specific steps, many of which operate on and modify a `SingleCellExperiment` instance.

Figure 5.1: Schematic of a typical scRNA-seq analysis workflow. Each stage (separated by dashed lines) consists of a number of specific steps, many of which operate on and modify a SingleCellExperiment instance.

5.2 Experimental Design

Before starting the analysis itself, some comments on experimental design may be helpful. The most obvious question is the choice of technology, which can be roughly divided into:

  • Droplet-based: 10X Genomics, inDrop, Drop-seq
  • Plate-based with unique molecular identifiers (UMIs): CEL-seq, MARS-seq
  • Plate-based with reads: Smart-seq2
  • Other: sci-RNA-seq, Seq-Well

Each of these methods have their own advantages and weaknesses that are discussed extensively elsewhere (Mereu et al. 2019; Ziegenhain et al. 2017). In practical terms, droplet-based technologies are the current de facto standard due to their throughput and low cost per cell. Plate-based methods can capture other phenotypic information (e.g., morphology) and are more amenable to customization. Read-based methods provide whole-transcript coverage, which is useful in some applications (e.g., splicing, exome mutations); otherwise, UMI-based methods are more popular as they mitigate the effects of PCR amplification noise. The choice of method is left to the reader’s circumstances - we will simply note that most aspects of our analysis pipeline are technology-agnostic.

The next question is how many cells should be captured, and to what depth they should be sequenced. The short answer is “as much as you can afford to spend”. The long answer is that it depends on the aim of the analysis. If we are aiming to discover rare cell subpopulations, then we need more cells. If we are aiming to characterize subtle differences, then we need more sequencing depth. As of time of writing, an informal survey of the literature suggests that typical droplet-based experiments would capture anywhere from 10,000 to 100,000 cells, sequenced at anywhere from 1,000 to 10,000 UMIs per cell (usually in inverse proportion to the number of cells). Droplet-based methods also have a trade-off between throughput and doublet rate that affects the true efficiency of sequencing.

For studies involving multiple samples or conditions, the design considerations are the same as those for bulk RNA-seq experiments. There should be multiple biological replicates for each condition and conditions should not be confounded with batch. Note that individual cells are not replicates; rather, we are referring to samples derived from replicate donors or cultures.

5.3 Obtaining a count matrix

Sequencing data from scRNA-seq experiments must be converted into a matrix of expression values that can be used for statistical analysis. Given the discrete nature of sequencing data, this is usually a count matrix containing the number of UMIs or reads mapped to each gene in each cell. The exact procedure for quantifying expression tends to be technology-dependent:

  • For 10X Genomics data, the CellRanger software suite provides a custom pipeline to obtain a count matrix. This uses STAR to align reads to the reference genome and then counts the number of unique UMIs mapped to each gene.
  • Pseudo-alignment methods such as alevin can be used to obtain a count matrix from the same data with greater efficiency. This avoids the need for explicit alignment, which reduces the compute time and memory usage.
  • For other highly multiplexed protocols, the scPipe package provides a more general pipeline for processing scRNA-seq data. This uses the Rsubread aligner to align reads and then counts UMIs per gene.
  • For CEL-seq or CEL-seq2 data, the scruff package provides a dedicated pipeline for quantification.
  • For read-based protocols, we can generally re-use the same pipelines for processing bulk RNA-seq data.
  • For any data involving spike-in transcripts, the spike-in sequences should be included as part of the reference genome during alignment and quantification.

After quantification, we import the count matrix into R and create a SingleCellExperiment object. This can be done with base methods (e.g., read.table()) followed by applying the SingleCellExperiment() constructor. Alternatively, for specific file formats, we can use dedicated methods from the DropletUtils (for 10X data) or tximport/tximeta packages (for pseudo-alignment methods). Depending on the origin of the data, this requires some vigilance:

  • Some feature-counting tools will report mapping statistics in the count matrix (e.g., the number of unaligned or unassigned reads). While these values can be useful for quality control, they would be misleading if treated as gene expression values. Thus, they should be removed (or at least moved to the colData) prior to further analyses.
  • Be careful of using the ^ERCC regular expression to detect spike-in rows in human data where the row names of the count matrix are gene symbols. An ERCC gene family actually exists in human annotation, so this would result in incorrect identification of genes as spike-in transcripts. This problem can be avoided by using count matrices with standard identifiers (e.g., Ensembl, Entrez).

5.4 Data processing and downstream analysis

In the simplest case, the workflow has the following form:

  1. We compute quality control metrics to remove low-quality cells that would interfere with downstream analyses. These cells may have been damaged during processing or may not have been fully captured by the sequencing protocol. Common metrics includes the total counts per cell, the proportion of spike-in or mitochondrial reads and the number of detected features.
  2. We convert the counts into normalized expression values to eliminate cell-specific biases (e.g., in capture efficiency). This allows us to perform explicit comparisons across cells in downstream steps like clustering. We also apply a transformation, typically log, to adjust for the mean-variance relationship.
  3. We perform feature selection to pick a subset of interesting features for downstream analysis. This is done by modelling the variance across cells for each gene and retaining genes that are highly variable. The aim is to reduce computational overhead and noise from uninteresting genes.
  4. We apply dimensionality reduction to compact the data and further reduce noise. Principal components analysis is typically used to obtain an initial low-rank representation for more computational work, followed by more aggressive methods like \(t\)-stochastic neighbor embedding for visualization purposes.
  5. We cluster cells into groups according to similarities in their (normalized) expression profiles. This aims to obtain groupings that serve as empirical proxies for distinct biological states. We typically interpret these groupings by identifying differentially expressed marker genes between clusters.

Additional steps such as data integration and cell annotation will be discussed in their respective chapters.

Session Info

R version 4.0.4 (2021-02-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /home/biocbuild/bbs-3.12-books/R/lib/libRblas.so
LAPACK: /home/biocbuild/bbs-3.12-books/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.5                scater_1.18.6              
 [3] ggplot2_3.3.3               scRNAseq_2.4.0             
 [5] SingleCellExperiment_1.12.0 SummarizedExperiment_1.20.0
 [7] Biobase_2.50.0              GenomicRanges_1.42.0       
 [9] GenomeInfoDb_1.26.4         IRanges_2.24.1             
[11] S4Vectors_0.28.1            BiocGenerics_0.36.0        
[13] MatrixGenerics_1.2.1        matrixStats_0.58.0         
[15] BiocStyle_2.18.1            rebook_1.0.0               

loaded via a namespace (and not attached):
  [1] AnnotationHub_2.22.0          BiocFileCache_1.14.0         
  [3] igraph_1.2.6                  lazyeval_0.2.2               
  [5] BiocParallel_1.24.1           digest_0.6.27                
  [7] ensembldb_2.14.0              htmltools_0.5.1.1            
  [9] viridis_0.5.1                 fansi_0.4.2                  
 [11] magrittr_2.0.1                memoise_2.0.0                
 [13] limma_3.46.0                  Biostrings_2.58.0            
 [15] askpass_1.1                   prettyunits_1.1.1            
 [17] colorspace_2.0-0              blob_1.2.1                   
 [19] rappdirs_0.3.3                xfun_0.22                    
 [21] dplyr_1.0.5                   callr_3.5.1                  
 [23] crayon_1.4.1                  RCurl_1.98-1.3               
 [25] jsonlite_1.7.2                graph_1.68.0                 
 [27] glue_1.4.2                    gtable_0.3.0                 
 [29] zlibbioc_1.36.0               XVector_0.30.0               
 [31] DelayedArray_0.16.2           BiocSingular_1.6.0           
 [33] scales_1.1.1                  DBI_1.1.1                    
 [35] edgeR_3.32.1                  Rcpp_1.0.6                   
 [37] viridisLite_0.3.0             xtable_1.8-4                 
 [39] progress_1.2.2                dqrng_0.2.1                  
 [41] bit_4.0.4                     rsvd_1.0.3                   
 [43] httr_1.4.2                    ellipsis_0.3.1               
 [45] farver_2.1.0                  pkgconfig_2.0.3              
 [47] XML_3.99-0.6                  scuttle_1.0.4                
 [49] uwot_0.1.10                   CodeDepends_0.6.5            
 [51] sass_0.3.1                    dbplyr_2.1.0                 
 [53] locfit_1.5-9.4                utf8_1.2.1                   
 [55] labeling_0.4.2                tidyselect_1.1.0             
 [57] rlang_0.4.10                  later_1.1.0.1                
 [59] AnnotationDbi_1.52.0          munsell_0.5.0                
 [61] BiocVersion_3.12.0            tools_4.0.4                  
 [63] cachem_1.0.4                  generics_0.1.0               
 [65] RSQLite_2.2.4                 ExperimentHub_1.16.0         
 [67] evaluate_0.14                 stringr_1.4.0                
 [69] fastmap_1.1.0                 yaml_2.2.1                   
 [71] processx_3.4.5                knitr_1.31                   
 [73] bit64_4.0.5                   purrr_0.3.4                  
 [75] AnnotationFilter_1.14.0       sparseMatrixStats_1.2.1      
 [77] mime_0.10                     xml2_1.3.2                   
 [79] biomaRt_2.46.3                compiler_4.0.4               
 [81] beeswarm_0.3.1                curl_4.3                     
 [83] interactiveDisplayBase_1.28.0 statmod_1.4.35               
 [85] tibble_3.1.0                  bslib_0.2.4                  
 [87] stringi_1.5.3                 highr_0.8                    
 [89] ps_1.6.0                      RSpectra_0.16-0              
 [91] GenomicFeatures_1.42.2        lattice_0.20-41              
 [93] bluster_1.0.0                 ProtGenerics_1.22.0          
 [95] Matrix_1.3-2                  vctrs_0.3.6                  
 [97] pillar_1.5.1                  lifecycle_1.0.0              
 [99] BiocManager_1.30.10           jquerylib_0.1.3              
[101] BiocNeighbors_1.8.2           cowplot_1.1.1                
[103] bitops_1.0-6                  irlba_2.3.3                  
[105] httpuv_1.5.5                  rtracklayer_1.50.0           
[107] R6_2.5.0                      bookdown_0.21                
[109] promises_1.2.0.1              gridExtra_2.3                
[111] vipor_0.4.5                   codetools_0.2-18             
[113] assertthat_0.2.1              openssl_1.4.3                
[115] withr_2.4.1                   GenomicAlignments_1.26.0     
[117] Rsamtools_2.6.0               GenomeInfoDbData_1.2.4       
[119] hms_1.0.0                     grid_4.0.4                   
[121] beachmat_2.6.4                rmarkdown_2.7                
[123] DelayedMatrixStats_1.12.3     shiny_1.6.0                  
[125] ggbeeswarm_0.6.0

Bibliography

Macosko, E. Z., A. Basu, R. Satija, J. Nemesh, K. Shekhar, M. Goldman, I. Tirosh, et al. 2015. “Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets.” Cell 161 (5): 1202–14.

Mereu, Elisabetta, Atefeh Lafzi, Catia Moutinho, Christoph Ziegenhain, Davis J. MacCarthy, Adrian Alvarez, Eduard Batlle, et al. 2019. “Benchmarking Single-Cell Rna Sequencing Protocols for Cell Atlas Projects.” bioRxiv. https://doi.org/10.1101/630087.

Ziegenhain, C., B. Vieth, S. Parekh, B. Reinius, A. Guillaumet-Adkins, M. Smets, H. Leonhardt, H. Heyn, I. Hellmann, and W. Enard. 2017. “Comparative Analysis of Single-Cell RNA Sequencing Methods.” Mol. Cell 65 (4): 631–43.