Chapter 8 Human pancreas (multiple technologies)
8.1 Introduction
For a period in 2016, there was a great deal of interest in using scRNA-seq to profile the human pancreas at cellular resolution (Muraro et al. 2016; Grun et al. 2016; Lawlor et al. 2017; Segerstolpe et al. 2016). As a consequence, we have a surplus of human pancreas datasets generated by different authors with different technologies, which provides an ideal use case for demonstrating more complex data integration strategies. This represents a more challenging application than the PBMC dataset in Chapter 1 as it involves different sequencing protocols and different patients, most likely with differences in cell type composition.
8.2 The good
We start by considering only two datasets from Muraro et al. (2016) and Grun et al. (2016). This is a relatively simple scenario involving very similar protocols (CEL-seq and CEL-seq2) and a similar set of authors.
#--- loading ---#
library(scRNAseq)
sce.grun <- GrunPancreasData()
#--- gene-annotation ---#
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 ---#
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]
#--- normalization ---#
library(scran)
set.seed(1000) # for irlba.
clusters <- quickCluster(sce.grun)
sce.grun <- computeSumFactors(sce.grun, clusters=clusters)
sce.grun <- logNormCounts(sce.grun)
#--- variance-modelling ---#
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)
## class: SingleCellExperiment
## dim: 17196 1063
## metadata(0):
## assays(2): counts logcounts
## rownames(17196): ENSG00000268895 ENSG00000121410 ... ENSG00000074755
## ENSG00000036549
## rowData names(2): symbol chr
## colnames(1063): D2ex_1 D2ex_2 ... D17TGFB_94 D17TGFB_95
## colData names(3): donor sample sizeFactor
## reducedDimNames(0):
## mainExpName: endogenous
## altExpNames(1): ERCC
#--- loading ---#
library(scRNAseq)
sce.muraro <- MuraroPancreasData()
#--- gene-annotation ---#
library(AnnotationHub)
edb <- AnnotationHub()[["AH73881"]]
gene.symb <- sub("__chr.*$", "", rownames(sce.muraro))
gene.ids <- mapIds(edb, keys=gene.symb,
keytype="SYMBOL", column="GENEID")
# Removing duplicated genes or genes without Ensembl IDs.
keep <- !is.na(gene.ids) & !duplicated(gene.ids)
sce.muraro <- sce.muraro[keep,]
rownames(sce.muraro) <- gene.ids[keep]
#--- quality-control ---#
library(scater)
stats <- perCellQCMetrics(sce.muraro)
qc <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent",
batch=sce.muraro$donor, subset=sce.muraro$donor!="D28")
sce.muraro <- sce.muraro[,!qc$discard]
#--- normalization ---#
library(scran)
set.seed(1000)
clusters <- quickCluster(sce.muraro)
sce.muraro <- computeSumFactors(sce.muraro, clusters=clusters)
sce.muraro <- logNormCounts(sce.muraro)
#--- variance-modelling ---#
block <- paste0(sce.muraro$plate, "_", sce.muraro$donor)
dec.muraro <- modelGeneVarWithSpikes(sce.muraro, "ERCC", block=block)
top.muraro <- getTopHVGs(dec.muraro, prop=0.1)
## class: SingleCellExperiment
## dim: 16940 2299
## metadata(0):
## assays(2): counts logcounts
## rownames(16940): ENSG00000268895 ENSG00000121410 ... ENSG00000159840
## ENSG00000074755
## rowData names(2): symbol chr
## colnames(2299): D28-1_1 D28-1_2 ... D30-8_93 D30-8_94
## colData names(4): label donor plate sizeFactor
## reducedDimNames(0):
## mainExpName: endogenous
## altExpNames(1): ERCC
We subset both batches to their common universe of genes; adjust their scaling to equalize sequencing coverage (not really necessary in this case, as the coverage is already similar, but we will do so anyway for consistency); and select those genes with positive average biological components for further use.
universe <- intersect(rownames(sce.grun), rownames(sce.muraro))
sce.grun2 <- sce.grun[universe,]
dec.grun2 <- dec.grun[universe,]
sce.muraro2 <- sce.muraro[universe,]
dec.muraro2 <- dec.muraro[universe,]
library(batchelor)
normed.pancreas <- multiBatchNorm(sce.grun2, sce.muraro2)
sce.grun2 <- normed.pancreas[[1]]
sce.muraro2 <- normed.pancreas[[2]]
library(scran)
combined.pan <- combineVar(dec.grun2, dec.muraro2)
chosen.genes <- rownames(combined.pan)[combined.pan$bio > 0]
We observe that rescaleBatches()
is unable to align cells from different batches in Figure 8.1.
This is attributable to differences in population composition between batches, with additional complications from non-linearities in the batch effect, e.g., when the magnitude or direction of the batch effect differs between cell types.
library(scater)
rescaled.pancreas <- rescaleBatches(sce.grun2, sce.muraro2)
set.seed(100101)
rescaled.pancreas <- runPCA(rescaled.pancreas, subset_row=chosen.genes,
exprs_values="corrected")
rescaled.pancreas <- runTSNE(rescaled.pancreas, dimred="PCA")
plotTSNE(rescaled.pancreas, colour_by="batch")
Here, we use fastMNN()
to merge together the two human pancreas datasets described earlier.
Clustering on the merged datasets yields fewer batch-specific clusters, which is recapitulated as greater intermingling between batches in Figure 8.2.
This improvement over Figure 8.1 represents the ability of fastMNN()
to adapt to more complex situations involving differences in population composition between batches.
set.seed(1011011)
mnn.pancreas <- fastMNN(sce.grun2, sce.muraro2, subset.row=chosen.genes)
snn.gr <- buildSNNGraph(mnn.pancreas, use.dimred="corrected")
clusters <- igraph::cluster_walktrap(snn.gr)$membership
tab <- table(Cluster=clusters, Batch=mnn.pancreas$batch)
tab
## Batch
## Cluster 1 2
## 1 241 279
## 2 324 256
## 3 205 842
## 4 57 194
## 5 24 108
## 6 107 392
## 7 63 78
## 8 19 0
## 9 18 114
## 10 0 17
## 11 5 19
8.3 The bad
Flushed with our previous success, we now attempt to merge the other datasets from Lawlor et al. (2017) and Segerstolpe et al. (2016). This is a more challenging task as it involves different technologies, mixtures of UMI and read count data and a more diverse set of authors (presumably with greater differences in the patient population).
#--- loading ---#
library(scRNAseq)
sce.lawlor <- LawlorPancreasData()
#--- gene-annotation ---#
library(AnnotationHub)
edb <- AnnotationHub()[["AH73881"]]
anno <- select(edb, keys=rownames(sce.lawlor), keytype="GENEID",
columns=c("SYMBOL", "SEQNAME"))
rowData(sce.lawlor) <- anno[match(rownames(sce.lawlor), anno[,1]),-1]
#--- quality-control ---#
library(scater)
stats <- perCellQCMetrics(sce.lawlor,
subsets=list(Mito=which(rowData(sce.lawlor)$SEQNAME=="MT")))
qc <- quickPerCellQC(stats, percent_subsets="subsets_Mito_percent",
batch=sce.lawlor$`islet unos id`)
sce.lawlor <- sce.lawlor[,!qc$discard]
#--- normalization ---#
library(scran)
set.seed(1000)
clusters <- quickCluster(sce.lawlor)
sce.lawlor <- computeSumFactors(sce.lawlor, clusters=clusters)
sce.lawlor <- logNormCounts(sce.lawlor)
#--- variance-modelling ---#
dec.lawlor <- modelGeneVar(sce.lawlor, block=sce.lawlor$`islet unos id`)
chosen.genes <- getTopHVGs(dec.lawlor, n=2000)
## class: SingleCellExperiment
## dim: 26616 604
## metadata(0):
## assays(2): counts logcounts
## rownames(26616): ENSG00000229483 ENSG00000232849 ... ENSG00000251576
## ENSG00000082898
## rowData names(2): SYMBOL SEQNAME
## colnames(604): 10th_C11_S96 10th_C13_S61 ... 9th-C96_S81 9th-C9_S13
## colData names(9): title age ... Sex sizeFactor
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
#--- loading ---#
library(scRNAseq)
sce.seger <- SegerstolpePancreasData()
#--- gene-annotation ---#
library(AnnotationHub)
edb <- AnnotationHub()[["AH73881"]]
symbols <- rowData(sce.seger)$symbol
ens.id <- mapIds(edb, keys=symbols, keytype="SYMBOL", column="GENEID")
ens.id <- ifelse(is.na(ens.id), symbols, ens.id)
# Removing duplicated rows.
keep <- !duplicated(ens.id)
sce.seger <- sce.seger[keep,]
rownames(sce.seger) <- ens.id[keep]
#--- sample-annotation ---#
emtab.meta <- colData(sce.seger)[,c("cell type", "disease",
"individual", "single cell well quality")]
colnames(emtab.meta) <- c("CellType", "Disease", "Donor", "Quality")
colData(sce.seger) <- emtab.meta
sce.seger$CellType <- gsub(" cell", "", sce.seger$CellType)
sce.seger$CellType <- paste0(
toupper(substr(sce.seger$CellType, 1, 1)),
substring(sce.seger$CellType, 2))
#--- quality-control ---#
low.qual <- sce.seger$Quality == "low quality cell"
library(scater)
stats <- perCellQCMetrics(sce.seger)
qc <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent",
batch=sce.seger$Donor,
subset=!sce.seger$Donor %in% c("HP1504901", "HP1509101"))
sce.seger <- sce.seger[,!(qc$discard | low.qual)]
#--- normalization ---#
library(scran)
clusters <- quickCluster(sce.seger)
sce.seger <- computeSumFactors(sce.seger, clusters=clusters)
sce.seger <- logNormCounts(sce.seger)
#--- variance-modelling ---#
for.hvg <- sce.seger[,librarySizeFactors(altExp(sce.seger)) > 0 & sce.seger$Donor!="AZ"]
dec.seger <- modelGeneVarWithSpikes(for.hvg, "ERCC", block=for.hvg$Donor)
chosen.hvgs <- getTopHVGs(dec.seger, n=2000)
## class: SingleCellExperiment
## dim: 25454 2090
## metadata(0):
## assays(2): counts logcounts
## rownames(25454): ENSG00000118473 ENSG00000142920 ... ENSG00000278306
## eGFP
## rowData names(2): symbol refseq
## colnames(2090): HP1502401_H13 HP1502401_J14 ... HP1526901T2D_N8
## HP1526901T2D_A8
## colData names(5): CellType Disease Donor Quality sizeFactor
## reducedDimNames(0):
## mainExpName: endogenous
## altExpNames(1): ERCC
We perform the usual routine to obtain re-normalized values and a set of HVGs. Here, we put all the objects into a list to avoid having to explicitly type their names separately.
all.sce <- list(Grun=sce.grun, Muraro=sce.muraro,
Lawlor=sce.lawlor, Seger=sce.seger)
all.dec <- list(Grun=dec.grun, Muraro=dec.muraro,
Lawlor=dec.lawlor, Seger=dec.seger)
universe <- Reduce(intersect, lapply(all.sce, rownames))
all.sce <- lapply(all.sce, "[", i=universe,)
all.dec <- lapply(all.dec, "[", i=universe,)
normed.pancreas <- do.call(multiBatchNorm, all.sce)
combined.pan <- do.call(combineVar, all.dec)
chosen.genes <- rownames(combined.pan)[combined.pan$bio > 0]
We observe that the merge is generally successful, with many clusters containing contributions from each batch (Figure 8.3). There are few clusters that are specific to the Segerstolpe dataset, and if we were naive, we might consider them to represent interesting subpopulations that are not present in the other datasets.
set.seed(1011110)
mnn.pancreas <- fastMNN(normed.pancreas)
# Bumping up 'k' to get broader clusters for this demonstration.
snn.gr <- buildSNNGraph(mnn.pancreas, use.dimred="corrected", k=20)
clusters <- igraph::cluster_walktrap(snn.gr)$membership
clusters <- factor(clusters)
tab <- table(Cluster=clusters, Batch=mnn.pancreas$batch)
tab
## Batch
## Cluster Grun Lawlor Muraro Seger
## 1 107 245 387 188
## 2 308 28 259 386
## 3 0 0 0 213
## 4 0 0 0 42
## 5 35 0 1 0
## 6 202 11 224 128
## 7 19 18 118 172
## 8 40 11 53 57
## 9 128 201 181 119
## 10 57 17 194 105
## 11 66 20 83 17
## 12 71 24 651 200
## 13 24 18 107 55
## 14 0 1 16 4
## 15 0 0 0 112
## 16 0 0 1 196
## 17 1 2 5 11
## 18 0 0 0 26
## 19 0 0 0 42
## 20 5 8 19 17
mnn.pancreas <- runTSNE(mnn.pancreas, dimred="corrected")
gridExtra::grid.arrange(
plotTSNE(mnn.pancreas, colour_by="batch", text_by=I(clusters)),
plotTSNE(mnn.pancreas, colour_by=I(clusters), text_by=I(clusters)),
ncol=2
)
Fortunately, we are battle-hardened and cynical, so we are sure to check for other sources of variation. The most obvious candidate is the donor of origin for each cell (Figure 8.4), which correlates strongly to these Segerstolpe-only clusters. This is not surprising given the large differences between humans in the wild, but donor-level variation is not interesting for the purposes of cell type characterization. (That said, preservation of the within-dataset donor effects is the technically correct course of action here, as a batch correction method should try to avoid removing heterogeneity within each of its defined batches.)
donors <- c(
normed.pancreas$Grun$donor,
normed.pancreas$Muraro$donor,
normed.pancreas$Lawlor$`islet unos id`,
normed.pancreas$Seger$Donor
)
seger.donors <- donors
seger.donors[mnn.pancreas$batch!="Seger"] <- NA
plotTSNE(mnn.pancreas, colour_by=I(seger.donors))
8.4 The ugly
Given these results, the most prudent course of action is to remove the donor effects within each dataset in addition to the batch effects across datasets.
This involves a bit more work to properly specify the two levels of unwanted heterogeneity.
To make our job a bit easier, we use the noCorrect()
utility to combine all batches into a single SingleCellExperiment
object.
We then call fastMNN()
on the combined object with our chosen HVGs, using the batch=
argument to specify which cells belong to which donors.
This will progressively merge cells from each donor in each batch until all cells are mapped onto a common coordinate space.
For some extra sophistication, we also set the weights=
argument to ensure that each batch contributes equally to the PCA, regardless of the number of donors present in that batch; see ?multiBatchPCA
for more details.
donors.per.batch <- split(combined$donor, combined$batch)
donors.per.batch <- lapply(donors.per.batch, unique)
donors.per.batch
## $Grun
## [1] "D2" "D3" "D7" "D10" "D17"
##
## $Lawlor
## [1] "ACIW009" "ACJV399" "ACCG268" "ACCR015A" "ACEK420A" "ACEL337" "ACHY057"
## [8] "ACIB065"
##
## $Muraro
## [1] "D28" "D29" "D31" "D30"
##
## $Seger
## [1] "HP1502401" "HP1504101T2D" "AZ" "HP1508501T2D" "HP1506401"
## [6] "HP1507101" "HP1509101" "HP1504901" "HP1525301T2D" "HP1526901T2D"
set.seed(1010100)
multiout <- fastMNN(combined, batch=combined$donor,
subset.row=chosen.genes, weights=donors.per.batch)
# Renaming metadata fields for easier communication later.
multiout$dataset <- combined$batch
multiout$donor <- multiout$batch
multiout$batch <- NULL
With this approach, we see that the Segerstolpe-only clusters have disappeared (Figure 8.5). Visually, there also seems to be much greater mixing between cells from different Segerstolpe donors. This suggests that we have removed most of the donor effect, which simplifies the interpretation of our clusters.
library(scater)
g <- buildSNNGraph(multiout, use.dimred=1, k=20)
clusters <- igraph::cluster_walktrap(g)$membership
tab <- table(clusters, multiout$dataset)
tab
##
## clusters Grun Lawlor Muraro Seger
## 1 248 20 278 186
## 2 337 27 259 388
## 3 172 256 473 292
## 4 200 237 835 865
## 5 58 17 193 108
## 6 24 18 107 55
## 7 5 9 19 17
## 8 0 1 17 4
## 9 19 19 118 175
multiout <- runTSNE(multiout, dimred="corrected")
gridExtra::grid.arrange(
plotTSNE(multiout, colour_by="dataset", text_by=I(clusters)),
plotTSNE(multiout, colour_by=I(seger.donors)),
ncol=2
)
Our clusters compare well to the published annotations, indicating that we did not inadvertently discard important factors of variation during correction. (Though in this case, the cell types are so well defined, it would be quite a feat to fail to separate them!)
proposed <- c(rep(NA, ncol(sce.grun)),
sce.muraro$label,
sce.lawlor$`cell type`,
sce.seger$CellType)
proposed <- tolower(proposed)
proposed[proposed=="gamma/pp"] <- "gamma"
proposed[proposed=="pp"] <- "gamma"
proposed[proposed=="duct"] <- "ductal"
proposed[proposed=="psc"] <- "stellate"
table(proposed, clusters)
## clusters
## proposed 1 2 3 4 5 6 7 8 9
## acinar 421 2 0 1 0 0 1 0 0
## alpha 1 7 24 1851 1 0 1 0 2
## beta 3 4 919 3 0 2 1 1 7
## co-expression 0 0 38 1 0 0 0 0 0
## delta 0 2 5 1 306 0 1 0 1
## ductal 6 614 4 0 0 5 0 10 1
## endothelial 0 0 0 0 0 1 33 0 0
## epsilon 0 0 1 0 0 0 0 0 7
## gamma 2 0 0 0 0 0 0 0 280
## mesenchymal 0 1 0 0 0 79 0 0 0
## mhc class ii 0 0 0 0 0 0 0 4 0
## nana 2 8 2 12 1 0 1 0 0
## none/other 0 3 1 2 0 0 4 1 1
## stellate 0 0 0 0 0 71 1 0 0
## unclassified 0 0 0 0 0 2 0 0 0
## unclassified endocrine 0 0 3 3 0 0 0 0 0
## unclear 0 4 0 0 0 0 0 0 0
Session Info
R version 4.3.2 Patched (2023-11-13 r85521)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS
Matrix products: default
BLAS: /home/biocbuild/bbs-3.18-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] bluster_1.12.0 scater_1.30.1
[3] ggplot2_3.4.4 scran_1.30.2
[5] scuttle_1.12.0 batchelor_1.18.1
[7] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
[9] Biobase_2.62.0 GenomicRanges_1.54.1
[11] GenomeInfoDb_1.38.6 IRanges_2.36.0
[13] S4Vectors_0.40.2 BiocGenerics_0.48.1
[15] MatrixGenerics_1.14.0 matrixStats_1.2.0
[17] BiocStyle_2.30.0 rebook_1.12.0
loaded via a namespace (and not attached):
[1] bitops_1.0-7 gridExtra_2.3
[3] CodeDepends_0.6.5 rlang_1.1.3
[5] magrittr_2.0.3 compiler_4.3.2
[7] dir.expiry_1.10.0 DelayedMatrixStats_1.24.0
[9] vctrs_0.6.5 pkgconfig_2.0.3
[11] crayon_1.5.2 fastmap_1.1.1
[13] XVector_0.42.0 labeling_0.4.3
[15] utf8_1.2.4 rmarkdown_2.25
[17] graph_1.80.0 ggbeeswarm_0.7.2
[19] xfun_0.42 zlibbioc_1.48.0
[21] cachem_1.0.8 beachmat_2.18.1
[23] jsonlite_1.8.8 highr_0.10
[25] DelayedArray_0.28.0 BiocParallel_1.36.0
[27] irlba_2.3.5.1 parallel_4.3.2
[29] cluster_2.1.6 R6_2.5.1
[31] bslib_0.6.1 limma_3.58.1
[33] jquerylib_0.1.4 Rcpp_1.0.12
[35] bookdown_0.37 knitr_1.45
[37] Matrix_1.6-5 igraph_2.0.1.1
[39] tidyselect_1.2.0 abind_1.4-5
[41] yaml_2.3.8 viridis_0.6.5
[43] codetools_0.2-19 lattice_0.22-5
[45] tibble_3.2.1 withr_3.0.0
[47] evaluate_0.23 Rtsne_0.17
[49] pillar_1.9.0 BiocManager_1.30.22
[51] filelock_1.0.3 generics_0.1.3
[53] RCurl_1.98-1.14 sparseMatrixStats_1.14.0
[55] munsell_0.5.0 scales_1.3.0
[57] glue_1.7.0 metapod_1.10.1
[59] tools_4.3.2 BiocNeighbors_1.20.2
[61] ScaledMatrix_1.10.0 locfit_1.5-9.8
[63] XML_3.99-0.16.1 cowplot_1.1.3
[65] grid_4.3.2 edgeR_4.0.15
[67] colorspace_2.1-0 GenomeInfoDbData_1.2.11
[69] beeswarm_0.4.0 BiocSingular_1.18.0
[71] vipor_0.4.7 cli_3.6.2
[73] rsvd_1.0.5 rappdirs_0.3.3
[75] fansi_1.0.6 S4Arrays_1.2.0
[77] viridisLite_0.4.2 dplyr_1.1.4
[79] ResidualMatrix_1.12.0 gtable_0.3.4
[81] sass_0.4.8 digest_0.6.34
[83] SparseArray_1.2.4 ggrepel_0.9.5
[85] dqrng_0.3.2 farver_2.1.1
[87] htmltools_0.5.7 lifecycle_1.0.4
[89] statmod_1.5.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.
Lawlor, N., J. George, M. Bolisetty, R. Kursawe, L. Sun, V. Sivakamasundari, I. Kycia, P. Robson, and M. L. Stitzel. 2017. “Single-cell transcriptomes identify human islet cell signatures and reveal cell-type-specific expression changes in type 2 diabetes.” Genome Res. 27 (2): 208–22.
Muraro, M. J., G. Dharmadhikari, D. Grun, N. Groen, T. Dielen, E. Jansen, L. van Gurp, et al. 2016. “A Single-Cell Transcriptome Atlas of the Human Pancreas.” Cell Syst 3 (4): 385–94.
Segerstolpe, A., A. Palasantza, P. Eliasson, E. M. Andersson, A. C. Andreasson, X. Sun, S. Picelli, et al. 2016. “Single-Cell Transcriptome Profiling of Human Pancreatic Islets in Health and Type 2 Diabetes.” Cell Metab. 24 (4): 593–607.