Chapter 9 Cross-annotating mouse brains
9.1 Loading the data
We load the classic Zeisel et al. (2015) dataset as our reference. Here, we’ll rely on the fact that the authors have already performed quality control.
We compute log-expression values for use in marker detection inside SingleR().
We examine the distribution of labels in this reference.
##
## (none) Astro1 Astro2 CA1Pyr1 CA1Pyr2 CA1PyrInt CA2Pyr2 Choroid
## 189 68 61 380 447 49 41 10
## ClauPyr Epend Int1 Int10 Int11 Int12 Int13 Int14
## 5 20 12 21 10 21 15 22
## Int15 Int16 Int2 Int3 Int4 Int5 Int6 Int7
## 18 20 24 10 15 20 22 23
## Int8 Int9 Mgl1 Mgl2 Oligo1 Oligo2 Oligo3 Oligo4
## 26 11 17 16 45 98 87 106
## Oligo5 Oligo6 Peric Pvm1 Pvm2 S1PyrDL S1PyrL23 S1PyrL4
## 125 359 21 32 33 81 74 26
## S1PyrL5 S1PyrL5a S1PyrL6 S1PyrL6b SubPyr Vend1 Vend2 Vsmc
## 16 28 39 21 22 32 105 62We load the Tasic et al. (2016) dataset as our test. While not strictly necessary, we remove putative low-quality cells to simplify later interpretation.
sceT <- TasicBrainData()
sceT <- addPerCellQC(sceT, subsets=list(mito=grep("^mt_", rownames(sceT))))
qc <- quickPerCellQC(colData(sceT),
percent_subsets=c("subsets_mito_percent", "altexps_ERCC_percent"))
sceT <- sceT[,which(!qc$discard)]The Tasic dataset was generated using read-based technologies so we need to adjust for the transcript length.
library(AnnotationHub)
mm.db <- AnnotationHub()[["AH73905"]]
mm.exons <- exonsBy(mm.db, by="gene")
mm.exons <- reduce(mm.exons)
mm.len <- sum(width(mm.exons))
mm.symb <- mapIds(mm.db, keys=names(mm.len), keytype="GENEID", column="SYMBOL")
names(mm.len) <- mm.symb
library(scater)
keep <- intersect(names(mm.len), rownames(sceT))
sceT <- sceT[keep,]
assay(sceT, "TPM") <- calculateTPM(sceT, lengths=mm.len[keep])9.2 Applying the annotation
We apply SingleR() with Wilcoxon rank sum test-based marker detection to annotate the Tasic dataset with the Zeisel labels.
library(SingleR)
pred.tasic <- SingleR(test=sceT, ref=sceZ, labels=sceZ$level2class,
assay.type.test="TPM", de.method="wilcox")We examine the distribution of predicted labels:
##
## Astro1 Astro2 CA1Pyr2 CA2Pyr2 ClauPyr Int1 Int10 Int11
## 1 5 1 3 1 152 98 2
## Int12 Int13 Int14 Int15 Int16 Int2 Int3 Int4
## 9 18 24 16 10 146 94 29
## Int6 Int7 Int8 Int9 Oligo1 Oligo2 Oligo3 Oligo4
## 14 2 35 31 8 1 7 1
## Oligo6 Peric S1PyrDL S1PyrL23 S1PyrL4 S1PyrL5 S1PyrL5a S1PyrL6
## 1 1 320 73 16 4 201 46
## S1PyrL6b SubPyr
## 72 5We can also examine the number of discarded cells for each label:
## Lost
## Label FALSE TRUE
## Astro1 1 0
## Astro2 5 0
## CA1Pyr2 1 0
## CA2Pyr2 3 0
## ClauPyr 1 0
## Int1 152 0
## Int10 98 0
## Int11 2 0
## Int12 9 0
## Int13 18 0
## Int14 23 1
## Int15 16 0
## Int16 10 0
## Int2 144 2
## Int3 94 0
## Int4 29 0
## Int6 14 0
## Int7 2 0
## Int8 35 0
## Int9 31 0
## Oligo1 8 0
## Oligo2 1 0
## Oligo3 7 0
## Oligo4 1 0
## Oligo6 1 0
## Peric 1 0
## S1PyrDL 304 16
## S1PyrL23 73 0
## S1PyrL4 16 0
## S1PyrL5 4 0
## S1PyrL5a 200 1
## S1PyrL6 45 1
## S1PyrL6b 72 0
## SubPyr 5 09.3 Diagnostics
We visualize the assignment scores for each label in Figure 9.1.
Figure 9.1: Heatmap of the (normalized) assignment scores for each cell (column) in the Tasic test dataset with respect to each label (row) in the Zeisel reference dataset. The final assignment for each cell is shown in the annotation bar at the top.
The delta for each cell is visualized in Figure 9.2.
Figure 9.2: Distributions of the deltas for each cell in the Tasic dataset assigned to each label in the Zeisel dataset. Each cell is represented by a point; low-quality assignments that were pruned out are colored in orange.
Finally, we visualize the heatmaps of the marker genes for the most frequent label in Figure 9.3. We could show these for all labels but I wouldn’t want to bore you with a parade of large heatmaps.
library(scater)
collected <- list()
all.markers <- metadata(pred.tasic)$de.genes
sceT <- logNormCounts(sceT)
top.label <- names(sort(table(pred.tasic$labels), decreasing=TRUE))[1]
per.label <- sumCountsAcrossCells(logcounts(sceT),
ids=pred.tasic$labels, average=TRUE)
per.label <- assay(per.label)[unique(unlist(all.markers[[top.label]])),]
pheatmap::pheatmap(per.label, main=top.label)
Figure 9.3: Heatmap of log-expression values in the Tasic dataset for all marker genes upregulated in the most frequent label from the Zeisel reference dataset.
9.4 Comparison to clusters
For comparison, we will perform a quick unsupervised analysis of the Grun dataset. We model the variances using the spike-in data and we perform graph-based clustering.
library(scran)
decT <- modelGeneVarWithSpikes(sceT, "ERCC")
set.seed(1000100)
sceT <- denoisePCA(sceT, decT, subset.row=getTopHVGs(decT, n=2500))
library(bluster)
sceT$cluster <- clusterRows(reducedDim(sceT, "PCA"), NNGraphParam())We do not observe a clean 1:1 mapping between clusters and labels in Figure 9.4, probably because many of the labels represent closely related cell types that are difficult to distinguish.
Figure 9.4: Heatmap of the log-transformed number of cells in each combination of label (column) and cluster (row) in the Tasic dataset.
We proceed to the most important part of the analysis. Yes, that’s right, the \(t\)-SNE plot (Figure 9.5).
set.seed(101010100)
sceT <- runTSNE(sceT, dimred="PCA")
plotTSNE(sceT, colour_by="cluster", text_colour="red",
text_by=I(pred.tasic$labels))
Figure 9.5: \(t\)-SNE plot of the Tasic dataset, where each point is a cell and is colored by the assigned cluster. Reference labels from the Zeisel dataset are also placed on the median coordinate across all cells assigned with that label.
Session information
R version 4.3.2 Patched (2023-11-13 r85521)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS
Matrix products: default
BLAS: /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
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
time zone: America/New_York
tzcode source: system (glibc)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] bluster_1.12.0 scran_1.30.0
[3] SingleR_2.4.0 ensembldb_2.26.0
[5] AnnotationFilter_1.26.0 GenomicFeatures_1.54.1
[7] AnnotationDbi_1.64.1 AnnotationHub_3.10.0
[9] BiocFileCache_2.10.1 dbplyr_2.4.0
[11] scater_1.30.1 ggplot2_3.4.4
[13] scuttle_1.12.0 scRNAseq_2.16.0
[15] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
[17] Biobase_2.62.0 GenomicRanges_1.54.1
[19] GenomeInfoDb_1.38.1 IRanges_2.36.0
[21] S4Vectors_0.40.2 BiocGenerics_0.48.1
[23] MatrixGenerics_1.14.0 matrixStats_1.1.0
[25] BiocStyle_2.30.0 rebook_1.12.0
loaded via a namespace (and not attached):
[1] RColorBrewer_1.1-3 jsonlite_1.8.8
[3] CodeDepends_0.6.5 magrittr_2.0.3
[5] ggbeeswarm_0.7.2 farver_2.1.1
[7] rmarkdown_2.25 BiocIO_1.12.0
[9] zlibbioc_1.48.0 vctrs_0.6.5
[11] memoise_2.0.1 Rsamtools_2.18.0
[13] DelayedMatrixStats_1.24.0 RCurl_1.98-1.13
[15] htmltools_0.5.7 S4Arrays_1.2.0
[17] progress_1.2.2 curl_5.1.0
[19] BiocNeighbors_1.20.0 SparseArray_1.2.2
[21] sass_0.4.7 bslib_0.6.1
[23] cachem_1.0.8 GenomicAlignments_1.38.0
[25] igraph_1.5.1 mime_0.12
[27] lifecycle_1.0.4 pkgconfig_2.0.3
[29] rsvd_1.0.5 Matrix_1.6-4
[31] R6_2.5.1 fastmap_1.1.1
[33] GenomeInfoDbData_1.2.11 shiny_1.8.0
[35] digest_0.6.33 colorspace_2.1-0
[37] dqrng_0.3.2 irlba_2.3.5.1
[39] ExperimentHub_2.10.0 RSQLite_2.3.3
[41] beachmat_2.18.0 labeling_0.4.3
[43] filelock_1.0.2 fansi_1.0.5
[45] httr_1.4.7 abind_1.4-5
[47] compiler_4.3.2 bit64_4.0.5
[49] withr_2.5.2 BiocParallel_1.36.0
[51] viridis_0.6.4 DBI_1.1.3
[53] highr_0.10 biomaRt_2.58.0
[55] rappdirs_0.3.3 DelayedArray_0.28.0
[57] rjson_0.2.21 tools_4.3.2
[59] vipor_0.4.5 beeswarm_0.4.0
[61] interactiveDisplayBase_1.40.0 httpuv_1.6.12
[63] glue_1.6.2 restfulr_0.0.15
[65] promises_1.2.1 grid_4.3.2
[67] Rtsne_0.16 cluster_2.1.6
[69] generics_0.1.3 gtable_0.3.4
[71] hms_1.1.3 metapod_1.10.0
[73] BiocSingular_1.18.0 ScaledMatrix_1.10.0
[75] xml2_1.3.6 utf8_1.2.4
[77] XVector_0.42.0 ggrepel_0.9.4
[79] BiocVersion_3.18.1 pillar_1.9.0
[81] stringr_1.5.1 limma_3.58.1
[83] later_1.3.1 dplyr_1.1.4
[85] lattice_0.22-5 rtracklayer_1.62.0
[87] bit_4.0.5 tidyselect_1.2.0
[89] locfit_1.5-9.8 Biostrings_2.70.1
[91] knitr_1.45 gridExtra_2.3
[93] bookdown_0.37 ProtGenerics_1.34.0
[95] edgeR_4.0.2 xfun_0.41
[97] statmod_1.5.0 pheatmap_1.0.12
[99] stringi_1.8.2 lazyeval_0.2.2
[101] yaml_2.3.7 evaluate_0.23
[103] codetools_0.2-19 tibble_3.2.1
[105] BiocManager_1.30.22 graph_1.80.0
[107] cli_3.6.1 xtable_1.8-4
[109] munsell_0.5.0 jquerylib_0.1.4
[111] Rcpp_1.0.11 dir.expiry_1.10.0
[113] png_0.1-8 XML_3.99-0.16
[115] parallel_4.3.2 ellipsis_0.3.2
[117] blob_1.2.4 prettyunits_1.2.0
[119] sparseMatrixStats_1.14.0 bitops_1.0-7
[121] viridisLite_0.4.2 scales_1.3.0
[123] purrr_1.0.2 crayon_1.5.2
[125] rlang_1.1.2 cowplot_1.1.1
[127] KEGGREST_1.42.0 Bibliography
Tasic, B., V. Menon, T. N. Nguyen, T. K. Kim, T. Jarsky, Z. Yao, B. Levi, et al. 2016. “Adult mouse cortical cell taxonomy revealed by single cell transcriptomics.” Nat. Neurosci. 19 (2): 335–46.
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.