单细胞分析实录(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 = te