Chapter 8 Cross-annotating human pancreas

8.1 Loading the data

We load the Muraro et al. (2016) dataset as our reference, removing unlabelled cells or cells without a clear label.

library(scRNAseq)
sceM <- MuraroPancreasData()
sceM <- sceM[,!is.na(sceM$label) & sceM$label!="unclear"] 

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

library(scater)
sceM <- logNormCounts(sceM)

We examine the distribution of labels in this reference.

table(sceM$label)
## 
##      acinar       alpha        beta       delta        duct endothelial 
##         219         812         448         193         245          21 
##     epsilon mesenchymal          pp 
##           3          80         101

We load the Grun et al. (2016) dataset as our test, applying some basic quality control to remove low-quality cells in some of the batches (see here for details).

sceG <- GrunPancreasData()

sceG <- addPerCellQC(sceG)
qc <- quickPerCellQC(colData(sceG), 
    percent_subsets="altexps_ERCC_percent",
    batch=sceG$donor,
    subset=sceG$donor %in% c("D17", "D7", "D2"))
sceG <- sceG[,!qc$discard]

Technically speaking, the test dataset does not need log-expression values but we compute them anyway for convenience.

sceG <- logNormCounts(sceG)

8.2 Applying the annotation

We apply SingleR() with Wilcoxon rank sum test-based marker detection to annotate the Grun dataset with the Muraro labels.

library(SingleR)
pred.grun <- SingleR(test=sceG, ref=sceM, labels=sceM$label, de.method="wilcox")

We examine the distribution of predicted labels:

table(pred.grun$labels)
## 
##      acinar       alpha        beta       delta        duct endothelial 
##         277         203         181          50         306           5 
##     epsilon mesenchymal          pp 
##           1          22          19

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

table(Label=pred.grun$labels,
    Lost=is.na(pred.grun$pruned.labels))
##              Lost
## Label         FALSE TRUE
##   acinar        251   26
##   alpha         198    5
##   beta          180    1
##   delta          49    1
##   duct          301    5
##   endothelial     4    1
##   epsilon         1    0
##   mesenchymal    22    0
##   pp             17    2

8.3 Diagnostics

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

plotScoreHeatmap(pred.grun)
Heatmap of the (normalized) assignment scores for each cell (column) in the Grun test dataset with respect to each label (row) in the Muraro reference dataset. The final assignment for each cell is shown in the annotation bar at the top.

Figure 8.1: Heatmap of the (normalized) assignment scores for each cell (column) in the Grun test dataset with respect to each label (row) in the Muraro 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 8.2.

plotDeltaDistribution(pred.grun)
Distributions of the deltas for each cell in the Grun dataset assigned to each label in the Muraro dataset. Each cell is represented by a point; low-quality assignments that were pruned out are colored in orange.

Figure 8.2: Distributions of the deltas for each cell in the Grun dataset assigned to each label in the Muraro 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 each label in Figure 8.3.

library(scater)
collected <- list()
all.markers <- metadata(pred.grun)$de.genes

sceG$labels <- pred.grun$labels
for (lab in unique(pred.grun$labels)) {
    collected[[lab]] <- plotHeatmap(sceG, silent=TRUE, 
        order_columns_by="labels", main=lab,
        features=unique(unlist(all.markers[[lab]])))[[4]] 
}
do.call(gridExtra::grid.arrange, collected)
Heatmaps of log-expression values in the Grun dataset for all marker genes upregulated in each label in the Muraro reference dataset. Assigned labels for each cell are shown at the top of each plot.

Figure 8.3: Heatmaps of log-expression values in the Grun dataset for all marker genes upregulated in each label in the Muraro reference dataset. Assigned labels for each cell are shown at the top of each plot.

8.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 (increasing the resolution by dropping k=5).

library(scran)
decG <- modelGeneVarWithSpikes(sceG, "ERCC")

set.seed(1000100)
sceG <- denoisePCA(sceG, decG)

library(bluster)
sceG$cluster <- clusterRows(reducedDim(sceG), NNGraphParam(k=5))

We see that the clusters map reasonably well to the labels in Figure 8.4.

tab <- table(cluster=sceG$cluster, label=pred.grun$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 Grun dataset.

Figure 8.4: Heatmap of the log-transformed number of cells in each combination of label (column) and cluster (row) in the Grun dataset.

We proceed to the most important part of the analysis. Yes, that’s right, the \(t\)-SNE plot (Figure 8.5).

set.seed(101010100)
sceG <- runTSNE(sceG, dimred="PCA")
plotTSNE(sceG, colour_by="cluster", text_colour="red",
    text_by=I(pred.grun$labels))
$t$-SNE plot of the Grun dataset, where each point is a cell and is colored by the assigned cluster. Reference labels from the Muraro dataset are also placed on the median coordinate across all cells assigned with that label.

Figure 8.5: \(t\)-SNE plot of the Grun dataset, where each point is a cell and is colored by the assigned cluster. Reference labels from the Muraro dataset are also placed on the median coordinate across all cells assigned with that label.

Session information

R Under development (unstable) (2023-11-11 r85510)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.3 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.13.0              scran_1.31.0               
 [3] SingleR_2.5.0               scater_1.31.1              
 [5] ggplot2_3.4.4               scuttle_1.13.0             
 [7] scRNAseq_2.17.0             SingleCellExperiment_1.25.0
 [9] SummarizedExperiment_1.33.1 Biobase_2.63.0             
[11] GenomicRanges_1.55.1        GenomeInfoDb_1.39.1        
[13] IRanges_2.37.0              S4Vectors_0.41.2           
[15] BiocGenerics_0.49.1         MatrixGenerics_1.15.0      
[17] matrixStats_1.1.0           BiocStyle_2.31.0           
[19] rebook_1.13.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              GenomicFeatures_1.55.1       
  [7] farver_2.1.1                  rmarkdown_2.25               
  [9] BiocIO_1.13.0                 zlibbioc_1.49.0              
 [11] vctrs_0.6.5                   memoise_2.0.1                
 [13] Rsamtools_2.19.2              DelayedMatrixStats_1.25.1    
 [15] RCurl_1.98-1.13               htmltools_0.5.7              
 [17] S4Arrays_1.3.1                progress_1.2.2               
 [19] AnnotationHub_3.11.0          curl_5.1.0                   
 [21] BiocNeighbors_1.21.1          SparseArray_1.3.1            
 [23] sass_0.4.7                    bslib_0.6.1                  
 [25] cachem_1.0.8                  GenomicAlignments_1.39.0     
 [27] igraph_1.5.1                  mime_0.12                    
 [29] lifecycle_1.0.4               pkgconfig_2.0.3              
 [31] rsvd_1.0.5                    Matrix_1.6-4                 
 [33] R6_2.5.1                      fastmap_1.1.1                
 [35] GenomeInfoDbData_1.2.11       shiny_1.8.0                  
 [37] digest_0.6.33                 colorspace_2.1-0             
 [39] AnnotationDbi_1.65.2          dqrng_0.3.2                  
 [41] irlba_2.3.5.1                 ExperimentHub_2.11.0         
 [43] RSQLite_2.3.3                 beachmat_2.19.0              
 [45] labeling_0.4.3                filelock_1.0.2               
 [47] fansi_1.0.5                   httr_1.4.7                   
 [49] abind_1.4-5                   compiler_4.4.0               
 [51] bit64_4.0.5                   withr_2.5.2                  
 [53] BiocParallel_1.37.0           viridis_0.6.4                
 [55] DBI_1.1.3                     highr_0.10                   
 [57] biomaRt_2.59.0                rappdirs_0.3.3               
 [59] DelayedArray_0.29.0           rjson_0.2.21                 
 [61] tools_4.4.0                   vipor_0.4.5                  
 [63] beeswarm_0.4.0                interactiveDisplayBase_1.41.0
 [65] httpuv_1.6.12                 glue_1.6.2                   
 [67] restfulr_0.0.15               promises_1.2.1               
 [69] grid_4.4.0                    Rtsne_0.16                   
 [71] cluster_2.1.6                 generics_0.1.3               
 [73] gtable_0.3.4                  ensembldb_2.27.1             
 [75] hms_1.1.3                     metapod_1.11.0               
 [77] BiocSingular_1.19.0           ScaledMatrix_1.11.0          
 [79] xml2_1.3.6                    utf8_1.2.4                   
 [81] XVector_0.43.0                ggrepel_0.9.4                
 [83] BiocVersion_3.19.1            pillar_1.9.0                 
 [85] stringr_1.5.1                 limma_3.59.1                 
 [87] later_1.3.1                   dplyr_1.1.4                  
 [89] BiocFileCache_2.11.1          lattice_0.22-5               
 [91] rtracklayer_1.63.0            bit_4.0.5                    
 [93] tidyselect_1.2.0              locfit_1.5-9.8               
 [95] Biostrings_2.71.1             knitr_1.45                   
 [97] gridExtra_2.3                 bookdown_0.37                
 [99] ProtGenerics_1.35.0           edgeR_4.1.2                  
[101] xfun_0.41                     statmod_1.5.0                
[103] pheatmap_1.0.12               stringi_1.8.2                
[105] lazyeval_0.2.2                yaml_2.3.7                   
[107] evaluate_0.23                 codetools_0.2-19             
[109] tibble_3.2.1                  BiocManager_1.30.22          
[111] graph_1.81.0                  cli_3.6.1                    
[113] xtable_1.8-4                  munsell_0.5.0                
[115] jquerylib_0.1.4               Rcpp_1.0.11                  
[117] dir.expiry_1.11.0             dbplyr_2.4.0                 
[119] png_0.1-8                     XML_3.99-0.16                
[121] parallel_4.4.0                ellipsis_0.3.2               
[123] blob_1.2.4                    prettyunits_1.2.0            
[125] AnnotationFilter_1.27.0       sparseMatrixStats_1.15.0     
[127] bitops_1.0-7                  viridisLite_0.4.2            
[129] scales_1.3.0                  purrr_1.0.2                  
[131] crayon_1.5.2                  rlang_1.1.2                  
[133] cowplot_1.1.1                 KEGGREST_1.43.0              

Bibliography

Grun, D., M. J. Muraro, J. C. Boisset, K. Wiebrands, A. Lyubimova, G. Dharmadhikari, M. van den Born, et al. 2016. “De Novo Prediction of Stem Cell Identity using Single-Cell Transcriptome Data.” Cell Stem Cell 19 (2): 266–77.

Muraro, M. J., G. Dharmadhikari, D. Grun, N. Groen, T. Dielen, E. Jansen, L. van Gurp, et al. 2016. “A Single-Cell Transcriptome Atlas of the Human Pancreas.” Cell Syst 3 (4): 385–94.