模拟细胞的分化轨迹的软件，最常用的软件为Monocle2。该算法不仅能模拟细胞的发育轨迹，同时也能对细胞进行聚类（t-SNE）。通过聚类获得不同状态下的差异基因，分析影响分支形成的关键基因及其功能，对研究相关生物学问题有指导性的作用。
Monocle2主要基于关键基因的表达模式，通过学习每个细胞必须经历的基因表达变化的序列，根据拟时间值中对单个细胞进行排序，模拟出时间发育过程的动态变化。而这个排序技术表现是一种在低维空间排布高维数据的降维技术。（具体描叙请参考： https://www.meiwen.com.cn/subject/zlgsvqtx.html ），目前monocle已经升级至monocle3，但在结果分析上monocle3的可读性不如monocle2，因此本研究中主推monocle2这个版本。
在做拟时序分析之前，先要明白几个研究目的：
先保证服务器上装有R ，我这里安装的版本是version 3.5.1。
对于输入文件，有3种类型的格式：
一种是单细胞运行获得的3个结果文件。格式如下表所示：
准备需要进行拟时间分析数据的三个文件：
1)表达量文件：exprs：基因在所有细胞中的count值矩阵。
格式示例：
注意：该文件的表头必须以“%%**\n%”的形式出现，否则就会出现报错。
2)表型文件 phenoData:
格式示例，该文件即为每个细胞的barcode信息。：
由于seurat软件包的升级，大量的函数命名发生了变化，对2和3这两种类型的版本的操作也是不一样。
Seurat2.4的版本：
Seurat3.0及更高版本：
如果你没有3个标准格式的文件，也没有生成rds文件，仅仅只有表达量矩阵，具体原理和有rds文件一致，可按照如下方法进行处理。
如果你采用的是第二种输入数据进行分析，即使你原先的Seurat对象已经归一化过了，官方推荐在monocle中重新进行归一化处理。
在进行降维聚类之前，我们先获得高表达的基因集作为每个聚类的用来order的Feature基因。当然我们可以使用所有的基因，但一些表达量特别低的基因提供的聚类信号往往会制造分析噪音，Feature基因的选择性很多，一种是我们可以根据基因的平均表达水平来进行筛选，另外我们也可以选择细胞间异常变异的基因。这些基因往往能较好地反映不同细胞的状态。以平均表达量高于0.1来进行筛选。
绘制用于order基因和非order的平均表达量的分布点图。这里不做展示。
对细胞进行聚类，在Seurat中采用的是分辨率来确定cluster的数目。而monocle中可以直接指定聚类数目。主要指出的这里所聚类获得的cluster数目比我们赋值的要少一个。即当num_clusters=3时，你只获得了2个cluster。具体解释有些难懂，在这里不做过多的解释。
拟时间分析不仅是要对分析的细胞群进行排序，还要获得觉得细胞排序的关键基因集。这种基因集有两种情况，在有先验知识的情况下，我们可以从系统生物学的角度获得一系列与细胞发育相关的基因集，从而对细胞进行排序，这种方式是最为保险的，但局限性是对未知的发育情况不能进行分析。另外一种情况就是使用无监督聚类方式计算关键基因集。接下来我们采用differentialGeneTest方式获得clustering_DEG_genes（与聚类相关的差异基因集）
上一步过程是对所有的细胞进行无监督训练，运行时间与细胞数和基因数相关，一般会花很长的时间。可以根据cores的数目进行并行。
differentialGeneTest这个函数测试每个基因的差异表达，取决于伪时间或根据指定的其他协变量。 “ differentialGeneTest”是Monocle的主要差异分析常规， 它接受一个CellDataSet和两个模型公式作为输入，指定由实现的广义谱系模型“ VGAM”包。 也就是说我们可以根据指定’~cluster’或者拟时间值来获得差异基因。差异基因的结果如下表所示：
在这个表格中，我们先看一下表头对应的关键列。第一列没有列名，为基因名称。pval,qval 为差异基因的显著检验打分。num_cells_expressed为表达这个基因的细胞数。use_for_ordering该基因是否是作为排序的差异基因。
根据排序好细胞进行结果可视化。
命令行如下所示：
上述部分结果如下图所示（不包括分面图）
上图表示在主成分中的细胞聚类分布的tsne图。不同颜色代表细胞群的不同细胞命运的分支。
上图表示依据seurat的cluster ID映射到拟时间的排序上。
接下来可视化展示细胞state相关的差异基因的表达量分布情况，可以根据这些基因来确定细胞的发育方向，下图仅展示qvalue值排在前6个基因。横纵左边意义见坐标轴的描述。
在monocle中的分析中发现，细胞群能从一个轨迹分叉成不同的分支点，在科研中，我们会比较关注发生分支的原因是什么，即分析分支点之间的差异。Monocle采用分支表达式分析建模，主要是BEAM函数，可以将分叉过程重构为一个分支轨迹，从而分析不同细胞命运下的差异。
命令行如下：
举个例子：我们分析branch_point = 1这个分支处的细胞命名分叉是如何进行的。
即下图中所示。下图只有1个分支点，即分析state1,2,3 这三个state的差异。
对影响分析的基因根据qvalue进行排序。
绘制与分支相关的基因热图。
关键参数为subset (BEAM_res, qval <1e-4))，挑选基因进行热图绘制，也可以设置成其他的阈值。
branch_point=1，分支点选为1。num_clusters=4，将基因根据表达相似性分成4个模块。结果如下图所示：
-getAnywhere(tobit)-help.search("funciton name") / ??function name / getAnywhere(function name)
-A single object matching ‘tobit’ was found.It was found in the following places namespace:VGAM(find namespace)
-HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix),
-+phenoData = pd,
-+featureData = fd,
-+lowerDetectionLimit=0.1,
-+expressionFamily=VGAM:::tobit(Lower=0.1))
OK!
1 If you have no clue about the package, you can use findFn in the sos package as explained in this answer .
2 RSiteSearch("some.function") or searching with rdocumentation or rseek are alternative ways to find the function.
