本文首发于“生信补给站”公众号 https://mp.weixin.qq.com/s/u50VeedGKnsCGIB-XsM4cw
拟时（pseudotime）分析，又称细胞轨迹（cell trajectory）分析，根据不同细胞亚群基因表达量随时间的变化情况通过拟时分析可以来推断发育过程细胞的分化轨迹或细胞亚型的演化过程。并非一定要不同时间段做实验的结果，因为细胞本身存在拟时序变化，细胞是有变化的，可以做拟时序分析。
本期分享2020年的Nature Communications文章 Single-cell analysis reveals new evolutionary complexity in uveal melanoma中的拟时序分析结果图，常规可视化的基础上添加了树形图，更直观。
一 加载数据 R包
本次介绍monocle2进行拟时序分析
#BiocManager::install("monocle") library(monocle) #载入注释后的数据 load('pbmc_tutorial_singleR.RData') pbmc #An object of class Seurat #13714 features across 2638 samples within 1 assay #Active assay: RNA (13714 features, 2000 variable features) # 3 dimensional reductions calculated: pca, umap, tsne
二 Monocle2分析
2.1 通过seurat结果导入数据
#Extract data, phenotype data, and feature data from the SeuratObject data <- as(as.matrix(pbmc@assays$RNA@counts), 'sparseMatrix') pd <- new('AnnotatedDataFrame', data = pbmc@meta.data) fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data)) fd <- new('AnnotatedDataFrame', data = fData) #构建S4对象，CellDataSet HSMM <- newCellDataSet(data, phenoData = pd, featureData = fd, lowerDetectionLimit = 0.5, expressionFamily = negbinomial.size())
注意expressionFamily的选择：
单细胞稀疏矩阵的话用negbinomial.size()，如果是UMI的话不要标准化；
FPKM值用tobit()；
logFPKM值用gaussianff()
2.2 估计size factors 和 dispersions （需要）
## HSMM <- estimateSizeFactors(HSMM) HSMM <- estimateDispersions(HSMM)
2.3 过滤低质量细胞
HSMM <- detectGenes(HSMM, min_expr = 3 ) print(head(fData(HSMM))) expressed_genes <- row.names(subset(fData(HSMM), num_cells_expressed >= 10)) head(pData(HSMM))
2.4 选择基因
diff_test_res <- differentialGeneTest(HSMM[expressed_genes,], fullModelFormulaStr = "~ seurat_clusters") ordering_genes <- row.names (subset(diff_test_res, qval < 0.1)) ## 不要也写0.1 ，而是要写0.01。 HSMM <- setOrderingFilter(HSMM, ordering_genes) plot_ordering_genes(HSMM)
这里选择基因的方式有很多，说明文档中建议以下4种选择基因的方式
（1）选择发育差异表达基因
（2）选择clusters差异表达基因
（3）选择离散程度高的基因
（4）自定义发育marker基因
2.5 降维 & 排序
HSMM <- reduceDimension(HSMM, max_components = 3, num_dim = 20, method = 'DDRTree') # DDRTree方式 HSMM <- orderCells(HSMM)
三 Monocle2可视化
3.1 基于各种“类型”可视化
有了树结构后，分类颜色是可以自己指定的。而且可以直接使用ggplot2的color设置，是不是觉得多了解一下ggplot2很有必要
colour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB","#98FB98","#F08080","#1E90FF","#7CFC00","#FFFF00", "#808000","#FF00FF","#FA8072","#7B68EE","#9400D3","#800080","#A0522D","#D2B48C","#D2691E","#87CEEB","#40E0D0","#5F9EA0", "#FF1493","#0000CD","#008B8B","#FFE4B5","#8A2BE2","#228B22","#E9967A","#4682B4","#32CD32","#F0E68C","#FFFFE0","#EE82EE", "#FF6347","#6A5ACD","#9932CC","#8B008B","#8B4513","#DEB887") a1 <- plot_cell_trajectory(HSMM, color_by = "seurat_clusters") + scale_color_manual(values = colour) a2 <- plot_cell_trajectory(HSMM, color_by = "State") + scale_color_manual(values = colour) a3 <- plot_cell_trajectory(HSMM, color_by = "Pseudotime") a4 <- plot_cell_trajectory(HSMM, color_by = "labels") + scale_color_manual(values = colour) (a1 + a2) / (a3 + a4)
分别为根据 seurat cluster ，State ，Pseudotime 和 singleR注释后的cell type 着色。
3.2 分面展示
可以使用分面单独查看各单一celltype的时序状态
#其他的类似 plot_cell_trajectory(HSMM, color_by = "labels") + facet_wrap(~labels, nrow = 2) #设置几行几列展示
3.3 添加“树形图”
基本展示完了，那如何得到文章中的结果呢？
此处需要plot_complex_cell_trajectory函数添加“树形图”即可
p1 <- plot_cell_trajectory(HSMM, x = 1, y = 2, color_by = "labels") + theme(legend.position='none',panel.border = element_blank()) + #去掉第一个的legend scale_color_manual(values = colour) p2 <- plot_complex_cell_trajectory(HSMM, x = 1, y = 2, color_by = "labels")+ scale_color_manual(values = colour) + theme(legend.title = element_blank())
p1 / p2
其中：p1中没有展示legend， p2去掉legend的legend.title 。ggplot2的函数可以无缝链接，不想看一下吗？
ggplot2|theme主题设置，详解绘图优化-“精雕细琢”
