Chapter 5 Grun human pancreas (CEL-seq2)

5.1 Introduction

This workflow performs an analysis of the Grun et al. (2016) CEL-seq2 dataset consisting of human pancreas cells from various donors.

5.2 Data loading

library(scRNAseq)
sce.grun <- GrunPancreasData()

We convert to Ensembl identifiers, and we remove duplicated genes or genes without Ensembl IDs.

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]

5.3 Quality control

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 5.1), we compute an appropriate threshold using the other donors as specified in the subset= argument.

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]
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.

Figure 5.1: 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.

colSums(as.matrix(qc), na.rm=TRUE)
##              low_lib_size            low_n_features high_altexps_ERCC_percent 
##                       447                       511                       605 
##                   discard 
##                       664

5.4 Normalization

library(scran)
set.seed(1000) # for irlba. 
clusters <- quickCluster(sce.grun)
sce.grun <- computeSumFactors(sce.grun, clusters=clusters)
sce.grun <- logNormCounts(sce.grun)
summary(sizeFactors(sce.grun))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0949  0.5015  0.7962  1.0000  1.2263  9.4546
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.

Figure 5.2: Relationship between the library size factors and the deconvolution size factors in the Grun pancreas dataset.

5.5 Variance modelling

We block on a combined plate and donor factor.

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.

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

Figure 1.4: 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.

5.6 Data integration

library(batchelor)
set.seed(1001010)
merged.grun <- fastMNN(sce.grun, subset.row=top.grun, batch=sce.grun$donor)
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

5.7 Dimensionality reduction

set.seed(100111)
merged.grun <- runTSNE(merged.grun, dimred="corrected")

5.8 Clustering

snn.gr <- buildSNNGraph(merged.grun, use.dimred="corrected")
colLabels(merged.grun) <- factor(igraph::cluster_walktrap(snn.gr)$membership)
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
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).

Figure 5.3: 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          

References

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.