Chapter 13 Messmer human ESC (Smart-seq2)
13.1 Introduction
This performs an analysis of the human embryonic stem cell (hESC) dataset generated with Smart-seq2 (Messmer et al. 2019), which contains several plates of naive and primed hESCs. The chapter’s code is based on the steps in the paper’s GitHub repository, with some additional steps for cell cycle effect removal contributed by Philippe Boileau.
13.2 Data loading
Converting the batch to a factor, to make life easier later on.
13.3 Quality control
Let’s have a look at the QC statistics.
## low_lib_size low_n_features high_subsets_Mito_percent
## 107 99 22
## high_altexps_ERCC_percent discard
## 117 156gridExtra::grid.arrange(
plotColData(original, x="experiment batch", y="sum",
colour_by=I(filtered$discard), other_field="phenotype") +
facet_wrap(~phenotype) + scale_y_log10(),
plotColData(original, x="experiment batch", y="detected",
colour_by=I(filtered$discard), other_field="phenotype") +
facet_wrap(~phenotype) + scale_y_log10(),
plotColData(original, x="experiment batch", y="subsets_Mito_percent",
colour_by=I(filtered$discard), other_field="phenotype") +
facet_wrap(~phenotype),
plotColData(original, x="experiment batch", y="altexps_ERCC_percent",
colour_by=I(filtered$discard), other_field="phenotype") +
facet_wrap(~phenotype),
ncol=1
)
Figure 13.1: Distribution of QC metrics across batches (x-axis) and phenotypes (facets) for cells in the Messmer hESC dataset. Each point is a cell and is colored by whether it was discarded.
13.4 Normalization
library(scran)
set.seed(10000)
clusters <- quickCluster(sce.mess)
sce.mess <- computeSumFactors(sce.mess, cluster=clusters)
sce.mess <- logNormCounts(sce.mess)par(mfrow=c(1,2))
plot(sce.mess$sum, sizeFactors(sce.mess), log = "xy", pch=16,
xlab = "Library size (millions)", ylab = "Size factor",
col = ifelse(sce.mess$phenotype == "naive", "black", "grey"))
spike.sf <- librarySizeFactors(altExp(sce.mess, "ERCC"))
plot(sizeFactors(sce.mess), spike.sf, log = "xy", pch=16,
ylab = "Spike-in size factor", xlab = "Deconvolution size factor",
col = ifelse(sce.mess$phenotype == "naive", "black", "grey"))
Figure 13.2: Deconvolution size factors plotted against the library size (left) and spike-in size factors plotted against the deconvolution size factors (right). Each point is a cell and is colored by its phenotype.
13.5 Cell cycle phase assignment
Here, we use multiple cores to speed up the processing.
set.seed(10001)
hs_pairs <- readRDS(system.file("exdata", "human_cycle_markers.rds", package="scran"))
assigned <- cyclone(sce.mess, pairs=hs_pairs,
gene.names=rownames(sce.mess),
BPPARAM=BiocParallel::MulticoreParam(10))
sce.mess$phase <- assigned$phases##
## G1 G2M S
## 460 406 322
Figure 13.3: G1 cyclone() phase scores against the G2/M phase scores for each cell in the Messmer hESC dataset.
13.6 Feature selection
dec <- modelGeneVarWithSpikes(sce.mess, "ERCC", block = sce.mess$`experiment batch`)
top.hvgs <- getTopHVGs(dec, prop = 0.1)par(mfrow=c(1,3))
for (i in seq_along(dec$per.block)) {
current <- dec$per.block[[i]]
plot(current$mean, current$total, xlab="Mean log-expression",
ylab="Variance", pch=16, cex=0.5, main=paste("Batch", i))
fit <- metadata(current)
points(fit$mean, fit$var, col="red", pch=16)
curve(fit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
}
Figure 13.4: Per-gene variance of the log-normalized expression values in the Messmer hESC dataset, plotted against the mean for each batch. Each point represents a gene with spike-ins shown in red and the fitted trend shown in blue.
13.7 Batch correction
We eliminate the obvious batch effect between batches with linear regression, which is possible due to the replicated nature of the experimental design.
We set keep=1:2 to retain the effect of the first two coefficients in design corresponding to our phenotype of interest.
13.8 Dimensionality Reduction
We could have set d= and subset.row= in correctExperiments() to automatically perform a PCA on the the residual matrix with the subset of HVGs,
but we’ll just explicitly call runPCA() here to keep things simple.
set.seed(1101001)
sce.mess <- runPCA(sce.mess, subset_row = top.hvgs, exprs_values = "corrected")
sce.mess <- runTSNE(sce.mess, dimred = "PCA", perplexity = 40)From a naive PCA, the cell cycle appears to be a major source of biological variation within each phenotype.
gridExtra::grid.arrange(
plotTSNE(sce.mess, colour_by = "phenotype") + ggtitle("By phenotype"),
plotTSNE(sce.mess, colour_by = "experiment batch") + ggtitle("By batch "),
plotTSNE(sce.mess, colour_by = "CDK1", swap_rownames="SYMBOL") + ggtitle("By CDK1"),
plotTSNE(sce.mess, colour_by = "phase") + ggtitle("By phase"),
ncol = 2
)
Figure 13.5: Obligatory \(t\)-SNE plots of the Messmer hESC dataset, where each point is a cell and is colored by various attributes.
We perform contrastive PCA (cPCA) and sparse cPCA (scPCA) on the corrected log-expression data to obtain the same number of PCs. Given that the naive hESCs are actually reprogrammed primed hESCs, we will use the single batch of primed-only hESCs as the “background” dataset to remove the cell cycle effect.
library(scPCA)
is.bg <- sce.mess$`experiment batch`=="3"
target <- sce.mess[,!is.bg]
background <- sce.mess[,is.bg]
mat.target <- t(assay(target, "corrected")[top.hvgs,])
mat.background <- t(assay(background, "corrected")[top.hvgs,])
set.seed(1010101001)
con_out <- scPCA(
target = mat.target,
background = mat.background,
penalties = 0, # no penalties = non-sparse cPCA.
n_eigen = 50,
contrasts = 100
)
reducedDim(target, "cPCA") <- con_out$xset.seed(101010101)
sparse_con_out <- scPCA(
target = mat.target,
background = mat.background,
penalties = 1e-4,
n_eigen = 50,
contrasts = 100,
alg = "rand_var_proj" # for speed.
)
reducedDim(target, "scPCA") <- sparse_con_out$xWe see greater intermingling between phases within both the naive and primed cells after cPCA and scPCA.
set.seed(1101001)
target <- runTSNE(target, dimred = "cPCA", perplexity = 40, name="cPCA+TSNE")
target <- runTSNE(target, dimred = "scPCA", perplexity = 40, name="scPCA+TSNE")gridExtra::grid.arrange(
plotReducedDim(target, "cPCA+TSNE", colour_by = "phase") + ggtitle("After cPCA"),
plotReducedDim(target, "scPCA+TSNE", colour_by = "phase") + ggtitle("After scPCA"),
ncol=2
)
Figure 13.6: More \(t\)-SNE plots of the Messmer hESC dataset after cPCA and scPCA, where each point is a cell and is colored by its assigned cell cycle phase.
We can quantify the change in the separation between phases within each phenotype using the silhouette coefficient.
library(bluster)
naive <- target[,target$phenotype=="naive"]
primed <- target[,target$phenotype=="primed"]
N <- approxSilhouette(reducedDim(naive, "PCA"), naive$phase)
P <- approxSilhouette(reducedDim(primed, "PCA"), primed$phase)
c(naive=mean(N$width), primed=mean(P$width))## naive primed
## 0.02032 0.03025cN <- approxSilhouette(reducedDim(naive, "cPCA"), naive$phase)
cP <- approxSilhouette(reducedDim(primed, "cPCA"), primed$phase)
c(naive=mean(cN$width), primed=mean(cP$width))## naive primed
## 0.007696 0.011941scN <- approxSilhouette(reducedDim(naive, "scPCA"), naive$phase)
scP <- approxSilhouette(reducedDim(primed, "scPCA"), primed$phase)
c(naive=mean(scN$width), primed=mean(scP$width))## naive primed
## 0.006614 0.014601Session Info
R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.1 LTS
Matrix products: default
BLAS: /home/biocbuild/bbs-3.20-bioc/R/lib/libRblas.so
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.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] bluster_1.16.0 scPCA_1.20.0
[3] batchelor_1.22.0 scran_1.34.0
[5] scater_1.34.0 ggplot2_3.5.1
[7] scuttle_1.16.0 AnnotationHub_3.14.0
[9] BiocFileCache_2.14.0 dbplyr_2.5.0
[11] ensembldb_2.30.0 AnnotationFilter_1.30.0
[13] GenomicFeatures_1.58.0 AnnotationDbi_1.68.0
[15] scRNAseq_2.19.1 SingleCellExperiment_1.28.0
[17] SummarizedExperiment_1.36.0 Biobase_2.66.0
[19] GenomicRanges_1.58.0 GenomeInfoDb_1.42.0
[21] IRanges_2.40.0 S4Vectors_0.44.0
[23] BiocGenerics_0.52.0 MatrixGenerics_1.18.0
[25] matrixStats_1.4.1 BiocStyle_2.34.0
[27] rebook_1.16.0
loaded via a namespace (and not attached):
[1] BiocIO_1.16.0 bitops_1.0-9
[3] filelock_1.0.3 tibble_3.2.1
[5] CodeDepends_0.6.6 graph_1.84.0
[7] XML_3.99-0.17 lifecycle_1.0.4
[9] httr2_1.0.5 Rdpack_2.6.1
[11] edgeR_4.4.0 globals_0.16.3
[13] lattice_0.22-6 alabaster.base_1.6.0
[15] magrittr_2.0.3 limma_3.62.0
[17] sass_0.4.9 rmarkdown_2.28
[19] jquerylib_0.1.4 yaml_2.3.10
[21] metapod_1.14.0 cowplot_1.1.3
[23] DBI_1.2.3 ResidualMatrix_1.16.0
[25] abind_1.4-8 zlibbioc_1.52.0
[27] Rtsne_0.17 purrr_1.0.2
[29] RCurl_1.98-1.16 rappdirs_0.3.3
[31] GenomeInfoDbData_1.2.13 ggrepel_0.9.6
[33] irlba_2.3.5.1 listenv_0.9.1
[35] alabaster.sce_1.6.0 RSpectra_0.16-2
[37] parallelly_1.38.0 dqrng_0.4.1
[39] DelayedMatrixStats_1.28.0 codetools_0.2-20
[41] DelayedArray_0.32.0 tidyselect_1.2.1
[43] UCSC.utils_1.2.0 farver_2.1.2
[45] ScaledMatrix_1.14.0 viridis_0.6.5
[47] GenomicAlignments_1.42.0 jsonlite_1.8.9
[49] BiocNeighbors_2.0.0 tools_4.4.1
[51] Rcpp_1.0.13 glue_1.8.0
[53] gridExtra_2.3 SparseArray_1.6.0
[55] xfun_0.48 dplyr_1.1.4
[57] HDF5Array_1.34.0 gypsum_1.2.0
[59] withr_3.0.2 BiocManager_1.30.25
[61] fastmap_1.2.0 sparsepca_0.1.2
[63] rhdf5filters_1.18.0 fansi_1.0.6
[65] digest_0.6.37 rsvd_1.0.5
[67] R6_2.5.1 mime_0.12
[69] colorspace_2.1-1 RSQLite_2.3.7
[71] utf8_1.2.4 generics_0.1.3
[73] data.table_1.16.2 rtracklayer_1.66.0
[75] httr_1.4.7 S4Arrays_1.6.0
[77] pkgconfig_2.0.3 gtable_0.3.6
[79] blob_1.2.4 XVector_0.46.0
[81] htmltools_0.5.8.1 bookdown_0.41
[83] ProtGenerics_1.38.0 scales_1.3.0
[85] alabaster.matrix_1.6.0 png_0.1-8
[87] knitr_1.48 rjson_0.2.23
[89] curl_5.2.3 cachem_1.1.0
[91] rhdf5_2.50.0 stringr_1.5.1
[93] BiocVersion_3.20.0 KernSmooth_2.23-24
[95] parallel_4.4.1 vipor_0.4.7
[97] restfulr_0.0.15 pillar_1.9.0
[99] grid_4.4.1 alabaster.schemas_1.6.0
[101] vctrs_0.6.5 origami_1.0.7
[103] BiocSingular_1.22.0 beachmat_2.22.0
[105] cluster_2.1.6 beeswarm_0.4.0
[107] evaluate_1.0.1 cli_3.6.3
[109] locfit_1.5-9.10 compiler_4.4.1
[111] Rsamtools_2.22.0 rlang_1.1.4
[113] crayon_1.5.3 future.apply_1.11.3
[115] labeling_0.4.3 ggbeeswarm_0.7.2
[117] stringi_1.8.4 alabaster.se_1.6.0
[119] viridisLite_0.4.2 BiocParallel_1.40.0
[121] assertthat_0.2.1 munsell_0.5.1
[123] Biostrings_2.74.0 coop_0.6-3
[125] lazyeval_0.2.2 Matrix_1.7-1
[127] dir.expiry_1.14.0 ExperimentHub_2.14.0
[129] future_1.34.0 sparseMatrixStats_1.18.0
[131] bit64_4.5.2 Rhdf5lib_1.28.0
[133] KEGGREST_1.46.0 statmod_1.5.0
[135] alabaster.ranges_1.6.0 highr_0.11
[137] kernlab_0.9-33 rbibutils_2.3
[139] igraph_2.1.1 memoise_2.0.1
[141] bslib_0.8.0 bit_4.5.0 References
Messmer, T., F. von Meyenn, A. Savino, F. Santos, H. Mohammed, A. T. L. Lun, J. C. Marioni, and W. Reik. 2019. “Transcriptional heterogeneity in naive and primed human pluripotent stem cells at single-cell resolution.” Cell Rep 26 (4): 815–24.