Chapter 4 Human PBMC with surface proteins (10X Genomics)
4.1 Introduction
Here, we describe a brief analysis of yet another peripheral blood mononuclear cell (PBMC) dataset from 10X Genomics (Zheng et al. 2017). Data are publicly available from the 10X Genomics website, from which we download the filtered gene/barcode count matrices for gene expression and cell surface proteins.
4.2 Data loading
library(BiocFileCache)
bfc <- BiocFileCache(ask=FALSE)
exprs.data <- bfcrpath(bfc, file.path(
"http://cf.10xgenomics.com/samples/cell-vdj/3.1.0",
"vdj_v1_hs_pbmc3",
"vdj_v1_hs_pbmc3_filtered_feature_bc_matrix.tar.gz"))
untar(exprs.data, exdir=tempdir())
library(DropletUtils)
sce.pbmc <- read10xCounts(file.path(tempdir(), "filtered_feature_bc_matrix"))
sce.pbmc <- splitAltExps(sce.pbmc, rowData(sce.pbmc)$Type)
4.3 Quality control
We discard cells with high mitochondrial proportions and few detectable ADT counts.
library(scater)
is.mito <- grep("^MT-", rowData(sce.pbmc)$Symbol)
stats <- perCellQCMetrics(sce.pbmc, subsets=list(Mito=is.mito))
high.mito <- isOutlier(stats$subsets_Mito_percent, type="higher")
low.adt <- stats$`altexps_Antibody Capture_detected` < nrow(altExp(sce.pbmc))/2
discard <- high.mito | low.adt
sce.pbmc <- sce.pbmc[,!discard]
We examine some of the statistics:
## Mode FALSE TRUE
## logical 6660 571
## Mode FALSE
## logical 7231
## Mode FALSE TRUE
## logical 6660 571
We examine the distribution of each QC metric (Figure 4.1).
colData(unfiltered) <- cbind(colData(unfiltered), stats)
unfiltered$discard <- discard
gridExtra::grid.arrange(
plotColData(unfiltered, y="sum", colour_by="discard") +
scale_y_log10() + ggtitle("Total count"),
plotColData(unfiltered, y="detected", colour_by="discard") +
scale_y_log10() + ggtitle("Detected features"),
plotColData(unfiltered, y="subsets_Mito_percent",
colour_by="discard") + ggtitle("Mito percent"),
plotColData(unfiltered, y="altexps_Antibody Capture_detected",
colour_by="discard") + ggtitle("ADT detected"),
ncol=2
)
We also plot the mitochondrial proportion against the total count for each cell, as one does (Figure 4.2).
4.4 Normalization
Computing size factors for the gene expression and ADT counts.
library(scran)
set.seed(1000)
clusters <- quickCluster(sce.pbmc)
sce.pbmc <- computeSumFactors(sce.pbmc, cluster=clusters)
altExp(sce.pbmc) <- computeMedianFactors(altExp(sce.pbmc))
sce.pbmc <- logNormCounts(sce.pbmc, use_altexps=TRUE)
We generate some summary statistics for both sets of size factors:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.074 0.719 0.908 1.000 1.133 8.858
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.10 0.70 0.83 1.00 1.03 227.36
We also look at the distribution of size factors compared to the library size for each set of features (Figure 4.3).
par(mfrow=c(1,2))
plot(librarySizeFactors(sce.pbmc), sizeFactors(sce.pbmc), pch=16,
xlab="Library size factors", ylab="Deconvolution factors",
main="Gene expression", log="xy")
plot(librarySizeFactors(altExp(sce.pbmc)), sizeFactors(altExp(sce.pbmc)), pch=16,
xlab="Library size factors", ylab="Median-based factors",
main="Antibody capture", log="xy")
4.5 Dimensionality reduction
We omit the PCA step for the ADT expression matrix, given that it is already so low-dimensional, and progress directly to \(t\)-SNE and UMAP visualizations.
4.6 Clustering
We perform graph-based clustering on the ADT data and use the assignments as the column labels of the alternative Experiment.
g.adt <- buildSNNGraph(altExp(sce.pbmc), k=10, d=NA)
clust.adt <- igraph::cluster_walktrap(g.adt)$membership
colLabels(altExp(sce.pbmc)) <- factor(clust.adt)
We examine some basic statistics about the size of each cluster, their separation (Figure 4.4) and their distribution in our \(t\)-SNE plot (Figure 4.5).
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 160 507 662 39 691 1415 32 650 76 1037 121 47 68 25 15 562
## 17 18 19 20 21 22 23 24
## 139 32 44 120 84 65 52 17
library(bluster)
mod <- pairwiseModularity(g.adt, clust.adt, as.ratio=TRUE)
library(pheatmap)
pheatmap::pheatmap(log10(mod + 10), cluster_row=FALSE, cluster_col=FALSE,
color=colorRampPalette(c("white", "blue"))(101))
We perform some additional subclustering using the expression data to mimic an in silico FACS experiment.
set.seed(1010010)
subclusters <- quickSubCluster(sce.pbmc, clust.adt,
prepFUN=function(x) {
dec <- modelGeneVarByPoisson(x)
top <- getTopHVGs(dec, prop=0.1)
denoisePCA(x, dec, subset.row=top)
},
clusterFUN=function(x) {
g.gene <- buildSNNGraph(x, k=10, use.dimred = 'PCA')
igraph::cluster_walktrap(g.gene)$membership
}
)
We counting the number of gene expression-derived subclusters in each ADT-derived parent cluster.
data.frame(
Cluster=names(subclusters),
Ncells=vapply(subclusters, ncol, 0L),
Nsub=vapply(subclusters, function(x) length(unique(x$subcluster)), 0L)
)
## Cluster Ncells Nsub
## 1 1 160 3
## 2 2 507 4
## 3 3 662 5
## 4 4 39 1
## 5 5 691 5
## 6 6 1415 7
## 7 7 32 1
## 8 8 650 7
## 9 9 76 2
## 10 10 1037 8
## 11 11 121 2
## 12 12 47 1
## 13 13 68 2
## 14 14 25 1
## 15 15 15 1
## 16 16 562 9
## 17 17 139 3
## 18 18 32 1
## 19 19 44 1
## 20 20 120 4
## 21 21 84 3
## 22 22 65 2
## 23 23 52 3
## 24 24 17 1
Session Info
R version 4.1.1 (2021-08-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.3 LTS
Matrix products: default
BLAS: /home/biocbuild/bbs-3.14-bioc/R/lib/libRblas.so
LAPACK: /home/biocbuild/bbs-3.14-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] pheatmap_1.0.12 bluster_1.4.0
[3] scran_1.22.0 scater_1.22.0
[5] ggplot2_3.3.5 scuttle_1.4.0
[7] DropletUtils_1.14.0 SingleCellExperiment_1.16.0
[9] SummarizedExperiment_1.24.0 Biobase_2.54.0
[11] GenomicRanges_1.46.0 GenomeInfoDb_1.30.0
[13] IRanges_2.28.0 S4Vectors_0.32.0
[15] BiocGenerics_0.40.0 MatrixGenerics_1.6.0
[17] matrixStats_0.61.0 BiocFileCache_2.2.0
[19] dbplyr_2.1.1 BiocStyle_2.22.0
[21] rebook_1.4.0
loaded via a namespace (and not attached):
[1] Rtsne_0.15 ggbeeswarm_0.6.0
[3] colorspace_2.0-2 ellipsis_0.3.2
[5] XVector_0.34.0 BiocNeighbors_1.12.0
[7] farver_2.1.0 ggrepel_0.9.1
[9] bit64_4.0.5 RSpectra_0.16-0
[11] fansi_0.5.0 codetools_0.2-18
[13] R.methodsS3_1.8.1 sparseMatrixStats_1.6.0
[15] cachem_1.0.6 knitr_1.36
[17] jsonlite_1.7.2 cluster_2.1.2
[19] R.oo_1.24.0 uwot_0.1.10
[21] graph_1.72.0 HDF5Array_1.22.0
[23] BiocManager_1.30.16 compiler_4.1.1
[25] httr_1.4.2 dqrng_0.3.0
[27] assertthat_0.2.1 Matrix_1.3-4
[29] fastmap_1.1.0 limma_3.50.0
[31] BiocSingular_1.10.0 htmltools_0.5.2
[33] tools_4.1.1 igraph_1.2.7
[35] rsvd_1.0.5 gtable_0.3.0
[37] glue_1.4.2 GenomeInfoDbData_1.2.7
[39] dplyr_1.0.7 rappdirs_0.3.3
[41] Rcpp_1.0.7 jquerylib_0.1.4
[43] vctrs_0.3.8 rhdf5filters_1.6.0
[45] DelayedMatrixStats_1.16.0 xfun_0.27
[47] stringr_1.4.0 beachmat_2.10.0
[49] lifecycle_1.0.1 irlba_2.3.3
[51] statmod_1.4.36 XML_3.99-0.8
[53] edgeR_3.36.0 zlibbioc_1.40.0
[55] scales_1.1.1 parallel_4.1.1
[57] rhdf5_2.38.0 RColorBrewer_1.1-2
[59] yaml_2.2.1 curl_4.3.2
[61] memoise_2.0.0 gridExtra_2.3
[63] sass_0.4.0 stringi_1.7.5
[65] RSQLite_2.2.8 highr_0.9
[67] ScaledMatrix_1.2.0 filelock_1.0.2
[69] BiocParallel_1.28.0 rlang_0.4.12
[71] pkgconfig_2.0.3 bitops_1.0-7
[73] evaluate_0.14 lattice_0.20-45
[75] purrr_0.3.4 Rhdf5lib_1.16.0
[77] CodeDepends_0.6.5 labeling_0.4.2
[79] cowplot_1.1.1 bit_4.0.4
[81] tidyselect_1.1.1 RcppAnnoy_0.0.19
[83] magrittr_2.0.1 bookdown_0.24
[85] R6_2.5.1 generics_0.1.1
[87] metapod_1.2.0 DelayedArray_0.20.0
[89] DBI_1.1.1 pillar_1.6.4
[91] withr_2.4.2 RCurl_1.98-1.5
[93] tibble_3.1.5 dir.expiry_1.2.0
[95] crayon_1.4.1 utf8_1.2.2
[97] rmarkdown_2.11 viridis_0.6.2
[99] locfit_1.5-9.4 grid_4.1.1
[101] blob_1.2.2 digest_0.6.28
[103] R.utils_2.11.0 munsell_0.5.0
[105] beeswarm_0.4.0 viridisLite_0.4.0
[107] vipor_0.4.5 bslib_0.3.1
References
Zheng, G. X., J. M. Terry, P. Belgrader, P. Ryvkin, Z. W. Bent, R. Wilson, S. B. Ziraldo, et al. 2017. “Massively parallel digital transcriptional profiling of single cells.” Nat Commun 8 (January): 14049.