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
We convert to Ensembl identifiers, and we remove duplicated genes or genes without Ensembl IDs.
5.3 Quality control
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
)
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.
## low_lib_size low_n_features high_altexps_ERCC_percent
## 452 510 606
## discard
## 6655.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)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.094 0.507 0.794 1.000 1.235 10.953plot(librarySizeFactors(sce.grun), sizeFactors(sce.grun), pch=16,
xlab="Library size factors", ylab="Deconvolution factors", log="xy")
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.
## 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
## 16 83
## live sorted cells, library 4_D3 live sorted cells, library 4_D7
## 29 83par(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)
}
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)## D10 D17 D2 D3 D7
## [1,] 0.029789 0.031754 0.000000 0.00000 0.00000
## [2,] 0.008008 0.012371 0.039101 0.00000 0.00000
## [3,] 0.004108 0.005397 0.008157 0.05204 0.00000
## [4,] 0.013393 0.016061 0.016364 0.01510 0.055225.8 Clustering
snn.gr <- buildSNNGraph(merged.grun, use.dimred="corrected")
colLabels(merged.grun) <- factor(igraph::cluster_walktrap(snn.gr)$membership)## Donor
## Cluster D10 D17 D2 D3 D7
## 1 32 71 33 80 29
## 2 11 119 0 0 55
## 3 2 7 3 3 6
## 4 3 43 0 0 11
## 5 5 14 0 0 10
## 6 4 4 2 4 2
## 7 11 69 29 3 69
## 8 16 37 12 10 46
## 9 14 31 3 2 66
## 10 1 9 0 0 7
## 11 4 13 0 0 1
## 12 5 17 0 2 33gridExtra::grid.arrange(
plotTSNE(merged.grun, colour_by="label"),
plotTSNE(merged.grun, colour_by="batch"),
ncol=2
)
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.3.0 RC (2023-04-13 r84269)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.2 LTS
Matrix products: default
BLAS: /home/biocbuild/bbs-3.17-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] batchelor_1.16.0 scran_1.28.0
[3] scater_1.28.0 ggplot2_3.4.2
[5] scuttle_1.10.0 org.Hs.eg.db_3.17.0
[7] AnnotationDbi_1.62.0 scRNAseq_2.13.0
[9] SingleCellExperiment_1.22.0 SummarizedExperiment_1.30.0
[11] Biobase_2.60.0 GenomicRanges_1.52.0
[13] GenomeInfoDb_1.36.0 IRanges_2.34.0
[15] S4Vectors_0.38.0 BiocGenerics_0.46.0
[17] MatrixGenerics_1.12.0 matrixStats_0.63.0
[19] BiocStyle_2.28.0 rebook_1.10.0
loaded via a namespace (and not attached):
[1] jsonlite_1.8.4 CodeDepends_0.6.5
[3] magrittr_2.0.3 ggbeeswarm_0.7.1
[5] GenomicFeatures_1.52.0 farver_2.1.1
[7] rmarkdown_2.21 BiocIO_1.10.0
[9] zlibbioc_1.46.0 vctrs_0.6.2
[11] memoise_2.0.1 Rsamtools_2.16.0
[13] DelayedMatrixStats_1.22.0 RCurl_1.98-1.12
[15] htmltools_0.5.5 progress_1.2.2
[17] AnnotationHub_3.8.0 curl_5.0.0
[19] BiocNeighbors_1.18.0 sass_0.4.5
[21] bslib_0.4.2 cachem_1.0.7
[23] ResidualMatrix_1.10.0 GenomicAlignments_1.36.0
[25] igraph_1.4.2 mime_0.12
[27] lifecycle_1.0.3 pkgconfig_2.0.3
[29] rsvd_1.0.5 Matrix_1.5-4
[31] R6_2.5.1 fastmap_1.1.1
[33] GenomeInfoDbData_1.2.10 shiny_1.7.4
[35] digest_0.6.31 colorspace_2.1-0
[37] dqrng_0.3.0 irlba_2.3.5.1
[39] ExperimentHub_2.8.0 RSQLite_2.3.1
[41] beachmat_2.16.0 labeling_0.4.2
[43] filelock_1.0.2 fansi_1.0.4
[45] httr_1.4.5 compiler_4.3.0
[47] bit64_4.0.5 withr_2.5.0
[49] BiocParallel_1.34.0 viridis_0.6.2
[51] DBI_1.1.3 highr_0.10
[53] biomaRt_2.56.0 rappdirs_0.3.3
[55] DelayedArray_0.26.0 bluster_1.10.0
[57] rjson_0.2.21 tools_4.3.0
[59] vipor_0.4.5 beeswarm_0.4.0
[61] interactiveDisplayBase_1.38.0 httpuv_1.6.9
[63] glue_1.6.2 restfulr_0.0.15
[65] promises_1.2.0.1 grid_4.3.0
[67] Rtsne_0.16 cluster_2.1.4
[69] generics_0.1.3 gtable_0.3.3
[71] ensembldb_2.24.0 hms_1.1.3
[73] metapod_1.8.0 BiocSingular_1.16.0
[75] ScaledMatrix_1.8.0 xml2_1.3.3
[77] utf8_1.2.3 XVector_0.40.0
[79] ggrepel_0.9.3 BiocVersion_3.17.1
[81] pillar_1.9.0 stringr_1.5.0
[83] limma_3.56.0 later_1.3.0
[85] dplyr_1.1.2 BiocFileCache_2.8.0
[87] lattice_0.21-8 rtracklayer_1.60.0
[89] bit_4.0.5 tidyselect_1.2.0
[91] locfit_1.5-9.7 Biostrings_2.68.0
[93] knitr_1.42 gridExtra_2.3
[95] bookdown_0.33 ProtGenerics_1.32.0
[97] edgeR_3.42.0 xfun_0.39
[99] statmod_1.5.0 stringi_1.7.12
[101] lazyeval_0.2.2 yaml_2.3.7
[103] evaluate_0.20 codetools_0.2-19
[105] tibble_3.2.1 BiocManager_1.30.20
[107] graph_1.78.0 cli_3.6.1
[109] xtable_1.8-4 munsell_0.5.0
[111] jquerylib_0.1.4 Rcpp_1.0.10
[113] dir.expiry_1.8.0 dbplyr_2.3.2
[115] png_0.1-8 XML_3.99-0.14
[117] parallel_4.3.0 ellipsis_0.3.2
[119] blob_1.2.4 prettyunits_1.1.1
[121] AnnotationFilter_1.24.0 sparseMatrixStats_1.12.0
[123] bitops_1.0-7 viridisLite_0.4.1
[125] scales_1.2.1 purrr_1.0.1
[127] crayon_1.5.2 rlang_1.1.0
[129] cowplot_1.1.1 KEGGREST_1.40.0 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.