Chapter 2 Zeisel mouse brain (STRT-Seq)
2.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.
2.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.
library(scRNAseq)
sce.zeisel <- ZeiselBrainData()
library(scater)
sce.zeisel <- aggregateAcrossFeatures(sce.zeisel,
id=sub("_loc[0-9]+$", "", rownames(sce.zeisel)))
We also fetch the Ensembl gene IDs, just in case we need them later.
2.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.
stats <- perCellQCMetrics(sce.zeisel, subsets=list(
Mt=rowData(sce.zeisel)$featureType=="mito"))
qc <- quickPerCellQC(stats, percent_subsets=c("altexps_ERCC_percent",
"subsets_Mt_percent"))
sce.zeisel <- sce.zeisel[,!qc$discard]
colData(unfiltered) <- cbind(colData(unfiltered), stats)
unfiltered$discard <- qc$discard
gridExtra::grid.arrange(
plotColData(unfiltered, y="sum", colour_by="discard") +
scale_y_log10() + ggtitle("Total count"),
plotColData(unfiltered, y="detected", colour_by="discard") +
scale_y_log10() + ggtitle("Detected features"),
plotColData(unfiltered, y="altexps_ERCC_percent",
colour_by="discard") + ggtitle("ERCC percent"),
plotColData(unfiltered, y="subsets_Mt_percent",
colour_by="discard") + ggtitle("Mito percent"),
ncol=2
)
gridExtra::grid.arrange(
plotColData(unfiltered, x="sum", y="subsets_Mt_percent",
colour_by="discard") + scale_x_log10(),
plotColData(unfiltered, x="altexps_ERCC_percent", y="subsets_Mt_percent",
colour_by="discard"),
ncol=2
)
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
2.4 Normalization
library(scran)
set.seed(1000)
clusters <- quickCluster(sce.zeisel)
sce.zeisel <- computeSumFactors(sce.zeisel, cluster=clusters)
sce.zeisel <- logNormCounts(sce.zeisel)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.119 0.486 0.831 1.000 1.321 4.509
plot(librarySizeFactors(sce.zeisel), sizeFactors(sce.zeisel), pch=16,
xlab="Library size factors", ylab="Deconvolution factors", log="xy")
2.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.
dec.zeisel <- modelGeneVarWithSpikes(sce.zeisel, "ERCC")
top.hvgs <- getTopHVGs(dec.zeisel, prop=0.1)
We see from Figure 2.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.
plot(dec.zeisel$mean, dec.zeisel$total, pch=16, cex=0.5,
xlab="Mean of log-expression", ylab="Variance of log-expression")
curfit <- metadata(dec.zeisel)
points(curfit$mean, curfit$var, col="red", pch=16)
curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
2.6 Dimensionality reduction
library(BiocSingular)
set.seed(101011001)
sce.zeisel <- denoisePCA(sce.zeisel, technical=dec.zeisel, subset.row=top.hvgs)
sce.zeisel <- runTSNE(sce.zeisel, dimred="PCA")
We have a look at the number of PCs retained by denoisePCA()
.
## [1] 50
2.7 Clustering
snn.gr <- buildSNNGraph(sce.zeisel, use.dimred="PCA")
colLabels(sce.zeisel) <- factor(igraph::cluster_walktrap(snn.gr)$membership)
##
## 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
2.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.
markers <- findMarkers(sce.zeisel, direction="up")
marker.set <- markers[["1"]]
head(marker.set[,1:8], 10) # only first 8 columns, for brevity
## 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 2.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).
top.markers <- rownames(marker.set)[marker.set$Top <= 10]
plotHeatmap(sce.zeisel, features=top.markers, order_columns_by="label")
An alternative visualization approach is to plot the log-fold changes to all other clusters directly (Figure 2.7). This is more concise and is useful in situations involving many clusters that contain different numbers of cells.
library(pheatmap)
logFCs <- getMarkerEffects(marker.set[1:50,])
pheatmap(logFCs, breaks=seq(-5, 5, length.out=101))
Session Info
R version 4.1.1 (2021-08-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.3 LTS
Matrix products: default
BLAS: /home/biocbuild/bbs-3.14-bioc/R/lib/libRblas.so
LAPACK: /home/biocbuild/bbs-3.14-bioc/R/lib/libRlapack.so
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_GB 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] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] pheatmap_1.0.12 BiocSingular_1.10.0
[3] scran_1.22.0 org.Mm.eg.db_3.14.0
[5] AnnotationDbi_1.56.0 scater_1.22.0
[7] ggplot2_3.3.5 scuttle_1.4.0
[9] scRNAseq_2.7.2 SingleCellExperiment_1.16.0
[11] SummarizedExperiment_1.24.0 Biobase_2.54.0
[13] GenomicRanges_1.46.0 GenomeInfoDb_1.30.0
[15] IRanges_2.28.0 S4Vectors_0.32.0
[17] BiocGenerics_0.40.0 MatrixGenerics_1.6.0
[19] matrixStats_0.61.0 BiocStyle_2.22.0
[21] rebook_1.4.0
loaded via a namespace (and not attached):
[1] AnnotationHub_3.2.0 BiocFileCache_2.2.0
[3] igraph_1.2.7 lazyeval_0.2.2
[5] BiocParallel_1.28.0 digest_0.6.28
[7] ensembldb_2.18.0 htmltools_0.5.2
[9] viridis_0.6.2 fansi_0.5.0
[11] magrittr_2.0.1 memoise_2.0.0
[13] ScaledMatrix_1.2.0 cluster_2.1.2
[15] limma_3.50.0 Biostrings_2.62.0
[17] prettyunits_1.1.1 colorspace_2.0-2
[19] blob_1.2.2 rappdirs_0.3.3
[21] ggrepel_0.9.1 xfun_0.27
[23] dplyr_1.0.7 crayon_1.4.1
[25] RCurl_1.98-1.5 jsonlite_1.7.2
[27] graph_1.72.0 glue_1.4.2
[29] gtable_0.3.0 zlibbioc_1.40.0
[31] XVector_0.34.0 DelayedArray_0.20.0
[33] scales_1.1.1 edgeR_3.36.0
[35] DBI_1.1.1 Rcpp_1.0.7
[37] viridisLite_0.4.0 xtable_1.8-4
[39] progress_1.2.2 dqrng_0.3.0
[41] bit_4.0.4 rsvd_1.0.5
[43] metapod_1.2.0 httr_1.4.2
[45] RColorBrewer_1.1-2 dir.expiry_1.2.0
[47] ellipsis_0.3.2 pkgconfig_2.0.3
[49] XML_3.99-0.8 farver_2.1.0
[51] CodeDepends_0.6.5 sass_0.4.0
[53] dbplyr_2.1.1 locfit_1.5-9.4
[55] utf8_1.2.2 labeling_0.4.2
[57] tidyselect_1.1.1 rlang_0.4.12
[59] later_1.3.0 munsell_0.5.0
[61] BiocVersion_3.14.0 tools_4.1.1
[63] cachem_1.0.6 generics_0.1.1
[65] RSQLite_2.2.8 ExperimentHub_2.2.0
[67] evaluate_0.14 stringr_1.4.0
[69] fastmap_1.1.0 yaml_2.2.1
[71] knitr_1.36 bit64_4.0.5
[73] purrr_0.3.4 KEGGREST_1.34.0
[75] AnnotationFilter_1.18.0 sparseMatrixStats_1.6.0
[77] mime_0.12 xml2_1.3.2
[79] biomaRt_2.50.0 compiler_4.1.1
[81] beeswarm_0.4.0 filelock_1.0.2
[83] curl_4.3.2 png_0.1-7
[85] interactiveDisplayBase_1.32.0 statmod_1.4.36
[87] tibble_3.1.5 bslib_0.3.1
[89] stringi_1.7.5 highr_0.9
[91] GenomicFeatures_1.46.0 bluster_1.4.0
[93] lattice_0.20-45 ProtGenerics_1.26.0
[95] Matrix_1.3-4 vctrs_0.3.8
[97] pillar_1.6.4 lifecycle_1.0.1
[99] BiocManager_1.30.16 jquerylib_0.1.4
[101] BiocNeighbors_1.12.0 cowplot_1.1.1
[103] bitops_1.0-7 irlba_2.3.3
[105] httpuv_1.6.3 rtracklayer_1.54.0
[107] R6_2.5.1 BiocIO_1.4.0
[109] bookdown_0.24 promises_1.2.0.1
[111] gridExtra_2.3 vipor_0.4.5
[113] codetools_0.2-18 assertthat_0.2.1
[115] rjson_0.2.20 withr_2.4.2
[117] GenomicAlignments_1.30.0 Rsamtools_2.10.0
[119] GenomeInfoDbData_1.2.7 parallel_4.1.1
[121] hms_1.1.1 grid_4.1.1
[123] beachmat_2.10.0 rmarkdown_2.11
[125] DelayedMatrixStats_1.16.0 Rtsne_0.15
[127] shiny_1.7.1 ggbeeswarm_0.6.0
[129] restfulr_0.0.13
References
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.