# Grun human pancreas (CEL-seq2) ## Introduction This workflow performs an analysis of the @grun2016denovo CEL-seq2 dataset consisting of human pancreas cells from various donors. ## Data loading ``` r library(scRNAseq) sce.grun <- GrunPancreasData() ``` We convert to Ensembl identifiers, and we remove duplicated genes or genes without Ensembl IDs. ``` r library(org.Hs.eg.db) gene.ids <- mapIds(org.Hs.eg.db, keys=rowData(sce.grun)$symbol, keytype="SYMBOL", column="ENSEMBL") keep <- !is.na(gene.ids) & !duplicated(gene.ids) sce.grun <- sce.grun[keep,] rownames(sce.grun) <- gene.ids[keep] ``` ## Quality control ``` r unfiltered <- sce.grun ``` This dataset lacks mitochondrial genes so we will do without them for quality control. We compute the median and MAD while blocking on the donor; for donors where the assumption of a majority of high-quality cells seems to be violated (Figure \@ref(fig:unref-grun-qc-dist)), we compute an appropriate threshold using the other donors as specified in the `subset=` argument. ``` r library(scater) stats <- perCellQCMetrics(sce.grun) qc <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent", batch=sce.grun$donor, subset=sce.grun$donor %in% c("D17", "D7", "D2")) sce.grun <- sce.grun[,!qc$discard] ``` ``` r colData(unfiltered) <- cbind(colData(unfiltered), stats) unfiltered$discard <- qc$discard gridExtra::grid.arrange( plotColData(unfiltered, x="donor", y="sum", colour_by="discard") + scale_y_log10() + ggtitle("Total count"), plotColData(unfiltered, x="donor", y="detected", colour_by="discard") + scale_y_log10() + ggtitle("Detected features"), plotColData(unfiltered, x="donor", y="altexps_ERCC_percent", colour_by="discard") + ggtitle("ERCC percent"), ncol=2 ) ```
Distribution of each QC metric across cells from each donor of the Grun pancreas dataset. Each point represents a cell and is colored according to whether that cell was discarded.

(\#fig:unref-grun-qc-dist)Distribution of each QC metric across cells from each donor of the Grun pancreas dataset. Each point represents a cell and is colored according to whether that cell was discarded.

``` r colSums(as.matrix(qc), na.rm=TRUE) ``` ``` ## low_lib_size low_n_features high_altexps_ERCC_percent ## 447 511 605 ## discard ## 664 ``` ## Normalization ``` r library(scran) set.seed(1000) # for irlba. clusters <- quickCluster(sce.grun) sce.grun <- computeSumFactors(sce.grun, clusters=clusters) sce.grun <- logNormCounts(sce.grun) ``` ``` r summary(sizeFactors(sce.grun)) ``` ``` ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.0949 0.5015 0.7962 1.0000 1.2263 9.4546 ``` ``` r plot(librarySizeFactors(sce.grun), sizeFactors(sce.grun), pch=16, xlab="Library size factors", ylab="Deconvolution factors", log="xy") ```
Relationship between the library size factors and the deconvolution size factors in the Grun pancreas dataset.

(\#fig:unref-grun-norm)Relationship between the library size factors and the deconvolution size factors in the Grun pancreas dataset.

## Variance modelling We block on a combined plate and donor factor. ``` r block <- paste0(sce.grun$sample, "_", sce.grun$donor) dec.grun <- modelGeneVarWithSpikes(sce.grun, spikes="ERCC", block=block) top.grun <- getTopHVGs(dec.grun, prop=0.1) ``` We examine the number of cells in each level of the blocking factor. ``` r table(block) ``` ``` ## block ## CD13+ sorted cells_D17 CD24+ CD44+ live sorted cells_D17 ## 86 87 ## CD63+ sorted cells_D10 TGFBR3+ sorted cells_D17 ## 40 90 ## exocrine fraction, live sorted cells_D2 exocrine fraction, live sorted cells_D3 ## 82 7 ## live sorted cells, library 1_D10 live sorted cells, library 1_D17 ## 33 88 ## live sorted cells, library 1_D3 live sorted cells, library 1_D7 ## 25 85 ## live sorted cells, library 2_D10 live sorted cells, library 2_D17 ## 35 83 ## live sorted cells, library 2_D3 live sorted cells, library 2_D7 ## 27 84 ## live sorted cells, library 3_D3 live sorted cells, library 3_D7 ## 17 83 ## live sorted cells, library 4_D3 live sorted cells, library 4_D7 ## 29 83 ``` ``` r par(mfrow=c(6,3)) blocked.stats <- dec.grun$per.block for (i in colnames(blocked.stats)) { current <- blocked.stats[[i]] plot(current$mean, current$total, main=i, pch=16, cex=0.5, xlab="Mean of log-expression", ylab="Variance of log-expression") curfit <- metadata(current) points(curfit$mean, curfit$var, col="red", pch=16) curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2) } ```
Per-gene variance as a function of the mean for the log-expression values in the Grun pancreas dataset. Each point represents a gene (black) with the mean-variance trend (blue) fitted to the spike-in transcripts (red) separately for each donor.

(\#fig:unref-416b-variance)Per-gene variance as a function of the mean for the log-expression values in the Grun pancreas dataset. Each point represents a gene (black) with the mean-variance trend (blue) fitted to the spike-in transcripts (red) separately for each donor.

## Data integration ``` r library(batchelor) set.seed(1001010) merged.grun <- fastMNN(sce.grun, subset.row=top.grun, batch=sce.grun$donor) ``` ``` r metadata(merged.grun)$merge.info$lost.var ``` ``` ## D10 D17 D2 D3 D7 ## [1,] 0.029320 0.030835 0.000000 0.00000 0.00000 ## [2,] 0.007790 0.011984 0.038270 0.00000 0.00000 ## [3,] 0.004175 0.005056 0.008115 0.05288 0.00000 ## [4,] 0.014759 0.017827 0.016893 0.01534 0.05513 ``` ## Dimensionality reduction ``` r set.seed(100111) merged.grun <- runTSNE(merged.grun, dimred="corrected") ``` ## Clustering ``` r snn.gr <- buildSNNGraph(merged.grun, use.dimred="corrected") colLabels(merged.grun) <- factor(igraph::cluster_walktrap(snn.gr)$membership) ``` ``` r table(Cluster=colLabels(merged.grun), Donor=merged.grun$batch) ``` ``` ## Donor ## Cluster D10 D17 D2 D3 D7 ## 1 12 128 0 0 62 ## 2 32 70 31 81 29 ## 3 11 73 31 2 70 ## 4 14 33 3 2 68 ## 5 2 8 3 3 6 ## 6 4 4 2 4 2 ## 7 3 41 0 0 9 ## 8 16 34 12 11 45 ## 9 5 13 0 0 10 ## 10 4 13 0 0 1 ## 11 5 17 0 2 33 ``` ``` r gridExtra::grid.arrange( plotTSNE(merged.grun, colour_by="label"), plotTSNE(merged.grun, colour_by="batch"), ncol=2 ) ```
Obligatory $t$-SNE plots of the Grun pancreas dataset. Each point represents a cell that is colored by cluster (left) or batch (right).

(\#fig:unref-grun-tsne)Obligatory $t$-SNE plots of the Grun pancreas dataset. Each point represents a cell that is colored by cluster (left) or batch (right).

## Session Info {-}
``` R version 4.5.0 RC (2025-04-04 r88126) Platform: x86_64-pc-linux-gnu Running under: Ubuntu 24.04.2 LTS Matrix products: default BLAS: /home/biocbuild/bbs-3.21-bioc/R/lib/libRblas.so LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.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] batchelor_1.24.0 scran_1.36.0 [3] scater_1.36.0 ggplot2_3.5.2 [5] scuttle_1.18.0 org.Hs.eg.db_3.21.0 [7] AnnotationDbi_1.70.0 scRNAseq_2.21.1 [9] SingleCellExperiment_1.30.0 SummarizedExperiment_1.38.0 [11] Biobase_2.68.0 GenomicRanges_1.60.0 [13] GenomeInfoDb_1.44.0 IRanges_2.42.0 [15] S4Vectors_0.46.0 BiocGenerics_0.54.0 [17] generics_0.1.3 MatrixGenerics_1.20.0 [19] matrixStats_1.5.0 BiocStyle_2.36.0 [21] rebook_1.18.0 loaded via a namespace (and not attached): [1] jsonlite_2.0.0 CodeDepends_0.6.6 [3] magrittr_2.0.3 ggbeeswarm_0.7.2 [5] GenomicFeatures_1.60.0 gypsum_1.4.0 [7] farver_2.1.2 rmarkdown_2.29 [9] BiocIO_1.18.0 vctrs_0.6.5 [11] DelayedMatrixStats_1.30.0 memoise_2.0.1 [13] Rsamtools_2.24.0 RCurl_1.98-1.17 [15] htmltools_0.5.8.1 S4Arrays_1.8.0 [17] AnnotationHub_3.16.0 curl_6.2.2 [19] BiocNeighbors_2.2.0 Rhdf5lib_1.30.0 [21] SparseArray_1.8.0 rhdf5_2.52.0 [23] sass_0.4.10 alabaster.base_1.8.0 [25] bslib_0.9.0 alabaster.sce_1.8.0 [27] httr2_1.1.2 cachem_1.1.0 [29] ResidualMatrix_1.18.0 GenomicAlignments_1.44.0 [31] igraph_2.1.4 lifecycle_1.0.4 [33] pkgconfig_2.0.3 rsvd_1.0.5 [35] Matrix_1.7-3 R6_2.6.1 [37] fastmap_1.2.0 GenomeInfoDbData_1.2.14 [39] digest_0.6.37 colorspace_2.1-1 [41] dqrng_0.4.1 irlba_2.3.5.1 [43] ExperimentHub_2.16.0 RSQLite_2.3.9 [45] beachmat_2.24.0 labeling_0.4.3 [47] filelock_1.0.3 httr_1.4.7 [49] abind_1.4-8 compiler_4.5.0 [51] bit64_4.6.0-1 withr_3.0.2 [53] BiocParallel_1.42.0 viridis_0.6.5 [55] DBI_1.2.3 HDF5Array_1.36.0 [57] alabaster.ranges_1.8.0 alabaster.schemas_1.8.0 [59] rappdirs_0.3.3 DelayedArray_0.34.0 [61] bluster_1.18.0 rjson_0.2.23 [63] tools_4.5.0 vipor_0.4.7 [65] beeswarm_0.4.0 glue_1.8.0 [67] h5mread_1.0.0 restfulr_0.0.15 [69] rhdf5filters_1.20.0 grid_4.5.0 [71] Rtsne_0.17 cluster_2.1.8.1 [73] gtable_0.3.6 ensembldb_2.32.0 [75] metapod_1.16.0 BiocSingular_1.24.0 [77] ScaledMatrix_1.16.0 XVector_0.48.0 [79] ggrepel_0.9.6 BiocVersion_3.21.1 [81] pillar_1.10.2 limma_3.64.0 [83] dplyr_1.1.4 BiocFileCache_2.16.0 [85] lattice_0.22-7 rtracklayer_1.68.0 [87] bit_4.6.0 tidyselect_1.2.1 [89] locfit_1.5-9.12 Biostrings_2.76.0 [91] knitr_1.50 gridExtra_2.3 [93] bookdown_0.43 ProtGenerics_1.40.0 [95] edgeR_4.6.0 xfun_0.52 [97] statmod_1.5.0 UCSC.utils_1.4.0 [99] lazyeval_0.2.2 yaml_2.3.10 [101] evaluate_1.0.3 codetools_0.2-20 [103] tibble_3.2.1 alabaster.matrix_1.8.0 [105] BiocManager_1.30.25 graph_1.86.0 [107] cli_3.6.4 munsell_0.5.1 [109] jquerylib_0.1.4 Rcpp_1.0.14 [111] dir.expiry_1.16.0 dbplyr_2.5.0 [113] png_0.1-8 XML_3.99-0.18 [115] parallel_4.5.0 blob_1.2.4 [117] AnnotationFilter_1.32.0 sparseMatrixStats_1.20.0 [119] bitops_1.0-9 viridisLite_0.4.2 [121] alabaster.se_1.8.0 scales_1.3.0 [123] crayon_1.5.3 rlang_1.1.6 [125] cowplot_1.1.3 KEGGREST_1.48.0 ```