单细胞分析实录(15): 基于monocle2的拟时序分析
关于什么是“拟时序分析”,可以参考本期推送的另一篇推文。这一篇直接演示代码
monocle2这个软件用得太多了,很多文章都是monocle2的图。因为只使用表达矩阵作为输入,相比于其他软件,已经很方便了。不过话说回来,我总感觉这种轨迹推断像玄学,大半的结果是调整出来的/事先已知的,比如拟时序里面的起始状态经常需要自己指定。
在做拟时序之前,最好先跑完seurat标准流程,并注释好细胞类型。我这里使用的数据都已经注释好细胞类型,并且事先已经知道哪一个细胞类型可能是初始状态。
1. 导入数据,创建对象,预处理
library(monocle) library(tidyverse) expr_matrix=read.table("count_mat.txt",header = T,row.names = 1,sep = "\t",stringsAsFactors = F) #10X的数据使用UMI count矩阵 cell_anno=read.table("cell_anno.txt",header = T,row.names = 1,sep = "\t",stringsAsFactors = F) gene_anno=read.table("gene_anno.txt",header = T,row.names = 1,sep = "\t",stringsAsFactors = F) #“gene_short_name”为列名 pd=new("AnnotatedDataFrame",data = cell_anno) fd=new("AnnotatedDataFrame",data = gene_anno) test=newCellDataSet(as(as.matrix(expr_matrix),'sparseMatrix'),phenoData = pd,featureData = fd) #大数据集使用稀疏矩阵,节省内存,加快运算 test <- estimateSizeFactors(test) test <- estimateDispersions(test) test=detectGenes(test,min_expr = 0.1) #计算每个基因在多少细胞中表达
2. 选择基因
选择研究的生物学过程涉及到的基因集,这一步对于轨迹形状的影响很大。 可以选择数据集中的高变基因,或者是在seurat中分析得到的marker基因列表。如果是时间序列数据,可以直接比较时间点选差异基因。总之选择的基因要含有尽可能多的信息。 我这里直接用的各种亚群差异基因的集合
marker_gene=read.table("seurat_marker_gene.txt",header=T,sep="\t",stringsAsFactors=F) test_ordering_genes=unique(marker_gene$gene) test=setOrderingFilter(test,ordering_genes = test_ordering_genes) #指明哪些基因用于后续的聚类/排序
3. 推断轨迹,并按照拟时序给细胞排序
test=reduceDimension(test,reduction_method = "DDRTree",max_components = 2, norm_method = 'log',residualModelFormulaStr = "~num_genes_expressed") #residualModelFormulaStr减少其他因素的影响,比如不同样本、不同批次 test=orderCells(test)
4. 绘图
plot_cell_trajectory(test,color_by = "celltype") ggsave("celltype.pdf",device = "pdf",width = 8,height = 9,units = c("cm")) plot_cell_trajectory(test,color_by = "State") ggsave("State.pdf",device = "pdf",width = 8,height = 9,units = c("cm")) plot_cell_trajectory(test,color_by = "Pseudotime") ggsave("Pseudotime.pdf",device = "pdf",width = 8,height = 9,units = c("cm")) plot_cell_trajectory(test,color_by = "celltype")+facet_wrap(~celltype,nrow=1) ggsave("celltypeb.pdf",device = "pdf",width = 21,height = 9,units = c("cm"))
有时候(大多数时候),拟时序的方向或是根节点弄错了,还需要手动更改
test=orderCells(test,root_state = 3) plot_cell_trajectory(test,color_by = "Pseudotime") ggsave("Pseudotime.pdf",device = "pdf",width = 8,height = 9,units = c("cm"))
5. 找差异基因
这里是指找随拟时序变化的差异基因,以及不同state之间的差异基因。这两个都是monocle里面的概念,与seurat里面找的差异基因不同。 如果在做monocle2的时候,想展示这种差异基因,就需要做这一步。
expressed_genes=row.names(subset(fData(test),num_cells_expressed>=10)) #在部分基因里面找 pseudotime_de <- differentialGeneTest(test[expressed_genes,], fullModelFormulaStr = "~sm.ns(Pseudotime)") pseudotime_de <- pseudotime_de[order(pseudotime_de$qval), ] states_de <- differentialGeneTest(test[expressed_genes,], fullModelFormulaStr = "~State") states_de <- states_de[order(states_de$qval), ] saveRDS(test, file = "test_monocle.rds") write.table(pseudotime_de, file = "pseudotime_de.rds", quote = FALSE, sep = '\t', row.names = FALSE, col.names = TRUE) write.table(states_de, file = "states_de.rds", quote = FALSE, sep = '\t', row.names = FALSE, col.names = TRUE)
6. 分支点的分析
解决的问题:沿着拟时序的方向,经过不同的branch,发生了哪些基因的变化? 经常在文献里面看到的monocle2热图就是这种分析。
BEAM_res=BEAM(test,branch_point = 1,cores = 1) #会返回每个基因的显著性,显著的基因就是那些随不同branch变化的基因 #这一步很慢 BEAM_res=BEAM_res[,c("gene_short_name","pval","qval")] saveRDS(BEAM_res, file = "BEAM_res.rds")
画图
tmp1=plot_genes_branched_heatmap(test[row.names(subset(BEAM_res,qval<1e-4)),], branch_point = 1, num_clusters = 4, #这些基因被分成几个group cores = 1, branch_labels = c("Cell fate 1", "Cell fate 2"), #hmcols = NULL, #默认值 hmcols = colorRampPalette(rev(brewer.pal(9, "PRGn")))(62), branch_colors = c("#979797", "#F05662", "#7990C8"), #pre-branch, Cell fate 1, Cell fate 2分别用什么颜色 use_gene_short_name = T, show_rownames = F, return_heatmap = T #是否返回一些重要信息 ) pdf("branched_heatmap.pdf",width = 5,height = 6) tmp1$ph_res dev.off()
monocle2的热图怎么看
- 从你想研究的节点(第4步图中的1节点)到根(一般是人为指定,图中是Ucell对应的状态)都是pre-branch。
- 然后沿着拟时序的方向,State 1对应的枝是fate 1(对应图中的Gcell),State 2对应的枝是fate 2(对应Ccell)。软件作者是这么说的:Cell fate 1 corresponds to the state with small id while cell fate 2 corresponds to state with bigger id
- 热图中的基因是BEAM函数找到的branch特异的基因。从热图中间向两边看,能看到不同的fate/branch基因的动态变化。热图中列是拟时序,行是基因,热图中间是拟时序的开始。
返回的列表tmp1包含热图对象,热图的数值矩阵,热图主题颜色,每一行的注释(基因属于哪一个group)等信息。 根据行注释提取出基因group,可以去做富集分析,后期加到热图的旁边。像这样:
gene_group=tmp1$annotation_row gene_group$gene=rownames(gene_group) library(clusterProfiler) library(org.Hs.eg.db) allcluster_go=data.frame() for (i in unique(gene_group$Cluster)) { small_gene_group=filter(gene_group,gene_group$Cluster==i) df_name=bitr(small_gene_group$gene, fromType="SYMBOL", toType=c("ENTREZID"), OrgDb="org.Hs.eg.db") go <- enrichGO(gene = unique(df_name$ENTREZID), OrgDb = org.Hs.eg.db, keyType = 'ENTREZID', ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.2, readable = TRUE) go_res=go@result if (dim(go_res)[1] != 0) { go_res$cluster=i allcluster_go=rbind(allcluster_go,go_res) } } head(allcluster_go[,c("ID","Description","qvalue","cluster")])
也可以换一种方式来展示具体的基因,这些基因可以是:
- 随拟时序、state而改变的基因(第5步)
- 细胞亚群的marker基因(seurat得到的差异基因)
- 分支点分析得到的分支特异的基因(第6步BEAM函数得到的基因)
test_genes=c("TFF3","GUCA2B") pdf("genes_branched_pseudotime.pdf",width = 9,height = 4) plot_genes_branched_pseudotime(test[test_genes,], branch_point = 1, color_by = "celltype", cell_size=2, ncol = 2) dev.off()
因水平有限,有错误的地方,欢迎批评指正! 在公众号可以获取本文的测试数据和全部代码。
- 基于堆叠卷积长短期神经网络【CNNLSTM】模型的时序数据预测分析
- 基于FPGA EEPROM读写实现及IIC总线协议和时序分析
- 基于visual c++之windows核心编程代码分析(15)使用Mutex同步线程
- 单细胞分析实录(8): 展示marker基因的4种图形(一)
- 基于LSTM的【气象数据+发电数据】多步时序数据建模预测分析实战
- 基于Verilog的VGA驱动设计(一)VGA时序分析
- 单细胞分析实录(5): Seurat标准流程
- 单细胞分析实录(6): 去除批次效应/整合数据
- 时序分析:DTW算法(基于模板)
- 基于建立/保持时间等的参数化时序分析
- 单细胞分析实录(2): 使用Cell Ranger得到表达矩阵
- 单细胞分析实录(9): 展示marker基因的4种图形(二)
- 单细胞分析实录(3): Cell Hashing数据拆分
- 基于某知名招聘网站的上海财务岗位数据分析
- OLAP引擎:基于Druid组件进行数据统计分析
- 基于Linux的集群系统(三) 关键技术分析
- 基于数据库的Struts Menu动态菜单分析
- HGSM——基于层级结构图的相似度分析
- 基于TCP网络通信的自动升级程序源码分析--生成升级文件相关的配置文件
- [个人推荐] Linux poll机制分析(基于内核3.10.0)