Chapter 9 Grun mouse HSC (CEL-seq)
9.1 Introduction
This performs an analysis of the mouse haematopoietic stem cell (HSC) dataset generated with CEL-seq (Grun et al. 2016). Despite its name, this dataset actually contains both sorted HSCs and a population of micro-dissected bone marrow cells.
9.2 Data loading
library(AnnotationHub)
ens.mm.v97 <- AnnotationHub()[["AH73905"]]
anno <- select(ens.mm.v97, keys=rownames(sce.grun.hsc),
keytype="GENEID", columns=c("SYMBOL", "SEQNAME"))
rowData(sce.grun.hsc) <- anno[match(rownames(sce.grun.hsc), anno$GENEID),]After loading and annotation, we inspect the resulting SingleCellExperiment object:
## class: SingleCellExperiment
## dim: 21817 1915
## metadata(0):
## assays(1): counts
## rownames(21817): ENSMUSG00000109644 ENSMUSG00000007777 ...
## ENSMUSG00000055670 ENSMUSG00000039068
## rowData names(3): GENEID SYMBOL SEQNAME
## colnames(1915): JC4_349_HSC_FE_S13_ JC4_350_HSC_FE_S13_ ...
## JC48P6_1203_HSC_FE_S8_ JC48P6_1204_HSC_FE_S8_
## colData names(2): sample protocol
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):9.3 Quality control
For some reason, no mitochondrial transcripts are available, and we have no spike-in transcripts, so we only use the number of detected genes and the library size for quality control. We block on the protocol used for cell extraction, ignoring the micro-dissected cells when computing this threshold. This is based on our judgement that a majority of micro-dissected plates consist of a majority of low-quality cells, compromising the assumptions of outlier detection.
library(scuttle)
stats <- perCellQCMetrics(sce.grun.hsc)
qc <- quickPerCellQC(stats, batch=sce.grun.hsc$protocol,
subset=grepl("sorted", sce.grun.hsc$protocol))
sce.grun.hsc <- sce.grun.hsc[,!qc$discard]We examine the number of cells discarded for each reason.
## low_lib_size low_n_features discard
## 465 482 488We create some diagnostic plots for each metric (Figure 9.1). The library sizes are unusually low for many plates of micro-dissected cells; this may be attributable to damage induced by the extraction protocol compared to cell sorting.
colData(unfiltered) <- cbind(colData(unfiltered), stats)
unfiltered$discard <- qc$discard
library(scater)
gridExtra::grid.arrange(
plotColData(unfiltered, y="sum", x="sample", colour_by="discard",
other_fields="protocol") + scale_y_log10() + ggtitle("Total count") +
facet_wrap(~protocol),
plotColData(unfiltered, y="detected", x="sample", colour_by="discard",
other_fields="protocol") + scale_y_log10() +
ggtitle("Detected features") + facet_wrap(~protocol),
ncol=1
)
Figure 9.1: Distribution of each QC metric across cells in the Grun HSC dataset. Each point represents a cell and is colored according to whether that cell was discarded.
9.4 Normalization
library(scran)
set.seed(101000110)
clusters <- quickCluster(sce.grun.hsc)
sce.grun.hsc <- computeSumFactors(sce.grun.hsc, clusters=clusters)
sce.grun.hsc <- logNormCounts(sce.grun.hsc)We examine some key metrics for the distribution of size factors, and compare it to the library sizes as a sanity check (Figure 9.2).
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.027 0.290 0.603 1.000 1.201 16.433plot(librarySizeFactors(sce.grun.hsc), sizeFactors(sce.grun.hsc), pch=16,
xlab="Library size factors", ylab="Deconvolution factors", log="xy")
Figure 9.2: Relationship between the library size factors and the deconvolution size factors in the Grun HSC dataset.
9.5 Variance modelling
We create a mean-variance trend based on the expectation that UMI counts have Poisson technical noise. We do not block on sample here as we want to preserve any difference between the micro-dissected cells and the sorted HSCs.
set.seed(00010101)
dec.grun.hsc <- modelGeneVarByPoisson(sce.grun.hsc)
top.grun.hsc <- getTopHVGs(dec.grun.hsc, prop=0.1)The lack of a typical “bump” shape in Figure 9.3 is caused by the low counts.
plot(dec.grun.hsc$mean, dec.grun.hsc$total, pch=16, cex=0.5,
xlab="Mean of log-expression", ylab="Variance of log-expression")
curfit <- metadata(dec.grun.hsc)
curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
Figure 9.3: Per-gene variance as a function of the mean for the log-expression values in the Grun HSC dataset. Each point represents a gene (black) with the mean-variance trend (blue) fitted to the simulated Poisson-distributed noise.
9.6 Dimensionality reduction
set.seed(101010011)
sce.grun.hsc <- denoisePCA(sce.grun.hsc, technical=dec.grun.hsc, subset.row=top.grun.hsc)
sce.grun.hsc <- runTSNE(sce.grun.hsc, dimred="PCA")We check that the number of retained PCs is sensible.
## [1] 99.7 Clustering
snn.gr <- buildSNNGraph(sce.grun.hsc, use.dimred="PCA")
colLabels(sce.grun.hsc) <- factor(igraph::cluster_walktrap(snn.gr)$membership)##
## 1 2 3 4 5 6 7 8 9 10 11 12
## 259 148 221 103 177 108 48 122 98 63 62 18short <- ifelse(grepl("micro", sce.grun.hsc$protocol), "micro", "sorted")
gridExtra:::grid.arrange(
plotTSNE(sce.grun.hsc, colour_by="label"),
plotTSNE(sce.grun.hsc, colour_by=I(short)),
ncol=2
)
Figure 9.4: Obligatory \(t\)-SNE plot of the Grun HSC dataset, where each point represents a cell and is colored according to the assigned cluster (left) or extraction protocol (right).
9.8 Marker gene detection
markers <- findMarkers(sce.grun.hsc, test.type="wilcox", direction="up",
row.data=rowData(sce.grun.hsc)[,"SYMBOL",drop=FALSE])To illustrate the manual annotation process, we examine the marker genes for one of the clusters. Upregulation of Camp, Lcn2, Ltf and lysozyme genes indicates that this cluster contains cells of neuronal origin.
chosen <- markers[['6']]
best <- chosen[chosen$Top <= 10,]
aucs <- getMarkerEffects(best, prefix="AUC")
rownames(aucs) <- best$SYMBOL
library(pheatmap)
pheatmap(aucs, color=viridis::plasma(100))
Figure 9.5: Heatmap of the AUCs for the top marker genes in cluster 6 compared to all other clusters in the Grun HSC dataset.
Session Info
R Under development (unstable) (2024-10-21 r87258)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.1 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
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] pheatmap_1.0.12 scran_1.35.0
[3] scater_1.35.0 ggplot2_3.5.1
[5] scuttle_1.17.0 AnnotationHub_3.15.0
[7] BiocFileCache_2.15.0 dbplyr_2.5.0
[9] ensembldb_2.31.0 AnnotationFilter_1.31.0
[11] GenomicFeatures_1.59.1 AnnotationDbi_1.69.0
[13] scRNAseq_2.21.0 SingleCellExperiment_1.29.1
[15] SummarizedExperiment_1.37.0 Biobase_2.67.0
[17] GenomicRanges_1.59.0 GenomeInfoDb_1.43.0
[19] IRanges_2.41.0 S4Vectors_0.45.1
[21] BiocGenerics_0.53.1 generics_0.1.3
[23] MatrixGenerics_1.19.0 matrixStats_1.4.1
[25] BiocStyle_2.35.0 rebook_1.17.0
loaded via a namespace (and not attached):
[1] RColorBrewer_1.1-3 jsonlite_1.8.9 CodeDepends_0.6.6
[4] magrittr_2.0.3 ggbeeswarm_0.7.2 gypsum_1.3.0
[7] farver_2.1.2 rmarkdown_2.29 BiocIO_1.17.0
[10] zlibbioc_1.53.0 vctrs_0.6.5 memoise_2.0.1
[13] Rsamtools_2.23.0 RCurl_1.98-1.16 htmltools_0.5.8.1
[16] S4Arrays_1.7.1 curl_6.0.0 BiocNeighbors_2.1.0
[19] Rhdf5lib_1.29.0 SparseArray_1.7.1 rhdf5_2.51.0
[22] sass_0.4.9 alabaster.base_1.7.1 bslib_0.8.0
[25] alabaster.sce_1.7.0 httr2_1.0.6 cachem_1.1.0
[28] GenomicAlignments_1.43.0 igraph_2.1.1 mime_0.12
[31] lifecycle_1.0.4 pkgconfig_2.0.3 rsvd_1.0.5
[34] Matrix_1.7-1 R6_2.5.1 fastmap_1.2.0
[37] GenomeInfoDbData_1.2.13 digest_0.6.37 colorspace_2.1-1
[40] dqrng_0.4.1 irlba_2.3.5.1 ExperimentHub_2.15.0
[43] RSQLite_2.3.7 beachmat_2.23.0 labeling_0.4.3
[46] filelock_1.0.3 fansi_1.0.6 httr_1.4.7
[49] abind_1.4-8 compiler_4.5.0 bit64_4.5.2
[52] withr_3.0.2 BiocParallel_1.41.0 viridis_0.6.5
[55] DBI_1.2.3 HDF5Array_1.35.1 alabaster.ranges_1.7.0
[58] alabaster.schemas_1.7.0 rappdirs_0.3.3 DelayedArray_0.33.1
[61] bluster_1.17.0 rjson_0.2.23 tools_4.5.0
[64] vipor_0.4.7 beeswarm_0.4.0 glue_1.8.0
[67] restfulr_0.0.15 rhdf5filters_1.19.0 grid_4.5.0
[70] Rtsne_0.17 cluster_2.1.6 gtable_0.3.6
[73] metapod_1.15.0 BiocSingular_1.23.0 ScaledMatrix_1.15.0
[76] utf8_1.2.4 XVector_0.47.0 ggrepel_0.9.6
[79] BiocVersion_3.21.1 pillar_1.9.0 limma_3.63.2
[82] dplyr_1.1.4 lattice_0.22-6 rtracklayer_1.67.0
[85] bit_4.5.0 tidyselect_1.2.1 locfit_1.5-9.10
[88] Biostrings_2.75.1 knitr_1.49 gridExtra_2.3
[91] bookdown_0.41 ProtGenerics_1.39.0 edgeR_4.5.0
[94] xfun_0.49 statmod_1.5.0 UCSC.utils_1.3.0
[97] lazyeval_0.2.2 yaml_2.3.10 evaluate_1.0.1
[100] codetools_0.2-20 tibble_3.2.1 alabaster.matrix_1.7.0
[103] BiocManager_1.30.25 graph_1.85.0 cli_3.6.3
[106] munsell_0.5.1 jquerylib_0.1.4 Rcpp_1.0.13-1
[109] dir.expiry_1.15.0 png_0.1-8 XML_3.99-0.17
[112] parallel_4.5.0 blob_1.2.4 bitops_1.0-9
[115] viridisLite_0.4.2 alabaster.se_1.7.0 scales_1.3.0
[118] purrr_1.0.2 crayon_1.5.3 rlang_1.1.4
[121] cowplot_1.1.3 KEGGREST_1.47.0