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.

library(scRNAseq)
sceZ <- ZeiselBrainData()

We compute log-expression values for use in marker detection inside SingleR().

library(scater)
sceZ <- logNormCounts(sceZ)

We examine the distribution of labels in this reference.

table(sceZ$level2class)
## 
##    (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        62

We 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:

table(pred.tasic$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        5

We can also examine the number of discarded cells for each label:

table(Label=pred.tasic$labels,
    Lost=is.na(pred.tasic$pruned.labels))
##           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    0

9.3 Diagnostics

We visualize the assignment scores for each label in Figure 9.1.

plotScoreHeatmap(pred.tasic)
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.

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.

plotDeltaDistribution(pred.tasic)
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.

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)
Heatmap of log-expression values in the Tasic dataset for all marker genes upregulated in the most frequent label from the Zeisel reference dataset.

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.

tab <- table(cluster=sceT$cluster, label=pred.tasic$labels) 
pheatmap::pheatmap(log10(tab+10))
Heatmap of the log-transformed number of cells in each combination of label (column) and cluster (row) in the Tasic dataset.

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))
$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.

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.4.0 beta (2024-04-15 r86425)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /home/biocbuild/bbs-3.19-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.14.0              scran_1.32.0               
 [3] SingleR_2.6.0               ensembldb_2.28.0           
 [5] AnnotationFilter_1.28.0     GenomicFeatures_1.56.0     
 [7] AnnotationDbi_1.66.0        AnnotationHub_3.12.0       
 [9] BiocFileCache_2.12.0        dbplyr_2.5.0               
[11] scater_1.32.0               ggplot2_3.5.1              
[13] scuttle_1.14.0              scRNAseq_2.17.10           
[15] SingleCellExperiment_1.26.0 SummarizedExperiment_1.34.0
[17] Biobase_2.64.0              GenomicRanges_1.56.0       
[19] GenomeInfoDb_1.40.0         IRanges_2.38.0             
[21] S4Vectors_0.42.0            BiocGenerics_0.50.0        
[23] MatrixGenerics_1.16.0       matrixStats_1.3.0          
[25] BiocStyle_2.32.0            rebook_1.14.0              

loaded via a namespace (and not attached):
  [1] RColorBrewer_1.1-3        jsonlite_1.8.8           
  [3] CodeDepends_0.6.6         magrittr_2.0.3           
  [5] ggbeeswarm_0.7.2          gypsum_1.0.0             
  [7] farver_2.1.1              rmarkdown_2.26           
  [9] BiocIO_1.14.0             zlibbioc_1.50.0          
 [11] vctrs_0.6.5               memoise_2.0.1            
 [13] Rsamtools_2.20.0          DelayedMatrixStats_1.26.0
 [15] RCurl_1.98-1.14           htmltools_0.5.8.1        
 [17] S4Arrays_1.4.0            curl_5.2.1               
 [19] BiocNeighbors_1.22.0      Rhdf5lib_1.26.0          
 [21] SparseArray_1.4.0         rhdf5_2.48.0             
 [23] sass_0.4.9                alabaster.base_1.4.0     
 [25] bslib_0.7.0               alabaster.sce_1.4.0      
 [27] httr2_1.0.1               cachem_1.0.8             
 [29] GenomicAlignments_1.40.0  igraph_2.0.3             
 [31] mime_0.12                 lifecycle_1.0.4          
 [33] pkgconfig_2.0.3           rsvd_1.0.5               
 [35] Matrix_1.7-0              R6_2.5.1                 
 [37] fastmap_1.1.1             GenomeInfoDbData_1.2.12  
 [39] digest_0.6.35             colorspace_2.1-0         
 [41] paws.storage_0.5.0        dqrng_0.3.2              
 [43] irlba_2.3.5.1             ExperimentHub_2.12.0     
 [45] RSQLite_2.3.6             beachmat_2.20.0          
 [47] labeling_0.4.3            filelock_1.0.3           
 [49] fansi_1.0.6               httr_1.4.7               
 [51] abind_1.4-5               compiler_4.4.0           
 [53] bit64_4.0.5               withr_3.0.0              
 [55] BiocParallel_1.38.0       viridis_0.6.5            
 [57] DBI_1.2.2                 highr_0.10               
 [59] HDF5Array_1.32.0          alabaster.ranges_1.4.0   
 [61] alabaster.schemas_1.4.0   rappdirs_0.3.3           
 [63] DelayedArray_0.30.0       rjson_0.2.21             
 [65] tools_4.4.0               vipor_0.4.7              
 [67] beeswarm_0.4.0            glue_1.7.0               
 [69] restfulr_0.0.15           rhdf5filters_1.16.0      
 [71] grid_4.4.0                Rtsne_0.17               
 [73] cluster_2.1.6             generics_0.1.3           
 [75] gtable_0.3.5              metapod_1.12.0           
 [77] BiocSingular_1.20.0       ScaledMatrix_1.12.0      
 [79] utf8_1.2.4                XVector_0.44.0           
 [81] ggrepel_0.9.5             BiocVersion_3.19.1       
 [83] pillar_1.9.0              limma_3.60.0             
 [85] dplyr_1.1.4               lattice_0.22-6           
 [87] rtracklayer_1.64.0        bit_4.0.5                
 [89] tidyselect_1.2.1          paws.common_0.7.2        
 [91] locfit_1.5-9.9            Biostrings_2.72.0        
 [93] knitr_1.46                gridExtra_2.3            
 [95] bookdown_0.39             ProtGenerics_1.36.0      
 [97] edgeR_4.2.0               xfun_0.43                
 [99] statmod_1.5.0             pheatmap_1.0.12          
[101] UCSC.utils_1.0.0          lazyeval_0.2.2           
[103] yaml_2.3.8                evaluate_0.23            
[105] codetools_0.2-20          tibble_3.2.1             
[107] alabaster.matrix_1.4.0    BiocManager_1.30.22      
[109] graph_1.82.0              cli_3.6.2                
[111] munsell_0.5.1             jquerylib_0.1.4          
[113] Rcpp_1.0.12               dir.expiry_1.12.0        
[115] png_0.1-8                 XML_3.99-0.16.1          
[117] parallel_4.4.0            blob_1.2.4               
[119] sparseMatrixStats_1.16.0  bitops_1.0-7             
[121] viridisLite_0.4.2         alabaster.se_1.4.0       
[123] scales_1.3.0              purrr_1.0.2              
[125] crayon_1.5.2              rlang_1.1.3              
[127] cowplot_1.1.3             KEGGREST_1.44.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.