Tutorial on using PhEMD to analyze multi-sample single-cell experiments
October 20, 2018
Contents
- 0.1 Overview
- 0.2 1. Installation
- 0.3 2. Preparing data for cell state definition and embedding
- 0.4 3. Generate Monocle 2 cell embedding with cell state definitions
- 0.5 4. Deconvolute single-cell samples and compare using Earth Mover’s Distance
- 0.6 5. Visualize single-cell samples based on PhEMD-based similarity
0.1 Overview
PhEMD is a package for comparing multiple heterogeneous single-cell samples to one another. It currently does so by first defining cell subtypes and relating them to one another using Monocle 2. It then computes Earth Mover’s Distance (EMD) between each pair of samples, incorporating information about intrinsic cell subtype dissimilarity (i.e. manifold distance between cell subtypes) and differences between samples with respect to relative abundance of each cell subtype. PhEMD uses these pairwise distances to construct a network of single-cell samples that may be visually related in 2D or 3D using a diffusion map and partitioned to identify groups of similar samples.
0.2 1. Installation
PhEMD requires R version >= 3.4.0 (recommended 3.5.0+), Bioconductor version >= 3.5 (recommended 3.7+), and Monocle 2 version >= 2.4.0 (recommended 2.8.0)
###Install from Bioconductor
BiocManager::install("phemd")###Install from Github (though direct installation from Bioconductor is preferred)
library(devtools)
install_github("wschen/phemd")###Load library after installation
library('phemd')
library('monocle')0.3 2. Preparing data for cell state definition and embedding
PhEMD expects single-cell data to be represented as an R list of samples. Each sample (i.e. list element) is expected to be a matrix of dimension num_cells x num_markers, where markers may represent genes or cytometry protein markers. For this vignette, we will be demonstrating our analysis pipeline on a melanoma dataset consisting of tumor-infiltrating immune cell scRNA-seq data (selected genes) that were log-transformed following TPM-normalization (first published by Tirosh et al., 2016).
We first start by creating a PhEMD data object, specifying the multi-sample expression data (R list), marker names (i.e. column names of the data matrices in the list of expression data), and sample names (in the order they appear in the list of expression data).
load('melanomaData.RData')
myobj <- createDataObj(all_expn_data, all_genes, as.character(snames))We can optionally remove samples in the PhEMD data object that have fewer than min_sz number of cells as follows:
myobj <- removeTinySamples(myobj, min_sz = 20)## [1] "Mel78 removed because only contains 3 cells"
## [1] "Mel59 removed because only contains 12 cells"Note that samples that don’t meet the meet the minimum cell yield criteria are removed from rawExpn(myobj) and from the list of sample names in sampleNames(myobj).
Next, aggregate data from all samples into a matrix that is stored in the PhEMD data object (in slot ‘data_aggregate’). This aggregated data will then be used for initial cell subtype definition and embedding. If there are more cells collectively in all samples than max_cells, an equal number of cells from each sample will be subsampled and stored in pooledCells(myobj).
myobj <- aggregateSamples(myobj, max_cells=12000)0.4 3. Generate Monocle 2 cell embedding with cell state definitions
Now that we have aggregated single-cell data from all samples, we are ready to perform cell subtype definition and dimensionality reduction to visually and spatially relate cells and cell subtypes. For this, we use Monocle 2. Before we begin, we first perform feature selection by selecting 44 important genes. Suggestions on how to choose important genes can be found here: http://cole-trapnell-lab.github.io/monocle-release/docs/#trajectory-step-1-choose-genes-that-define-a-cell-s-progress
myobj <- selectFeatures(myobj, selected_genes)We are now ready to generate a Monocle 2 embedding. Our embedCells() function is a wrapper function for the reduceDimension() function in Monocle 2. For our example dataset, we specify the expression distribution model as ‘gaussianff’ as is recommended in the Monocle 2 tutorial for log-transformed scRNA-seq TPM values (http://cole-trapnell-lab.github.io/monocle-release/docs/#choosing-a-distribution-for-your-data-required). ‘negbinomial_sz’ is the recommended data type for most unnormalized scRNA-seq data (raw read counts) and ‘gaussianff’ is recommended for log-transformed data or arcsin-transformed mass cytometry data. See above link for more details.
Additional parameters may be passed to Monocle 2 reduceDimension() as optional named parameters in embed_cells(). We found that Monocle 2 is robust to a range of parameters. Sigma can be thought of as a “noise” parameter and we empirically found that sigma in the range of [0.01, 0.1] often works well for log-transformed scRNA-seq data or normalized CyTOF data. Greater values of sigma generally result in fewer total number of clusters. See Monocle 2 publication (Qiu et al., 2017) for additional details on parameter selection.
# generate 2D cell embedding and cell subtype assignments
myobj <- embedCells(myobj, data_model = 'gaussianff', pseudo_expr=0, sigma=0.02)
# generate pseudotime ordering of cells
myobj <- orderCellsMonocle(myobj)The result of the code above is a Monocle 2 object stored in monocleInfo(myobj). This object contains cell subtype and pseudotime assignments for each cell in the aggregated data matrix (stored in pooledCells(myobj)). A 2D embedding of these cells has also been generated. We can visualize the embedding by writing them to file in this way:
cmap <- plotEmbeddings(myobj, cell_model='monocle2')
To visualize the expression profiles of the cell subtypes, we can plot a heatmap and save to file as such:
plotHeatmaps(myobj, cell_model='monocle2', selected_genes=heatmap_genes)
0.5 4. Deconvolute single-cell samples and compare using Earth Mover’s Distance
Now that we have identified a comprehensive set of cell subtypes across all single-cell samples and related them in a low-dimensional embedding by aggregating cells from all samples, we want to perform deconvolution to determine the abundance of each cell subtype on a per sample basis. To do so, we call this function:
# Determine cell subtype breakdown of each sample
myobj <- clusterIndividualSamples(myobj)The results of this process are stored in celltypeFreqs(myobj). Row i column j represents the fraction of all cells in sample i assigned to cell subtype j.
To compare single-cell samples, we use Earth Mover’s Distance, which is a metric that takes into account both the difference in relative frequencies of matching cell subtypes (e.g. % of all cells in each sample that are CD8+ T-cells) and the dissimilarity of the cell subtypes themselves (e.g. intrinsic dissimilarity between CD8+ and CD4+ T-cells). To compute the intrinsic dissimilarity between cell subtypes, we call the following function:
# Determine (dis)similarity of different cell subtypes
myobj <- generateGDM(myobj)generateGDM() stores the pairwise dissimilarity (i.e. “ground-distance” or “tree-distance”) between cell subtypes in GDM(myobj).
We are now ready to compare single-cell samples using EMD. To do so, we simply call the function compareSamples():
# Perform inter-sample comparisons using EMD
my_distmat <- compareSamples(myobj)compareSamples() returns a distance matrix representing the pairwise EMD between single-cell samples; my_distmat[i,j] represents the dissimilarity between samples i and j (i.e. samples represented by rows i and j in celltypeFreqs(myobj)). We can use this distance matrix to identify similar groups of samples as such:
## Identify similar groups of inhibitors
group_assignments <- groupSamples(my_distmat, distfun = 'hclust', ncluster=5)0.6 5. Visualize single-cell samples based on PhEMD-based similarity
We can also use the PhEMD-based distance matrix to generate an embedding of single-cell samples, colored by group assignments.
dmap_obj <- plotGroupedSamplesDmap(my_distmat, group_assignments, pt_sz = 1.5)
To retrieve the cell subtype distribution on a per-sample basis, use the following function. Histograms can be subsequently plotted for a given sample as desired.
# Plot cell subtype distribution (histogram) for each sample
sample.cellfreqs <- getSampleHistsByCluster(myobj, group_assignments, cell_model='monocle2')To plot cell subtype histograms summarizing groups of similar samples (bin-wise mean of each cell subtype across all samples assigned to a particular group), use the following plotting function:
# Plot representative cell subtype distribution for each group of samples
plotSummaryHistograms(myobj, group_assignments, cell_model='monocle2', cmap,
ncol.plot = 5, ax.lab.sz=1.3, title.sz=2)
To plot cell yield of each samples as a barplot, use the following function:
# Plot cell yield of each experimental condition
plotCellYield(myobj, group_assignments, font_sz = 0.7, w=8, h=5)
sessionInfo()## R version 4.1.3 (2022-03-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 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] splines stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] phemd_1.10.0 SeuratObject_4.0.4 Seurat_4.1.0
## [4] monocle_2.22.0 DDRTree_0.1.5 irlba_2.3.5
## [7] VGAM_1.1-6 ggplot2_3.3.5 Biobase_2.54.0
## [10] BiocGenerics_0.40.0 Matrix_1.4-1 BiocStyle_2.22.0
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 reticulate_1.24
## [3] tidyselect_1.1.2 htmlwidgets_1.5.4
## [5] ranger_0.13.1 grid_4.1.3
## [7] combinat_0.0-8 docopt_0.7.1
## [9] Rtsne_0.15 munsell_0.5.0
## [11] destiny_3.8.1 codetools_0.2-18
## [13] ica_1.0-2 future_1.24.0
## [15] miniUI_0.1.1.1 withr_2.5.0
## [17] spatstat.random_2.2-0 colorspace_2.0-3
## [19] fastICA_1.2-3 highr_0.9
## [21] knitr_1.38 SingleCellExperiment_1.16.0
## [23] ROCR_1.0-11 robustbase_0.95-0
## [25] vcd_1.4-9 tensor_1.5
## [27] VIM_6.1.1 TTR_0.24.3
## [29] listenv_0.8.0 labeling_0.4.2
## [31] MatrixGenerics_1.6.0 slam_0.1-50
## [33] GenomeInfoDbData_1.2.7 polyclip_1.10-0
## [35] farver_2.1.0 pheatmap_1.0.12
## [37] parallelly_1.30.0 vctrs_0.4.0
## [39] generics_0.1.2 xfun_0.30
## [41] ggthemes_4.2.4 phateR_1.0.7
## [43] R6_2.5.1 GenomeInfoDb_1.30.1
## [45] RcppEigen_0.3.3.9.1 cachem_1.0.6
## [47] bitops_1.0-7 spatstat.utils_2.3-0
## [49] DelayedArray_0.20.0 assertthat_0.2.1
## [51] promises_1.2.0.1 scales_1.1.1
## [53] nnet_7.3-17 gtable_0.3.0
## [55] globals_0.14.0 goftest_1.2-3
## [57] rlang_1.0.2 scatterplot3d_0.3-41
## [59] lazyeval_0.2.2 hexbin_1.28.2
## [61] spatstat.geom_2.4-0 BiocManager_1.30.16
## [63] yaml_2.3.5 reshape2_1.4.4
## [65] abind_1.4-5 httpuv_1.6.5
## [67] tools_4.1.3 bookdown_0.25
## [69] ellipsis_0.3.2 spatstat.core_2.4-2
## [71] jquerylib_0.1.4 RColorBrewer_1.1-3
## [73] proxy_0.4-26 ggridges_0.5.3
## [75] Rcpp_1.0.8.3 plyr_1.8.7
## [77] zlibbioc_1.40.0 purrr_0.3.4
## [79] RCurl_1.98-1.6 densityClust_0.3.2
## [81] rpart_4.1.16 deldir_1.0-6
## [83] pbapply_1.5-0 viridis_0.6.2
## [85] cowplot_1.1.1 S4Vectors_0.32.4
## [87] zoo_1.8-9 SummarizedExperiment_1.24.0
## [89] ggrepel_0.9.1 cluster_2.1.3
## [91] magrittr_2.0.3 magick_2.7.3
## [93] data.table_1.14.2 RSpectra_0.16-0
## [95] scattermore_0.8 lmtest_0.9-40
## [97] RANN_2.6.1 pcaMethods_1.86.0
## [99] fitdistrplus_1.1-8 matrixStats_0.61.0
## [101] patchwork_1.1.1 mime_0.12
## [103] evaluate_0.15 xtable_1.8-4
## [105] smoother_1.1 sparsesvd_0.2
## [107] IRanges_2.28.0 gridExtra_2.3
## [109] HSMMSingleCell_1.14.0 compiler_4.1.3
## [111] transport_0.12-2 tibble_3.1.6
## [113] KernSmooth_2.23-20 crayon_1.5.1
## [115] htmltools_0.5.2 mgcv_1.8-40
## [117] later_1.3.0 tidyr_1.2.0
## [119] DBI_1.1.2 MASS_7.3-56
## [121] boot_1.3-28 car_3.0-12
## [123] cli_3.2.0 parallel_4.1.3
## [125] igraph_1.3.0 GenomicRanges_1.46.1
## [127] pkgconfig_2.0.3 laeken_0.5.2
## [129] sp_1.4-6 plotly_4.10.0
## [131] spatstat.sparse_2.1-0 bslib_0.3.1
## [133] XVector_0.34.0 stringr_1.4.0
## [135] digest_0.6.29 sctransform_0.3.3
## [137] RcppAnnoy_0.0.19 pracma_2.3.8
## [139] spatstat.data_2.1-4 rmarkdown_2.13
## [141] leiden_0.3.9 uwot_0.1.11
## [143] curl_4.3.2 ggplot.multistats_1.0.0
## [145] shiny_1.7.1 lifecycle_1.0.1
## [147] nlme_3.1-157 jsonlite_1.8.0
## [149] carData_3.0-5 viridisLite_0.4.0
## [151] limma_3.50.1 fansi_1.0.3
## [153] pillar_1.7.0 lattice_0.20-45
## [155] fastmap_1.1.0 httr_1.4.2
## [157] DEoptimR_1.0-11 survival_3.3-1
## [159] xts_0.12.1 glue_1.6.2
## [161] qlcMatrix_0.9.7 FNN_1.1.3
## [163] png_0.1-7 class_7.3-20
## [165] stringi_1.7.6 sass_0.4.1
## [167] RcppHNSW_0.3.0 memoise_2.0.1
## [169] dplyr_1.0.8 maptree_1.4-7
## [171] e1071_1.7-9 future.apply_1.8.1