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 571We 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
)
Figure 4.1: Distribution of each QC metric in the PBMC dataset, where each point is a cell and is colored by whether or not it was discarded by the outlier-based QC approach.
We also plot the mitochondrial proportion against the total count for each cell, as one does (Figure 4.2).
Figure 4.2: Percentage of UMIs mapped to mitochondrial genes against the totalcount for each cell.
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 <- applySCE(sce.pbmc, logNormCounts)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.36We 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")
Figure 4.3: Plot of the deconvolution size factors for the gene expression values (left) or the median-based size factors for the ADT expression values (right) compared to the library size-derived factors for the corresponding set of features. Each point represents a cell.
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 17library(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))
Figure 4.4: Heatmap of the pairwise cluster modularity scores in the PBMC dataset, computed based on the shared nearest neighbor graph derived from the ADT expression values.
Figure 4.5: Obligatory \(t\)-SNE plot of PBMC dataset based on its ADT expression values, where each point is a cell and is colored by the cluster of origin. Cluster labels are also overlaid at the median coordinates across all cells in the cluster.
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 1Session Info
R version 4.3.0 RC (2023-04-13 r84269)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.2 LTS
Matrix products: default
BLAS: /home/biocbuild/bbs-3.17-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] pheatmap_1.0.12 bluster_1.10.0
[3] scran_1.28.0 scater_1.28.0
[5] ggplot2_3.4.2 scuttle_1.10.0
[7] DropletUtils_1.20.0 SingleCellExperiment_1.22.0
[9] SummarizedExperiment_1.30.0 Biobase_2.60.0
[11] GenomicRanges_1.52.0 GenomeInfoDb_1.36.0
[13] IRanges_2.34.0 S4Vectors_0.38.0
[15] BiocGenerics_0.46.0 MatrixGenerics_1.12.0
[17] matrixStats_0.63.0 BiocFileCache_2.8.0
[19] dbplyr_2.3.2 BiocStyle_2.28.0
[21] rebook_1.10.0
loaded via a namespace (and not attached):
[1] DBI_1.1.3 bitops_1.0-7
[3] gridExtra_2.3 CodeDepends_0.6.5
[5] rlang_1.1.0 magrittr_2.0.3
[7] RcppAnnoy_0.0.20 compiler_4.3.0
[9] RSQLite_2.3.1 dir.expiry_1.8.0
[11] DelayedMatrixStats_1.22.0 vctrs_0.6.2
[13] pkgconfig_2.0.3 fastmap_1.1.1
[15] XVector_0.40.0 labeling_0.4.2
[17] utf8_1.2.3 rmarkdown_2.21
[19] graph_1.78.0 ggbeeswarm_0.7.1
[21] purrr_1.0.1 bit_4.0.5
[23] xfun_0.39 zlibbioc_1.46.0
[25] cachem_1.0.7 beachmat_2.16.0
[27] jsonlite_1.8.4 blob_1.2.4
[29] highr_0.10 rhdf5filters_1.12.0
[31] DelayedArray_0.26.0 Rhdf5lib_1.22.0
[33] BiocParallel_1.34.0 cluster_2.1.4
[35] irlba_2.3.5.1 parallel_4.3.0
[37] R6_2.5.1 RColorBrewer_1.1-3
[39] bslib_0.4.2 limma_3.56.0
[41] jquerylib_0.1.4 Rcpp_1.0.10
[43] bookdown_0.33 knitr_1.42
[45] R.utils_2.12.2 igraph_1.4.2
[47] Matrix_1.5-4 tidyselect_1.2.0
[49] viridis_0.6.2 yaml_2.3.7
[51] codetools_0.2-19 curl_5.0.0
[53] lattice_0.21-8 tibble_3.2.1
[55] withr_2.5.0 Rtsne_0.16
[57] evaluate_0.20 pillar_1.9.0
[59] BiocManager_1.30.20 filelock_1.0.2
[61] generics_0.1.3 RCurl_1.98-1.12
[63] sparseMatrixStats_1.12.0 munsell_0.5.0
[65] scales_1.2.1 glue_1.6.2
[67] metapod_1.8.0 tools_4.3.0
[69] BiocNeighbors_1.18.0 ScaledMatrix_1.8.0
[71] locfit_1.5-9.7 XML_3.99-0.14
[73] cowplot_1.1.1 rhdf5_2.44.0
[75] grid_4.3.0 edgeR_3.42.0
[77] colorspace_2.1-0 GenomeInfoDbData_1.2.10
[79] beeswarm_0.4.0 BiocSingular_1.16.0
[81] HDF5Array_1.28.0 vipor_0.4.5
[83] cli_3.6.1 rsvd_1.0.5
[85] fansi_1.0.4 viridisLite_0.4.1
[87] dplyr_1.1.2 uwot_0.1.14
[89] gtable_0.3.3 R.methodsS3_1.8.2
[91] sass_0.4.5 digest_0.6.31
[93] ggrepel_0.9.3 dqrng_0.3.0
[95] farver_2.1.1 memoise_2.0.1
[97] htmltools_0.5.5 R.oo_1.25.0
[99] lifecycle_1.0.3 httr_1.4.5
[101] statmod_1.5.0 bit64_4.0.5 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.