Chapter 26 Zeisel mouse brain (STRT-Seq)

26.1 Introduction

Here, we examine a heterogeneous dataset from a study of cell types in the mouse brain (Zeisel et al. 2015). This contains approximately 3000 cells of varying types such as oligodendrocytes, microglia and neurons. Individual cells were isolated using the Fluidigm C1 microfluidics system (Pollen et al. 2014) and library preparation was performed on each cell using a UMI-based protocol. After sequencing, expression was quantified by counting the number of unique molecular identifiers (UMIs) mapped to each gene.

26.2 Data loading

We obtain a SingleCellExperiment object for this dataset using the relevant function from the scRNAseq package. The idiosyncrasies of the published dataset means that we need to do some extra work to merge together redundant rows corresponding to alternative genomic locations for the same gene.

We also fetch the Ensembl gene IDs, just in case we need them later.

26.3 Quality control

The original authors of the study have already removed low-quality cells prior to data publication. Nonetheless, we compute some quality control metrics to check whether the remaining cells are satisfactory.

Distribution of each QC metric across cells in the Zeisel brain dataset. Each point represents a cell and is colored according to whether that cell was discarded.

Figure 26.1: Distribution of each QC metric across cells in the Zeisel brain dataset. Each point represents a cell and is colored according to whether that cell was discarded.

Percentage of mitochondrial reads in each cell in the Zeisel brain dataset, compared to the total count (left) or the percentage of spike-in reads (right). Each point represents a cell and is colored according to whether that cell was discarded.

Figure 26.2: Percentage of mitochondrial reads in each cell in the Zeisel brain dataset, compared to the total count (left) or the percentage of spike-in reads (right). Each point represents a cell and is colored according to whether that cell was discarded.

We also examine the number of cells removed for each reason.

##              low_lib_size            low_n_features high_altexps_ERCC_percent 
##                         0                         3                        65 
##   high_subsets_Mt_percent                   discard 
##                       128                       189

26.4 Normalization

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.119   0.486   0.831   1.000   1.321   4.509
Relationship between the library size factors and the deconvolution size factors in the Zeisel brain dataset.

Figure 26.3: Relationship between the library size factors and the deconvolution size factors in the Zeisel brain dataset.

26.5 Variance modelling

In theory, we should block on the plate of origin for each cell. However, only 20-40 cells are available on each plate, and the population is also highly heterogeneous. This means that we cannot assume that the distribution of sampled cell types on each plate is the same. Thus, to avoid regressing out potential biology, we will not block on any factors in this analysis.

We see from Figure 26.4 that the technical and total variances are much smaller than those in the read-based datasets. This is due to the use of UMIs, which reduces the noise caused by variable PCR amplification. Furthermore, the spike-in trend is consistently lower than the variances of the endogenous gene, which reflects the heterogeneity in gene expression across cells of different types.

Per-gene variance as a function of the mean for the log-expression values in the Zeisel brain dataset. Each point represents a gene (black) with the mean-variance trend (blue) fitted to the spike-in transcripts (red).

Figure 26.4: Per-gene variance as a function of the mean for the log-expression values in the Zeisel brain dataset. Each point represents a gene (black) with the mean-variance trend (blue) fitted to the spike-in transcripts (red).

26.7 Clustering

## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14 
## 283 451 114 143 599 167 191 128 350  70 199  58  39  24
Obligatory $t$-SNE plot of the Zeisel brain dataset, where each point represents a cell and is colored according to the assigned cluster.

Figure 26.5: Obligatory \(t\)-SNE plot of the Zeisel brain dataset, where each point represents a cell and is colored according to the assigned cluster.

26.8 Interpretation

We focus on upregulated marker genes as these can quickly provide positive identification of cell type in a heterogeneous population. We examine the table for cluster 1, in which log-fold changes are reported between cluster 1 and every other cluster. The same output is provided for each cluster in order to identify genes that discriminate between clusters.

## DataFrame with 10 rows and 8 columns
##                Top      p.value          FDR summary.logFC   logFC.2   logFC.3
##          <integer>    <numeric>    <numeric>     <numeric> <numeric> <numeric>
## Atp1a3           1 1.45982e-282 7.24035e-279       3.45669 0.0398568 0.0893943
## Celf4            1 2.27030e-246 4.50404e-243       3.10465 0.3886716 0.6145023
## Gad1             1 7.44925e-232 1.34351e-228       4.57719 4.5392751 4.3003280
## Gad2             1 2.88086e-207 3.57208e-204       4.25393 4.2322487 3.8884654
## Mllt11           1 1.72982e-249 3.81309e-246       2.88363 0.5782719 1.4933128
## Ndrg4            1  0.00000e+00  0.00000e+00       3.84337 0.8887239 1.0183408
## Slc32a1          1 2.38276e-110 4.04030e-108       1.92859 1.9196173 1.8252062
## Syngr3           1 3.68257e-143 1.30462e-140       2.55531 1.0981258 1.1994793
## Atp6v1g2         2 3.04451e-204 3.55295e-201       2.50875 0.0981706 0.5203760
## Napb             2 1.10402e-231 1.82522e-228       2.81533 0.1774508 0.3046901
##            logFC.4   logFC.5
##          <numeric> <numeric>
## Atp1a3    1.241388   3.45669
## Celf4     0.869334   3.10465
## Gad1      4.050305   4.47236
## Gad2      3.769556   4.16902
## Mllt11    0.951649   2.88363
## Ndrg4     1.140041   3.84337
## Slc32a1   1.804311   1.92426
## Syngr3    1.188856   2.47696
## Atp6v1g2  0.616391   2.50875
## Napb      0.673772   2.81533

Figure 26.6 indicates that most of the top markers are strongly DE in cells of cluster 1 compared to some or all of the other clusters. We can use these markers to identify cells from cluster 1 in validation studies with an independent population of cells. A quick look at the markers suggest that cluster 1 represents interneurons based on expression of Gad1 and Slc6a1 (Zeng et al. 2012).

Heatmap of the log-expression of the top markers for cluster 1 compared to each other cluster. Cells are ordered by cluster and the color is scaled to the log-expression of each gene in each cell.

Figure 26.6: Heatmap of the log-expression of the top markers for cluster 1 compared to each other cluster. Cells are ordered by cluster and the color is scaled to the log-expression of each gene in each cell.

An alternative visualization approach is to plot the log-fold changes to all other clusters directly (Figure 26.7). This is more concise and is useful in situations involving many clusters that contain different numbers of cells.

Heatmap of the log-fold changes of the top markers for cluster 1 compared to each other cluster.

Figure 26.7: Heatmap of the log-fold changes of the top markers for cluster 1 compared to each other cluster.

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] pheatmap_1.0.12             BiocSingular_1.6.0         
 [3] scran_1.18.5                org.Mm.eg.db_3.12.0        
 [5] AnnotationDbi_1.52.0        scater_1.18.6              
 [7] ggplot2_3.3.3               scRNAseq_2.4.0             
 [9] SingleCellExperiment_1.12.0 SummarizedExperiment_1.20.0
[11] Biobase_2.50.0              GenomicRanges_1.42.0       
[13] GenomeInfoDb_1.26.4         IRanges_2.24.1             
[15] S4Vectors_0.28.1            BiocGenerics_0.36.0        
[17] MatrixGenerics_1.2.1        matrixStats_0.58.0         
[19] 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           scales_1.1.1                 
 [33] edgeR_3.32.1                  DBI_1.1.1                    
 [35] Rcpp_1.0.6                    viridisLite_0.3.0            
 [37] xtable_1.8-4                  progress_1.2.2               
 [39] dqrng_0.2.1                   bit_4.0.4                    
 [41] rsvd_1.0.3                    httr_1.4.2                   
 [43] RColorBrewer_1.1-2            ellipsis_0.3.1               
 [45] pkgconfig_2.0.3               XML_3.99-0.6                 
 [47] farver_2.1.0                  scuttle_1.0.4                
 [49] CodeDepends_0.6.5             sass_0.3.1                   
 [51] dbplyr_2.1.0                  locfit_1.5-9.4               
 [53] utf8_1.2.1                    tidyselect_1.1.0             
 [55] labeling_0.4.2                rlang_0.4.10                 
 [57] later_1.1.0.1                 munsell_0.5.0                
 [59] BiocVersion_3.12.0            tools_4.0.4                  
 [61] cachem_1.0.4                  generics_0.1.0               
 [63] RSQLite_2.2.4                 ExperimentHub_1.16.0         
 [65] evaluate_0.14                 stringr_1.4.0                
 [67] fastmap_1.1.0                 yaml_2.2.1                   
 [69] processx_3.4.5                knitr_1.31                   
 [71] bit64_4.0.5                   purrr_0.3.4                  
 [73] AnnotationFilter_1.14.0       sparseMatrixStats_1.2.1      
 [75] mime_0.10                     xml2_1.3.2                   
 [77] biomaRt_2.46.3                compiler_4.0.4               
 [79] beeswarm_0.3.1                curl_4.3                     
 [81] interactiveDisplayBase_1.28.0 statmod_1.4.35               
 [83] tibble_3.1.0                  bslib_0.2.4                  
 [85] stringi_1.5.3                 highr_0.8                    
 [87] ps_1.6.0                      GenomicFeatures_1.42.2       
 [89] lattice_0.20-41               bluster_1.0.0                
 [91] ProtGenerics_1.22.0           Matrix_1.3-2                 
 [93] vctrs_0.3.6                   pillar_1.5.1                 
 [95] lifecycle_1.0.0               BiocManager_1.30.10          
 [97] jquerylib_0.1.3               BiocNeighbors_1.8.2          
 [99] cowplot_1.1.1                 bitops_1.0-6                 
[101] irlba_2.3.3                   httpuv_1.5.5                 
[103] rtracklayer_1.50.0            R6_2.5.0                     
[105] bookdown_0.21                 promises_1.2.0.1             
[107] gridExtra_2.3                 vipor_0.4.5                  
[109] codetools_0.2-18              assertthat_0.2.1             
[111] openssl_1.4.3                 withr_2.4.1                  
[113] GenomicAlignments_1.26.0      Rsamtools_2.6.0              
[115] GenomeInfoDbData_1.2.4        hms_1.0.0                    
[117] grid_4.0.4                    beachmat_2.6.4               
[119] rmarkdown_2.7                 DelayedMatrixStats_1.12.3    
[121] Rtsne_0.15                    shiny_1.6.0                  
[123] ggbeeswarm_0.6.0

Bibliography

Pollen, A. A., T. J. Nowakowski, J. Shuga, X. Wang, A. A. Leyrat, J. H. Lui, N. Li, et al. 2014. “Low-coverage single-cell mRNA sequencing reveals cellular heterogeneity and activated signaling pathways in developing cerebral cortex.” Nat. Biotechnol. 32 (10): 1053–8.

Zeisel, A., A. B. Munoz-Manchado, S. Codeluppi, P. Lonnerberg, G. La Manno, A. Jureus, S. Marques, et al. 2015. “Brain structure. Cell types in the mouse cortex and hippocampus revealed by single-cell RNA-seq.” Science 347 (6226): 1138–42.

Zeng, H., E. H. Shen, J. G. Hohmann, S. W. Oh, A. Bernard, J. J. Royall, K. J. Glattfelder, et al. 2012. “Large-scale cellular-resolution gene profiling in human neocortex reveals species-specific molecular signatures.” Cell 149 (2): 483–96.