单细胞分析实录(15): 基于monocle2的拟时序分析

关于什么是“拟时序分析”，可以参考本期推送的另一篇推文。这一篇直接演示代码

monocle2这个软件用得太多了，很多文章都是monocle2的图。因为只使用表达矩阵作为输入，相比于其他软件，已经很方便了。不过话说回来，我总感觉这种轨迹推断像玄学，大半的结果是调整出来的/事先已知的，比如拟时序里面的起始状态经常需要自己指定。

在做拟时序之前，最好先跑完seurat标准流程，并注释好细胞类型。我这里使用的数据都已经注释好细胞类型，并且事先已经知道哪一个细胞类型可能是初始状态。

1. 导入数据，创建对象，预处理

library(monocle)
library(tidyverse)

expr_matrix=read.table("count_mat.txt",header = T,row.names = 1,sep = "	",stringsAsFactors = F)
#10X的数据使用UMI count矩阵
cell_anno=read.table("cell_anno.txt",header = T,row.names = 1,sep = "	",stringsAsFactors = F)
gene_anno=read.table("gene_anno.txt",header = T,row.names = 1,sep = "	",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="	",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 = '	', row.names = FALSE, col.names = TRUE)
write.table(states_de, file = "states_de.rds", quote = FALSE, sep = '	', 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()

因水平有限，有错误的地方，欢迎批评指正！
在公众号可以获取本文的测试数据和全部代码。