# 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.1 (2025-06-13) Platform: x86_64-pc-linux-gnu Running under: Ubuntu 24.04.3 LTS Matrix products: default BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; 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.25.0 scran_1.37.0 [3] scater_1.37.0 ggplot2_3.5.2 [5] scuttle_1.19.0 org.Hs.eg.db_3.21.0 [7] AnnotationDbi_1.71.1 scRNAseq_2.23.0 [9] SingleCellExperiment_1.31.1 SummarizedExperiment_1.39.1 [11] Biobase_2.69.0 GenomicRanges_1.61.1 [13] Seqinfo_0.99.2 IRanges_2.43.0 [15] S4Vectors_0.47.0 BiocGenerics_0.55.1 [17] generics_0.1.4 MatrixGenerics_1.21.0 [19] matrixStats_1.5.0 BiocStyle_2.37.1 [21] rebook_1.19.0 loaded via a namespace (and not attached): [1] RColorBrewer_1.1-3 jsonlite_2.0.0 [3] CodeDepends_0.6.6 magrittr_2.0.3 [5] ggbeeswarm_0.7.2 GenomicFeatures_1.61.6 [7] gypsum_1.5.0 farver_2.1.2 [9] rmarkdown_2.29 BiocIO_1.19.0 [11] vctrs_0.6.5 DelayedMatrixStats_1.31.0 [13] memoise_2.0.1 Rsamtools_2.25.2 [15] RCurl_1.98-1.17 htmltools_0.5.8.1 [17] S4Arrays_1.9.1 AnnotationHub_3.99.6 [19] curl_6.4.0 BiocNeighbors_2.3.1 [21] Rhdf5lib_1.31.0 SparseArray_1.9.1 [23] rhdf5_2.53.4 sass_0.4.10 [25] alabaster.base_1.9.5 bslib_0.9.0 [27] alabaster.sce_1.9.0 httr2_1.2.1 [29] cachem_1.1.0 ResidualMatrix_1.19.0 [31] GenomicAlignments_1.45.2 igraph_2.1.4 [33] lifecycle_1.0.4 pkgconfig_2.0.3 [35] rsvd_1.0.5 Matrix_1.7-3 [37] R6_2.6.1 fastmap_1.2.0 [39] digest_0.6.37 dqrng_0.4.1 [41] irlba_2.3.5.1 ExperimentHub_2.99.5 [43] RSQLite_2.4.2 beachmat_2.25.4 [45] labeling_0.4.3 filelock_1.0.3 [47] httr_1.4.7 abind_1.4-8 [49] compiler_4.5.1 bit64_4.6.0-1 [51] withr_3.0.2 BiocParallel_1.43.4 [53] viridis_0.6.5 DBI_1.2.3 [55] HDF5Array_1.37.0 alabaster.ranges_1.9.1 [57] alabaster.schemas_1.9.0 rappdirs_0.3.3 [59] DelayedArray_0.35.2 bluster_1.19.0 [61] rjson_0.2.23 tools_4.5.1 [63] vipor_0.4.7 beeswarm_0.4.0 [65] glue_1.8.0 h5mread_1.1.1 [67] restfulr_0.0.16 rhdf5filters_1.21.0 [69] grid_4.5.1 Rtsne_0.17 [71] cluster_2.1.8.1 gtable_0.3.6 [73] ensembldb_2.33.1 metapod_1.17.0 [75] BiocSingular_1.25.0 ScaledMatrix_1.17.0 [77] XVector_0.49.0 ggrepel_0.9.6 [79] BiocVersion_3.22.0 pillar_1.11.0 [81] limma_3.65.3 dplyr_1.1.4 [83] BiocFileCache_2.99.5 lattice_0.22-7 [85] rtracklayer_1.69.1 bit_4.6.0 [87] tidyselect_1.2.1 locfit_1.5-9.12 [89] Biostrings_2.77.2 knitr_1.50 [91] gridExtra_2.3 bookdown_0.43 [93] ProtGenerics_1.41.0 edgeR_4.7.3 [95] xfun_0.52 statmod_1.5.0 [97] UCSC.utils_1.5.0 lazyeval_0.2.2 [99] yaml_2.3.10 evaluate_1.0.4 [101] codetools_0.2-20 tibble_3.3.0 [103] alabaster.matrix_1.9.0 BiocManager_1.30.26 [105] graph_1.87.0 cli_3.6.5 [107] jquerylib_0.1.4 dichromat_2.0-0.1 [109] Rcpp_1.1.0 GenomeInfoDb_1.45.9 [111] dir.expiry_1.17.0 dbplyr_2.5.0 [113] png_0.1-8 XML_3.99-0.18 [115] parallel_4.5.1 blob_1.2.4 [117] AnnotationFilter_1.33.0 sparseMatrixStats_1.21.0 [119] bitops_1.0-9 viridisLite_0.4.2 [121] alabaster.se_1.9.0 scales_1.4.0 [123] crayon_1.5.3 rlang_1.1.6 [125] cowplot_1.2.0 KEGGREST_1.49.1 ```