PhEMD是一个包比较多个异质单细胞样品彼此。目前,它通过首先定义细胞亚型并使用Monocle 2将它们相互关联来实现这一点。然后计算每对样本之间的地球移动距离(EMD),结合有关内在细胞亚型不相似性的信息(即细胞亚型之间的流形距离)和样本之间相对于每个细胞亚型的相对丰度的差异。PhEMD使用这些两两距离来构建单细胞样本的网络,这些单细胞样本可能在2D或3D中具有视觉关联,并使用扩散图和分区来识别相似样本的组。
PhEMD需要R版本>= 3.4.0(建议3.5.0+),Bioconductor版本>= 3.5(建议3.7+),Monocle 2版本>= 2.4.0(建议2.8.0)
###从Bioconductor安装
BiocManager::安装(“phemd”)
###从Github安装(但最好直接从Bioconductor安装)
库(devtools) install_github(“wschen / phemd”)
###安装后加载库
库(phemd)库(单片眼镜)
PhEMD期望单细胞数据表示为样本的R列表。每个样本(即列表元素)应该是一个维度的矩阵num_cellsxnum_markers,其中标记可以代表基因或细胞计数蛋白标记。对于这篇小短文,我们将展示我们对黑色素瘤数据集的分析管道,该数据集由肿瘤浸润性免疫细胞scRNA-seq数据(选定的基因)组成,这些数据在tpm归一化后进行对数转换(由Tirosh等人首次发表,2016年)。
我们首先创建一个PhEMD数据对象,指定多样本表达式数据(R列表)、标记名(即表达式数据列表中数据矩阵的列名)和样本名(按照它们在表达式数据列表中出现的顺序)。
load('melanomaData.RData') myobj <- createDataObj(all_expn_data, all_genes, as.character(snames))
我们可以选择删除PhEMD数据对象中单元格数少于min_sz的样本,如下所示:
myobj <- removetinyssamples (myobj, min_sz = 20)
## [1] "Mel59因为只包含12个单元格而被删除"
请注意,不符合最小细胞产量标准的样本将从rawExpn(myobj)和sampleNames(myobj)中的样本名称列表中删除。
接下来,将来自所有样本的数据聚合到存储在PhEMD数据对象(在槽' data_aggregate '中)中的矩阵中。然后,这些聚合的数据将用于初始单元格子类型定义和嵌入。如果所有样本中的细胞总数大于max_cells,则每个样本中的细胞数量相等,并存储在pooledCells(myobj)中。
myobj <- aggregatessamples (myobj, max_cells=12000)
现在我们已经汇总了来自所有样本的单细胞数据,我们准备执行细胞亚型定义和降维,以在视觉上和空间上关联细胞和细胞亚型。为此,我们使用Monocle 2。在我们开始之前,我们首先通过选择44个重要基因进行特征选择。关于如何选择重要基因的建议可以在这里找到: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)
现在我们准备生成一个Monocle 2嵌入。我们的embedCells ()的包装器函数reduceDimension ()函数在Monocle 2。对于我们的示例数据集,我们指定表达式分布模型为“gaussianff”这是Monocle 2教程中关于日志转换scRNA-seq TPM值的建议(http://cole-trapnell-lab.github.io/monocle-release/docs/#choosing-a-distribution-for-your-data-required).“negbinomial_sz”对于大多数非规范化scRNA-seq数据(原始读取计数),推荐的数据类型是什么“gaussianff”推荐用于对数转换数据或arcsin转换的细胞计数数据。请参阅上面的链接了解更多细节。
附加的参数可以传递给Monocle 2reduceDimension ()中的可选命名参数embed_cells ().我们发现Monocle 2对一系列参数具有鲁棒性。Sigma可以被认为是一个“噪声”参数,我们的经验发现,在[0.01,0.1]范围内的Sigma通常适用于对数变换的scRNA-seq数据或归一化的CyTOF数据。西格玛值越大,通常导致集群总数越少。参见Monocle 2出版物(Qiu et al., 2017)了解参数选择的更多细节。
#生成二维细胞嵌入和细胞子类型分配myobj <- embedCells(myobj, data_model = 'gaussianff', pseudo_expr=0, sigma=0.02) #生成细胞的伪时间排序myobj <- orderCellsMonocle(myobj)
上面代码的结果是存储在中的Monocle 2对象monocleInfo (myobj).该对象包含聚合数据矩阵中每个单元格的子类型和伪时间赋值pooledCells (myobj)).这些细胞的二维嵌入也已生成。我们可以通过这样将它们写入文件来可视化嵌入:
cmap <- plotEmbeddings(myobj, cell_model='monocle2')
为了可视化细胞亚型的表达谱,我们可以绘制热图并保存到文件,如下所示:
plotHeatmaps(myobj, cell_model='monocle2', selected_genes= ' heatmap_genes ')
现在我们已经在所有单细胞样本中确定了一组全面的细胞亚型,并通过聚集来自所有样本的细胞将它们在低维嵌入中联系起来,我们想要执行反卷积以确定每个样本上每种细胞亚型的丰度。为此,我们调用这个函数:
myobj <- clusterIndividualSamples(myobj)
此过程的结果存储在celltypeFreqs (myobj).行我列j表示样本中所有细胞的比例我分配给细胞子类型j.
为了比较单细胞样本,我们使用了地球移动距离,这是一个考虑到匹配细胞亚型的相对频率差异(例如,每个样本中所有细胞中CD8+ t细胞的百分比)和细胞亚型本身的差异性(例如,CD8+和CD4+ t细胞之间的内在差异性)的度量。为了计算细胞亚型之间的内在不相似性,我们调用以下函数:
#确定(dis)不同细胞亚型的相似性myobj <- generateGDM(myobj)
generateGDM ()中存储单元格亚型之间的两两不相似性(即“地面距离”或“树距离”)GDM (myobj).
现在我们准备使用EMD来比较单细胞样品。为此,只需调用该函数compareSamples ():
使用EMD my_distmat <- compareSamples(myobj)执行样本间比较
compareSamples ()返回表示单细胞样本之间成对EMD的距离矩阵;my_distmat (i, j)表示样本间的不相似性我和j(即以行表示的样本我和j在celltypeFreqs (myobj))。我们可以使用这个距离矩阵来识别类似的样本组,如下所示:
##识别相似的抑制因子组group_assignments <- groupSamples(my_distmat, distfun = ' hcluster ', ncluster=5)
我们还可以使用基于phemd的距离矩阵来生成单细胞样本的嵌入,通过分组分配着色。
dmap_obj <- plotGroupedSamplesDmap(my_distmat, group_assignments, pt_sz = 1.5)
##警告:as(, "dsTMatrix")自Matrix 1.5-0起已弃用;做……,“TsparseMatrix”)
要在每个样本的基础上检索单元格子类型分布,请使用以下函数。随后可以根据需要绘制给定样本的直方图。
绘制每个样本的细胞亚型分布(直方图)。cellfreqs <- getSampleHistsByCluster(myobj, group_assignments, cell_model='monocle2')
要绘制总结相似样本组的细胞亚型直方图(分配给特定组的所有样本中每个细胞亚型的双象平均值),请使用以下绘图函数:
#绘制每组样本的代表性细胞子类型分布plotSummaryHistograms(myobj, group_assignments, cell_model='monocle2', cmap, ncol。Plot = 5, ax.lab.sz=1.3, title.sz=2)
要将每个样品的细胞产率绘制为柱状图,请使用以下函数:
Plot cell yield (myobj, group_assignments, font_sz = 0.7, w=8, h=5)
sessionInfo ()
## R版本4.2.3(2023-03-15)##平台:x86_64-pc-linux-gnu(64位)##运行在:Ubuntu 20.04.6 LTS ## ##矩阵产品:默认## BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas。/home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack。所以## ## 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 ## ##附加的基本包:## [1]splines stats4 stats graphics grDevices utils数据集## [8]methods base ## ##其他附加包:## [1]phemd_1.14.1 SeuratObject_4.1.3 Seurat_4.3.0 ## [4] monocle_2.26.0 DDRTree_0.1.5 irlba_2.3.5.1 ## [7] VGAM_1.1-8 ggplot2_1 .4.1 Biobase_2.58.0 ## [10] BiocGenerics_0.44.0 Matrix_1.5-3 BiocStyle_2.26.0 ## ##通过命名空间加载(未附加):# # # # [1] utf8_1.2.3 spatstat.explore_3.1-0 [3] reticulate_1.28 tidyselect_1.2.0 # # [5] htmlwidgets_1.6.2 grid_4.2.3 # # [7] combinat_0.0-8 ranger_0.14.1 # # [9] docopt_0.7.1 Rtsne_0.16 # # [11] munsell_0.5.0 destiny_3.12.0 # # [13] codetools_0.2-19 ica_1.0-3 # # [15] future_1.32.0 miniUI_0.1.1.1 # # [17] withr_2.5.0 spatstat.random_3.1-4 # # [19] colorspace_2.1-0 fastICA_1.2-3 # # [21] progressr_0.13.0 highr_0.10 # # [23] knitr_1.42 SingleCellExperiment_1.20.1 # # [25] ROCR_1.0-11 robustbase_0.95-0 # #[27] vcd_1.4-11 tensor_1.5 # # [29] VIM_6.2.2 TTR_0.24.3 # # [31] listenv_0.9.0 labeling_0.4.2 # # [33] MatrixGenerics_1.10.0 slam_0.1-50 # # [35] GenomeInfoDbData_1.2.9 polyclip_1.10-4 # # [37] farver_2.1.1 pheatmap_1.0.12 # # [39] parallelly_1.35.0 vctrs_0.6.1 # # [41] generics_0.1.3 xfun_0.38 # # [43] ggthemes_4.2.4 phateR_1.0.7 # # [45] R6_2.5.1 GenomeInfoDb_1.34.9 # # [47] RcppEigen_0.3.3.9.3 bitops_1.0-7 # # [49] spatstat.utils_3.0-2 cachem_1.0.7 # # [51] DelayedArray_0.24.0 promises_1.2.0.1 # # [53]scales_1.2.1 nnet_1 . 7.3-18 ## [55] gtable_1 .3.3 globals_0.16.2 ## [57] goftest_1. 3- 3 rlang_1.1.0 ## [57] scatterplot3d_0.3-43 lazyeval_1 .2.2 ## [61] hexbin_1.28.3 spatstat.geom_3.1-0 ## [63] BiocManager_1.30.20 yaml_2.3.7 ## [65] shape2_1.4.4 abind_1.4-5 ## [67] httpuv_1.6.9 tools_4.2.3 ## [69] bookdown_1 .33 ellipsis_0.3.2 ## [71] jquerylib_0.1.4 rcolorbrewer_1 .3- 3 ## [73] proxy_0.4-27 griridges_0.5.4 ## [75] Rcpp_1.0.10 plyr_1.8.8 ## [77] zlibbioc_1.44.0 purrr_1.0.1 ## [79] RCurl_1.98-1.12rpart_4.1.19 ## [81] deldir_1.0-6 pbapply_1.7-0 ## [83] viridis_0.6.2 cowplot_1.1.1 ## [85] S4Vectors_0.36.2 zoo_1.8-11 ## [87] SummarizedExperiment_1.28.0 ggrepel_0.9.3 ## [89] cluster_0.1.3 magrittr_1 .0.3 ## [91] magick_2.7.4 data.table_1.14.8 ## [93] RSpectra_0.16-1 scattermore_0.8 ## [95] lmtest_0.9-40 RANN_2.6.1 ## [97] pcaMethods_1.90.0 fitdistrplus_1.1-8 ## [99] matrixStats_0.63.0 patchwork_1.1.2 ## [101] mime_0.12 evaluate_0.20 ## [103] xtable_1.8-4 smoother_1.1 ## [105] sparsesvd_0.2-2IRanges_2.32.0 ## [107] gridextra2.3 HSMMSingleCell_1.18.0 ## [109] compiler_4.2.3 transport_0.13-0 ## [111] tibble_1 .2.1 KernSmooth_2.23-20 ## [113] htmltools_0.5.5 later_1.3.0 ## [115] tidyr_1.3.0 MASS_7.3-58.3 ## [117] boot_1.3-28.1 leidenbase_1 .1.14 ## [119] car_3.1-1 cli_3.6.1 ## [121] parallel_4.2.3 igraph_1.4.1 ## [123] genomicranges_1 . 1.50.2 pkgconfig_2.0.3 ## [125] laeken_0.5.2 sp_1. 0 ## [127] plotly_4.10.1 spatstat.sparse_3.0-1 ## [129] bslib_1 .4.2 xvector_1 .38.0 ## [131] stringr_1.5.0 digest_0.6.31 ## [133] sctransform_0.3.5 RcppAnnoy_0.0.20 ## [135] pracma_2.4.2 spatstat.data_3.0-1 ## [137] rmarkdown_2.21 leiden_0.4.3 ## [139] uwot_0.1.14 curl_5.0.0 ## [141] ggplot.multistats_1.0.0 shiny_1.7.4 ## [143] lifecycle_1.0.3 nlme_3.1-162 ## [145] jsonlite_1.8.4 carData_3.0-5 ## [147] viridisLite_0.4.1 limma_3.54.2 ## [149] fansi_1.0.4 pillar_1.9.0 ## [151] lattice_0.20-45 fastmap_1.1.1 ## [153] httr_1.4.5 DEoptimR_1.0-11 ## [155] survival_3.5-5 xts_0.13.0 ## [157] glue_1.6.2 qlcMatrix_0.9.7 ## [159] png_0.1-8 class_7.3-21 ## [161] stringi_1.7.12 sass_0.4.5 ## [163] RcppHNSW_0.4.1 memoise_2.0.1 ## [165] dplyr_1.1.1 maptree_1.4-8 ## [167] e1071_1.7-13 future.apply_1.10.0