Chapter 8 Segerstolpe human pancreas (Smart-seq2)
8.1 Introduction
This performs an analysis of the Segerstolpe et al. (2016) dataset, consisting of human pancreas cells from various donors.
8.2 Data loading
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]
We simplify the names of some of the relevant column metadata fields for ease of access. Some editing of the cell type labels is necessary for consistency with other data sets.
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))
8.3 Quality control
We remove low quality cells that were marked by the authors.
We then perform additional quality control as some of the remaining cells still have very low counts and numbers of detected features.
For some batches that seem to have a majority of low-quality cells (Figure 8.1), we use the other batches to define an appropriate threshold via subset=
.
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)]
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") +
theme(axis.text.x = element_text(angle = 90)),
plotColData(unfiltered, x="Donor", y="detected", colour_by="discard") +
scale_y_log10() + ggtitle("Detected features") +
theme(axis.text.x = element_text(angle = 90)),
plotColData(unfiltered, x="Donor", y="altexps_ERCC_percent",
colour_by="discard") + ggtitle("ERCC percent") +
theme(axis.text.x = element_text(angle = 90)),
ncol=2
)
## low_lib_size low_n_features high_altexps_ERCC_percent
## 788 1056 1031
## discard
## 1246
8.4 Normalization
We don’t normalize the spike-ins at this point as there are some cells with no spike-in counts.
library(scran)
clusters <- quickCluster(sce.seger)
sce.seger <- computeSumFactors(sce.seger, clusters=clusters)
sce.seger <- logNormCounts(sce.seger)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.014 0.390 0.708 1.000 1.332 11.182
plot(librarySizeFactors(sce.seger), sizeFactors(sce.seger), pch=16,
xlab="Library size factors", ylab="Deconvolution factors", log="xy")
8.5 Variance modelling
We do not use cells with no spike-ins for variance modelling. Donor AZ also has very low spike-in counts and is subsequently ignored.
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)
par(mfrow=c(3,3))
blocked.stats <- dec.seger$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)
}
8.6 Dimensionality reduction
We pick the first 25 PCs for downstream analyses, as it’s a nice square number.
8.7 Clustering
library(bluster)
clust.out <- clusterRows(reducedDim(sce.seger, "PCA"), NNGraphParam(), full=TRUE)
snn.gr <- clust.out$objects$graph
colLabels(sce.seger) <- clust.out$clusters
We see a strong donor effect in Figures 8.4 and 5.3.
This might be due to differences in cell type composition between donors,
but the more likely explanation is that of a technical difference in plate processing or uninteresting genotypic differences.
The implication is that we should have called fastMNN()
at some point.
tab <- table(Cluster=colLabels(sce.seger), Donor=sce.seger$Donor)
library(pheatmap)
pheatmap(log10(tab+10), color=viridis::viridis(100))
gridExtra::grid.arrange(
plotTSNE(sce.seger, colour_by="label"),
plotTSNE(sce.seger, colour_by="Donor"),
ncol=2
)
8.8 Data integration
We repeat the clustering after running fastMNN()
on the donors.
This yields a more coherent set of clusters in Figure 8.6 where each cluster contains contributions from all donors.
library(batchelor)
set.seed(10001010)
corrected <- fastMNN(sce.seger, batch=sce.seger$Donor, subset.row=chosen.hvgs)
set.seed(10000001)
corrected <- runTSNE(corrected, dimred="corrected")
colLabels(corrected) <- clusterRows(reducedDim(corrected, "corrected"), NNGraphParam())
tab <- table(Cluster=colLabels(corrected), Donor=corrected$batch)
tab
## Donor
## Cluster AZ HP1502401 HP1504101T2D HP1504901 HP1506401 HP1507101 HP1508501T2D
## 1 3 19 3 11 67 8 78
## 2 14 53 13 19 37 41 20
## 3 2 2 1 1 44 1 1
## 4 2 18 7 3 36 2 28
## 5 29 114 140 72 26 136 121
## 6 8 21 9 6 2 6 6
## 7 1 1 1 9 0 1 2
## 8 2 1 3 10 2 6 12
## 9 4 20 70 8 16 2 8
## Donor
## Cluster HP1509101 HP1525301T2D HP1526901T2D
## 1 27 124 46
## 2 14 11 70
## 3 0 1 4
## 4 2 23 9
## 5 49 85 96
## 6 11 5 34
## 7 2 2 1
## 8 3 13 4
## 9 1 10 34
gridExtra::grid.arrange(
plotTSNE(corrected, colour_by="label"),
plotTSNE(corrected, colour_by="batch"),
ncol=2
)
8.9 Multi-sample comparisons
This particular dataset contains both healthy donors and those with type II diabetes. It is thus of some interest to identify genes that are differentially expressed upon disease in each cell type. To keep things simple, we use the author-provided annotation rather than determining the cell type for each of our clusters.
## class: SingleCellExperiment
## dim: 25454 105
## metadata(0):
## assays(1): counts
## rownames(25454): ENSG00000118473 ENSG00000142920 ... ENSG00000278306
## eGFP
## rowData names(2): symbol refseq
## colnames: NULL
## colData names(9): CellType Disease ... CellType ncells
## reducedDimNames(2): PCA TSNE
## mainExpName: endogenous
## altExpNames(0):
Here, we will use the voom
pipeline from the limma package instead of the QL approach with edgeR.
This allows us to use sample weights to better account for the variation in the precision of each pseudo-bulk profile.
We see that insulin is downregulated in beta cells in the disease state, which is sensible enough.
summed.beta <- summed[,summed$CellType=="Beta"]
library(edgeR)
y.beta <- DGEList(counts(summed.beta), samples=colData(summed.beta),
genes=rowData(summed.beta)[,"symbol",drop=FALSE])
y.beta <- y.beta[filterByExpr(y.beta, group=y.beta$samples$Disease),]
y.beta <- calcNormFactors(y.beta)
design <- model.matrix(~Disease, y.beta$samples)
v.beta <- voomWithQualityWeights(y.beta, design)
fit.beta <- lmFit(v.beta)
fit.beta <- eBayes(fit.beta, robust=TRUE)
res.beta <- topTable(fit.beta, sort.by="p", n=Inf,
coef="Diseasetype II diabetes mellitus")
head(res.beta)
## symbol logFC AveExpr t P.Value adj.P.Val B
## ENSG00000254647 INS -2.728 16.680 -7.671 3.191e-06 0.03902 4.842
## ENSG00000137731 FXYD2 -2.595 7.265 -6.705 1.344e-05 0.08219 3.353
## ENSG00000169297 NR0B1 -2.092 6.790 -5.789 5.810e-05 0.09916 1.984
## ENSG00000181029 TRAPPC5 -2.127 7.046 -5.678 7.007e-05 0.09916 1.877
## ENSG00000105707 HPN -1.803 6.118 -5.654 7.298e-05 0.09916 1.740
## LOC284889 LOC284889 -2.113 6.652 -5.515 9.259e-05 0.09916 1.571
We also create some diagnostic plots to check for potential problems in the analysis.
The MA plots exhibit the expected shape (Figure 8.7)
while the differences in the sample weights in Figure 8.8 justify the use of voom()
in this context.
# Easier to just re-run it with plot=TRUE than
# to try to make the plot from 'v.beta'.
voomWithQualityWeights(y.beta, design, plot=TRUE)
Session Info
R version 4.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.5 LTS
Matrix products: default
BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas.so
LAPACK: /home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack.so
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
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] edgeR_3.40.0 limma_3.54.0
[3] batchelor_1.14.0 pheatmap_1.0.12
[5] bluster_1.8.0 BiocSingular_1.14.0
[7] scran_1.26.0 scater_1.26.0
[9] ggplot2_3.3.6 scuttle_1.8.0
[11] ensembldb_2.22.0 AnnotationFilter_1.22.0
[13] GenomicFeatures_1.50.0 AnnotationDbi_1.60.0
[15] AnnotationHub_3.6.0 BiocFileCache_2.6.0
[17] dbplyr_2.2.1 scRNAseq_2.11.0
[19] SingleCellExperiment_1.20.0 SummarizedExperiment_1.28.0
[21] Biobase_2.58.0 GenomicRanges_1.50.0
[23] GenomeInfoDb_1.34.0 IRanges_2.32.0
[25] S4Vectors_0.36.0 BiocGenerics_0.44.0
[27] MatrixGenerics_1.10.0 matrixStats_0.62.0
[29] BiocStyle_2.26.0 rebook_1.8.0
loaded via a namespace (and not attached):
[1] igraph_1.3.5 lazyeval_0.2.2
[3] BiocParallel_1.32.0 digest_0.6.30
[5] htmltools_0.5.3 viridis_0.6.2
[7] fansi_1.0.3 magrittr_2.0.3
[9] memoise_2.0.1 ScaledMatrix_1.6.0
[11] cluster_2.1.4 Biostrings_2.66.0
[13] prettyunits_1.1.1 colorspace_2.0-3
[15] blob_1.2.3 rappdirs_0.3.3
[17] ggrepel_0.9.1 xfun_0.34
[19] dplyr_1.0.10 crayon_1.5.2
[21] RCurl_1.98-1.9 jsonlite_1.8.3
[23] graph_1.76.0 glue_1.6.2
[25] gtable_0.3.1 zlibbioc_1.44.0
[27] XVector_0.38.0 DelayedArray_0.24.0
[29] scales_1.2.1 DBI_1.1.3
[31] Rcpp_1.0.9 viridisLite_0.4.1
[33] xtable_1.8-4 progress_1.2.2
[35] dqrng_0.3.0 bit_4.0.4
[37] rsvd_1.0.5 ResidualMatrix_1.8.0
[39] metapod_1.6.0 httr_1.4.4
[41] RColorBrewer_1.1-3 dir.expiry_1.6.0
[43] ellipsis_0.3.2 pkgconfig_2.0.3
[45] XML_3.99-0.12 farver_2.1.1
[47] CodeDepends_0.6.5 sass_0.4.2
[49] locfit_1.5-9.6 utf8_1.2.2
[51] labeling_0.4.2 tidyselect_1.2.0
[53] rlang_1.0.6 later_1.3.0
[55] munsell_0.5.0 BiocVersion_3.16.0
[57] tools_4.2.1 cachem_1.0.6
[59] cli_3.4.1 generics_0.1.3
[61] RSQLite_2.2.18 ExperimentHub_2.6.0
[63] evaluate_0.17 stringr_1.4.1
[65] fastmap_1.1.0 yaml_2.3.6
[67] knitr_1.40 bit64_4.0.5
[69] purrr_0.3.5 KEGGREST_1.38.0
[71] sparseMatrixStats_1.10.0 mime_0.12
[73] xml2_1.3.3 biomaRt_2.54.0
[75] compiler_4.2.1 beeswarm_0.4.0
[77] filelock_1.0.2 curl_4.3.3
[79] png_0.1-7 interactiveDisplayBase_1.36.0
[81] statmod_1.4.37 tibble_3.1.8
[83] bslib_0.4.0 stringi_1.7.8
[85] highr_0.9 lattice_0.20-45
[87] ProtGenerics_1.30.0 Matrix_1.5-1
[89] vctrs_0.5.0 pillar_1.8.1
[91] lifecycle_1.0.3 BiocManager_1.30.19
[93] jquerylib_0.1.4 BiocNeighbors_1.16.0
[95] cowplot_1.1.1 bitops_1.0-7
[97] irlba_2.3.5.1 httpuv_1.6.6
[99] rtracklayer_1.58.0 R6_2.5.1
[101] BiocIO_1.8.0 bookdown_0.29
[103] promises_1.2.0.1 gridExtra_2.3
[105] vipor_0.4.5 codetools_0.2-18
[107] assertthat_0.2.1 rjson_0.2.21
[109] withr_2.5.0 GenomicAlignments_1.34.0
[111] Rsamtools_2.14.0 GenomeInfoDbData_1.2.9
[113] parallel_4.2.1 hms_1.1.2
[115] grid_4.2.1 beachmat_2.14.0
[117] rmarkdown_2.17 DelayedMatrixStats_1.20.0
[119] Rtsne_0.16 shiny_1.7.3
[121] ggbeeswarm_0.6.0 restfulr_0.0.15
References
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.