近年来在Lancet、BMJ等杂志经常见到利用限制性立方样条来拟合非线性关系。
什么是立方样条?
回归样条(regression spline)本质上是一个分段多项式, 但它一般要求每个分段点上连续并且二阶可导，这样可以保证曲线的平滑性。而限制性立方样条是在回归样条的基础上附加要求：样条函数在自变量数据范围两端的两个区间内为线性函数。
在利用限制性立方样条绘制曲线关系时，通常需要设置样条函数节点的个数(k)和位置(ti)。绝大多数情况下, 节点的位置对限制性立方样条的拟合影响不大, 而节点的个数则决定曲线的形状, 或者说平滑程度。当节点的个数为2时, 得到的拟合曲线就是一条直线，大多数研究者推荐的节点为3-5个。
在《Regression Modeling Strategies》这本书中，Harrell建议节点数为4时，模型的拟合较好，同时可以兼顾曲线的平滑程度和避免过拟合造成的精度降低。而当样本量较大时，例如因变量为未删失的连续变量并且大于100时，5个节点是更好的选择。小样本(如n<30)可以选择3个节点。以下是Harrell推荐的节点数和相应的节点位置，大家可以参考。
案例说明(模拟数据)
目前SAS、STATA、R等软件都可以进行限制性立方样条分析。基于画图的方便，我们以R语言为例进行说明。首先参照rms包，生成一个模拟数据集，包括性别(sex)，年龄(age)以及生存时间(time)和结局变量(death)。
若想分析年龄和生存率之间关系，传统的方法可以在Cox回归中将年龄作为连续变量处理，也可以对年龄进行分组，这样的做法都无法更直观的呈现年龄与死亡风险之间的关联。以下我们用限制性立方样条来分析年龄与死亡风险之间的关系：
可以看到age整体是有意义的(包括线性或者非线性关联)，然后看P-Nonlinear =0.0168<0.05，这里我们可以说年龄与死亡风险之间存在非线性关联。
如果自变量与关注的结局变量存在非线性关系，如何在文章中对结果更详细的描述呢，建议大家可以参照上文中提到的Lancet的文章。
原文链接：http://tecdat.cn/?p=20882
1导言
这篇文章探讨了为什么使用广义相加模型 是一个不错的选择。为此，我们首先需要看一下线性回归，看看为什么在某些情况下它可能不是最佳选择。
2回归模型
假设我们有一些带有两个属性Y和X的数据。如果它们是线性相关的，则它们可能看起来像这样：
a<-ggplot(my_data, aes(x=X,y=Y))+geom_point()+
为了检查这种关系，我们可以使用回归模型。线性回归是一种使用X来预测变量Y的方法。将其应用于我们的数据将预测成红线的一组值：
a+geom_smooth(col="red", method="lm")+
这就是“直线方程式”。根据此等式，我们可以从直线在y轴上开始的位置（“截距”或α）开始描述，并且每个单位的x都增加了多少y（“斜率”），我们将它称为x的系数，或称为β）。还有一点自然的波动，如果没有的话，所有的点都将是完美的。我们将此称为“残差”（ϵ）。数学上是：
或者，如果我们用实际数字代替，则会得到以下结果：
这篇文章通过考虑每个数据点和线之间的差异（“残差）然后最小化这种差异来估算模型。我们在线的上方和下方都有正误差和负误差，因此，通过对它们进行平方并最小化“平方和”，使它们对于估计都为正。这称为“普通最小二乘法”或OLS。
3非线性关系如何？
因此，如果我们的数据看起来像这样，我们该怎么办：
我们刚刚看到的模型的关键假设之一是y和x线性相关。如果我们的y不是正态分布的，则使用广义线性模型 （Nelder＆Wedderburn，1972），其中y通过链接函数进行变换，但再次假设f（y）和x线性相关。如果不是这种情况，并且关系在x的范围内变化，则可能不是最合适的。我们在这里有一些选择：
我们可以使用线性拟合，但是如果这样做的话，我们会在数据的某些部分上面或者下面。
我们可以分为几类。我在下面的图中使用了三个，这是一个合理的选择。同样，我们可能处于数据某些部分之下或之上，而在类别之间的边界附近似乎是准确的。例如，如果x = 49时，与x = 50相比，y是否有很大不同？
我们可以使用多项式之类的变换。下面，我使用三次多项式，因此模型适合：。这些的组合使函数可以光滑地近似变化。这是一个很好的选择，但可能会极端波动，并可能在数据中引起相关性，从而降低拟合度。
请点击输入图片描述
请点击输入图片描述
4样条曲线
多项式的进一步细化是拟合“分段”多项式，我们在数据范围内将多项式链在一起以描述形状。“样条线”是分段多项式，以绘图员用来绘制曲线的工具命名。物理样条曲线是一种柔性条，可以弯曲成形，并由砝码固定。在构造数学样条曲线时，我们有多项式函数，二阶导数连续，固定在“结”点上。
下面是一个ggplot2 对象，该 对象的 geom_smooth 的公式包含ns 函数中的“自然三次样条” 。这种样条曲线为“三次”，并且使用10个结
请点击输入图片描述
请点击输入图片描述
5光滑函数
样条曲线可以是光滑的或“摇摆的”，这可以通过改变节点数（k）或使用光滑惩罚γ来控制。如果我们增加结的数目，它将更“摇摆”。这可能会更接近数据，而且误差也会更小，但我们开始“过度拟合”关系，并拟合我们数据中的噪声。当我们结合光滑惩罚时，我们会惩罚模型中的复杂度，这有助于减少过度拟合。
请点击输入图片描述
6广义相加模型（GAM）
广义加性模型（GAM）（Hastie，1984）使用光滑函数（如样条曲线）作为回归模型中的预测因子。这些模型是严格可加的，这意味着我们不能像正常回归那样使用交互项，但是我们可以通过重新参数化作为一个更光滑的模型来实现同样的效果。事实并非如此，但本质上，我们正转向一种模型，如：
请点击输入图片描述
摘自Wood （2017）的GAM的更正式示例 是：
请点击输入图片描述
其中：
μi≡E（Yi），Y的期望
Yi〜EF（μi，ϕi），Yi是一个响应变量，根据均值μi和形状参数ϕ的指数族分布。
Ai是任何严格参数化模型分量的模型矩阵的一行，其中θ为对应的参数向量。
fi是协变量xk的光滑函数，其中k是每个函数的基础。
如果您要建立回归模型，但怀疑光滑拟合会做得更好，那么GAM是一个不错的选择。它们适合于非线性或有噪声的数据。
7 gam拟合
那么，如何 为上述S型数据建立 GAM模型？在这里，我将使用三次样条回归 ：
gam(Y ~ s(X, bs="cr")上面的设置意味着：
s()指定光滑器。还有其他选项，但是s是一个很好的默认选项
bs=“cr”告诉它使用三次回归样条（'basis'）。
s函数计算出要使用的默认结数，但是您可以将其更改为k=10，例如10个结。
8模型输出：
查看模型摘要：
#### Family: gaussian## Link function: identity## Parametric coefficients:## Estimate Std. Error t value Pr(>|t|)## (Intercept) 43.9659 0.8305 52.94 <2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#### Approximate significance of smooth terms:## edf Ref.df F p-value## s(X) 6.087 7.143 296.3 <2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#### R-sq.(adj) = 0.876 Deviance explained = 87.9%## GCV = 211.94 Scale est. = 206.93 n = 300显示了我们截距的模型系数，所有非光滑参数将在此处显示
每个光滑项的总体含义如下。
这是基于“有效自由度”（edf）的，因为我们使用的样条函数可以扩展为许多参数，但我们也在惩罚它们并减少它们的影响。
9检查模型：
该 gam.check() 函数可用于查看残差图，但它也可以测试光滑器以查看是否有足够的结来描述数据。但是如果p值很低，则需要更多的结。
请点击输入图片描述
#### Method: GCV Optimizer: magic## Smoothing parameter selection converged after 4 iterations.## The RMS GCV score gradient at convergence was 1.107369e-05 .## The Hessian was positive definite.## Model rank = 10 / 10#### Basis dimension (k) checking results. Low p-value (k-index<1) may## indicate that k is too low, especially if edf is close to k'.#### k' edf k-index p-value## s(X) 9.00 6.09 1.1 0.9710它比线性模型好吗？
让我们对比具有相同数据的普通线性回归模型：
anova(my_lm, my_gam)## Analysis of Variance Table#### Model 1: Y ~ X## Model 2: Y ~ s(X, bs = "cr")## Res.Df RSS Df Sum of Sq F Pr(>F)## 1 298.00 88154## 2 292.91 60613 5.0873 27540 26.161 <2.2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
我们的方差分析函数在这里执行了f检验，我们的GAM模型明显优于线性回归。
11小结
所以，我们看了什么是回归模型，我们是如何解释一个变量y和另一个变量x的。其中一个基本假设是线性关系，但情况并非总是这样。当关系在x的范围内变化时，我们可以使用函数来改变这个形状。一个很好的方法是在“结”点处将光滑曲线链接在一起，我们称之为“样条曲线”
我们可以在常规回归中使用这些样条曲线，但是如果我们在GAM的背景中使用它们，我们同时估计了回归模型以及如何使我们的模型更光滑。
上面的示例显示了基于样条的GAM，其拟合度比线性回归模型好得多。
12参考：
NELDER, J. A. &WEDDERBURN, R. W. M. 1972. Generalized Linear Models. Journal of the Royal Statistical Society. Series A (General), 135, 370-384.
HARRELL, F. E., JR. 2001. Regression Modeling Strategies, New York, Springer-Verlag New York.
请点击输入图片描述
最受欢迎的见解
1.R语言多元Logistic逻辑回归 应用案例
2.面板平滑转移回归(PSTR)分析案例实现
3.matlab中的偏最小二乘回归（PLSR）和主成分回归（PCR）
4.R语言泊松Poisson回归模型分析案例
5.R语言回归中的Hosmer-Lemeshow拟合优度检验
6.r语言中对LASSO回归，Ridge岭回归和Elastic Net模型实现
7.在R语言中实现Logistic逻辑回归
8.python用线性回归预测股票价格
9.R语言如何在生存分析与Cox回归中计算IDI，NRI指标
伪时间是衡量单个细胞在细胞分化等过程中取得了多大进展的指标。在许多生物学过程中，细胞并不是完全同步的。在细胞分化等过程的单细胞表达研究中，捕获的细胞在分化方面可能分布广泛。也就是说，在同一时间捕获的细胞群中，有些细胞可能已经很长时间了，而有些细胞甚至还没有开始这个过程。当您想要了解在细胞从一种状态转换到另一种状态时所发生的调节更改的顺序时，这种异步性会产生主要问题。跟踪同时捕获的细胞间的表达可以产生对基因动力学一个大致的认识，该基因表达的明显变异性将非常高。Monocle根据每个cell在学习轨迹上的进展对其进行排序，从而缓解了由于异步而产生的问题。Monocle不是跟踪表达式随时间变化的函数，而是跟踪沿轨迹变化的函数，我们称之为伪时间。伪时间是一个抽象的分化单位：它只是一个cell到轨迹起点的距离，沿着最短路径测量。轨迹的总长度是由细胞从起始状态移动到结束状态所经历的总转录变化量来定义的。
还是老规矩，跟着官网学习Monocle2的使用方法。官网如下：
单细胞基因表达研究使人们能够在复杂的生物过程和高度异质性的细胞群中描述转录调控。这些研究有助于发现识别特定亚型细胞的基因，或标记生物过程中的中间状态，以及在两种不同的细胞命运之间过渡态的基因。在许多单细胞研究中，单个细胞是以不同步的方式执行基因表达过程。实际上，每个细胞都是正在研究的转录过程中的一个瞬间。Monocle包是分析单细胞测序的工具。 Monocle引入了在伪时间(拟时间)内对单个细胞排序的策略，利用单个细胞的非同步进程，将它们置于与细胞分化等生物学过程相对应的轨迹上。****Monocle利用先进的机器学习技术(反向图嵌入)从单细胞数据中学习显式的主图来对细胞进行排序，这可以强大而准确地解决复杂的生物过程。 Monocle也可以进行聚类(即使用t-SNE和密度峰值聚类)。Monocle也可以进行差异基因表达测试，使人们能够识别在不同状态下差异表达的基因，沿着生物过程以及不同的细胞命运时基因表达的变化。Monocle是专为单细胞RNA-Seq研究设计的，但也可以用于其他分析。
Monocle 2包括新的和改进的算法，用于对细胞进行分类和计数，执行细胞亚群之间的差异表达分析，以及重建细胞轨迹。Monocle 2也被重新设计，可以很好地完成包含数万个或更多细胞的大型单细胞RNA-Seq数据分析。
Monocle perform three main types of analysis:
首先，Monocle 2使用一种简单的、无偏的和高度可扩展的统计程序来选择具有轨迹进展特征的基因。然后，它采用了一类流形学习算法，旨在在高维单细胞RNA-seq数据中嵌入一个主图。以前的方法是通过启发式分析细胞之间的成对距离来推断分支结构，而Monocle 2可以使用这张图来直接识别发育的命运决定。我们已经通过广泛的基准测试证明，Monocle 2优于其他工具，如Wishbone，而不需要用户指定轨迹的结构。
Monocle 2的算法是再2017年发表在nature methods 上，文章链接如下
文章中的核心理论为：每个细胞都可以表示为高维空间中的一个点，在高维空间中，每个维对应着一个有序基因的表达水平。高维数据首先通过几种降维方法，如PCA(默认)、扩散映射等，投射到低维空间。 Monocle 2然后在自动选择的一组数据质心上构造一棵生成树。 然后， 该算法将细胞移动到它们最近的树的顶点，更新顶点的位置以适应细胞，学习新的生成树，并迭代地继续这个过程，直到树和细胞的位置已经收敛。在这个过程中，Monocle 2保持了高维空间和低维空间之间的可逆映射，从而既学习了轨迹，又降低了数据的维数。****一旦Monocle 2学会了树，用户就会选择一个tip作为根。计算每个单元的伪时间作为其沿树到根的测地线距离，并根据主图自动分配其分枝。 因为monocle2学习树结构，与其他方法相比，分支结构自动出现。当它更新细胞位置并细化树时，monocle2简化了轨迹的结构，修剪了小的分支，这样最终的轨迹只保留了描述细胞状态显著差异的分支。
Monocle需要在R语言环境下运行。官网给出的都是R-3.4比较老的版本，现在都是R-3.6或者R-4.0以上版本。我们可以直接在Bioconductor中安装Monocle 2 。
查看系统中安装的这个包的版本文档
当然也可以在GitHub上安装最新版的Monocle2。
有时我们添加的特性需要你安装某些额外的软件包。当您尝试上面的命令时，可能会看到错误。您可以通过输入(例如)来在错误消息中安装包。
最后检查一下是否安装成功
方法一: 将Seurat object中数据提取来创建
FPKM/TPM值通常是对数正态分布的，而UMIs或读计数使用负二项更好地建模。要处理计数数据，请将负二项分布指定为newCellDataSet的expressionFamily参数：
稀疏矩阵用negbinomial.size()，
FPKM值用tobit()，
logFPKM值用gaussianff()
方法二: 直接读取表达矩阵来创建
方法三: 将Seurat对象直接转化为CellDataSet对象
如果我们想要从Seurat对象或SCESet中导入所有的插槽，我们可以设置参数'import_all'为TRUE。#(默认为FALSE或只保留最小数据集)
size facotr帮助我们标准化细胞之间的mRNA的差异。
离散度值可以帮助我们进行后续的差异分析。
大多数单细胞工作流程至少会包含一些由死细胞或空孔组成的库。同样重要的是要删除doublets:由两个或多个细胞意外生成的库。这些细胞可以破坏下游步骤，如伪时间排序或聚类。要知道一个特定的基因有多少个表达，或者一个给定的细胞有多少个基因表达，通常是很方便的。Monocle提供了一个简单的函数来计算这些统计数据。
Monocle官网教程提供了4个分类方法：
1、Classifying cells by type
Monocle提供了一个简单的系统，根据您选择的标记基因的表达来标记细胞。您只需提供一组函数，Monocle就可以使用这些函数对每个单元格进行注释。例如，您可以为几种单元格类型中的每种类型提供一个函数。这些函数接受每个单元格的表达式数据作为输入，并返回TRUE以告诉Monocle某个单元格满足函数定义的条件。所以你可以有一个功能，对表达成肌细胞特异性基因的细胞来说是真的，另一个功能对成纤维细胞特异性基因来说是真的，等等。
2、Clustering cells without marker genes
Monocle提供了一种算法，可以用来计算“未知”单元格的类型。该算法在函数clusterCells中实现，根据全局表达式轮廓将单元格分组在一起。这样的话，如果你的细胞表达了许多特定于成肌细胞的基因，但恰好缺少MYF5，我们仍然可以识别它为成肌细胞。clusterCells可以在无监督的方式下使用，也可以在半监督的“”模式下使用，这允许使用一些专家知识来辅助算法。我们先来看看无监督模式。
第一步是决定哪些基因用于细胞聚类。我们可以使用所有的基因，但是我们会包含很多基因，这些基因的表达水平不足以提供有意义的信号。包括它们只会给系统增加噪音。我们可以根据平均表达水平筛选基因，我们还可以选择在细胞间异常可变的基因。这些基因往往对细胞状态有很高的信息量。
setOrderingFilter函数标记了在随后对clusterCells的调用中用于聚类的基因，尽管我们可以根据需要提供其他基因列表。plot_ordering_genes函数显示了基因表达的变异性（离散性）如何依赖于细胞间的平均表达。红线表示基于这种关系的单片体对色散的预期。我们标记用于聚类的基因显示为黑点，而其他基因显示为灰点。
Monocle使用t-SNE来给细胞聚类可视化
指定10个cluster
3、Clustering cells using marker genes
首先，我们将选择一组不同的基因用于细胞聚类。在我们挑选高表达和高度可变的基因之前。现在，我们将挑选与标记共变异的基因。从某种意义上说，我们将建立一个大的基因列表来作为标记，这样即使一个细胞没有MYF5，它也可以被其他基因识别为成肌细胞。
函数markerDiffTable接受CellDataSet和CellTypeHierarchy，并根据提供的函数将所有单元格分类为类型。然后在识别不同基因类型之间差异表达的基因之前，它去除了所有“未知”和“不明确”的功能。同样，您可以提供要从该测试中排除的影响的剩余模型。然后函数返回一个测试结果的数据帧，您可以使用它来选择要用于聚类的基因。通常情况下，最好挑选最适合每种细胞类型的前10或20个基因。这确保了聚类基因不受一种细胞类型标记的支配。如果可能的话，通常需要为每种类型创建一个平衡的标记面板。Monocle提供了一个方便的功能，可以根据基因在每种类型中的表达受限程度来对基因进行排序。
4、Imputing cell type
注意，我们已经减少了成肌细胞群中“污染”成纤维细胞的数量，反之亦然。但是“未知”细胞呢？如果为clusterCells提供CellTypeHierarcy，Monocle将使用它对整个簇进行分类，而不仅仅是单个细胞。基本上，clusterCells的工作原理与以前完全相同，只是在构建集群之后，它会统计每个集群中每种细胞类型的频率。当一个簇由超过一定百分比（在本例中为10%）的特定类型组成时，簇中的所有单元格都被设置为该类型。如果一个集群由多个单元类型组成，则整个集群都被标记为“不明确”。如果没有高于阈值的单元类型，则集群被标记为“未知”。因此，Monocle可以帮助您在缺少标记数据的情况下估算每个细胞的类型。
在发育过程中，为了对刺激做出反应，在整个生命周期中，细胞从一种功能“状态”过渡到另一种功能“状态”。处于不同状态的细胞表达不同的基因，产生蛋白质和代谢物的动态重复，从而完成它们的工作。当细胞在不同的状态间移动时，会经历一个转录重组的过程，一些基因被沉默，而另一些则被激活。这些瞬态通常很难描述，因为在更稳定的端点状态之间净化细胞可能很困难或不可能。单细胞RNA-Seq可以让你看到这些状态而不需要纯化。然而，要做到这一点，我们必须确定每个细胞的可能状态范围。
Monocle介绍了利用RNA-Seq进行单细胞轨迹分析的策略。Monocle不是通过实验将细胞净化成离散状态，而是使用一种算法来学习每个细胞作为动态生物过程的一部分必须经历的基因表达变化序列。一旦掌握了基因表达变化的整体“轨迹”，Monocle就可以将每个细胞放置在轨迹中的适当位置。然后，你可以使用Monocle的差异分析工具来寻找在轨迹过程中受调控的基因，如“寻找随假时间变化的基因”一节所述。如果这个过程有多个结果，Monocle将重建一个“分支”轨迹。这些分支与细胞的“决定”相对应，Monocle提供了强有力的工具来识别受它们影响并参与制造它们的基因。在“分析单细胞轨迹中的分支”部分中，可以看到如何分析分支。Monocle依靠一种称为反向图嵌入的机器学习技术来构造单细胞轨迹.
Step 1: 选择定义过程的基因
推断单细胞轨迹是一个机器学习问题。第一步是选择Monocle将用作机器学习方法输入的基因。这叫做特征选择，它对轨迹的形状有很大的影响。在单细胞RNA-Seq中，低水平表达的基因通常非常嘈杂，但有些基因包含有关细胞状态的重要信息。单核细胞通过检测这些基因在细胞群体中的表达模式来对细胞进行排序。Monocle寻找以“有趣的”（即不只是嘈杂的）方式变化的基因，并用它们来构造数据。Monocle为您提供了多种工具来选择基因，这些基因将产生一个健壮、准确和具有生物学意义的轨迹。你可以使用这些工具来执行一个完全“无监督”的分析，在这个分析中，Monocle不知道你认为哪个基因是重要的。或者，你可以利用marker gene知识来定义生物学进展，从而形成Monocle的轨迹。我们认为这种模式是“半监督”的分析。
Monocle官网教程提供了4个选择方法：
前三种都是无监督分析方法，细胞发育轨迹生成完全不受人工干预；最后一种是半监督分析方法，可以使用先验知识辅助分析。
理想情况下，我们希望尽可能少地使用正在研究的系统生物学的先验知识。我们希望从数据中发现重要的排序基因，而不是依赖于文献和教科书，因为这可能会在排序中引入偏见。我们将从一种更简单的方法开始，但是我们通常推荐一种更复杂的方法，称为“dpFeature”。
Step 2: 降维
下一步，我们将把空间缩小到一个二维空间，当Monocle对细胞进行排序时，我们将能够很容易地可视化和解释这些空间。
Step 3: 按照轨迹排序细胞
如果你的树上有一堆状态，那就很难弄清楚每个状态都落在树上的什么地方了。有时，可以方便地“分面”轨迹图，以便更容易地查看每个状态的位置.
如果你没有时间序列，你可能需要利用你对系统的生物学知识，根据特定标记基因的表达位置来设置根。例如，在这个实验中，一个高度增殖的祖细胞群正在产生两种类型的有丝分裂后细胞。因此，根部应该有表达高水平增殖标记的细胞。我们可以使用抖动图来判断哪种状态对应于快速扩散：
官方给出的差异分析有三大方法，我们重点关注第三个：****根据伪时间功能寻找差异基因
1、Basic Differential Analysis
2、Finding Genes that Distinguish Cell Type or State
3、Finding Genes that Change as a Function of Pseudotime
Monocle的主要工作是通过生物过程（如细胞分化）将细胞按顺序排列，而不知道要提前查看哪些基因。一旦这样做了，你就可以分析细胞，找到随着细胞进展而变化的基因。例如，你可以发现随着细胞“成熟”而显著上调的基因。让我们来看看一组对肌肉生成很重要的基因：
这个sm.ns函数指出Monocle应该通过表达式值拟合自然样条曲线，以帮助它将表达式的变化描述为进程的函数。让我们加入基因注释，这样就很容易看出哪些基因是重要的。
** 研究时间序列基因表达研究时出现的一个常见问题是：“哪些基因遵循相似的动力学趋势”？Monocle可以通过对具有相似趋势的基因进行分组来帮助你回答这个问题，所以你可以分析这些基因组，看看它们有什么共同点。Monocle提供了一种方便的方法来可视化所有假时间依赖性基因。函数plot_pseudotime_heatmap采用CellDataSet对象（通常只包含重要基因的子集），并生成平滑的表达曲线，非常类似于plot_genes_in_pseudotime。然后，它对这些基因进行聚类，并使用pheatmap软件包进行绘图。这让你可以想象出基因模块在不同的时间内共同变化。**
通常，单细胞轨迹包括分支。这些分支的产生是因为细胞执行不同的基因表达程序。在发育过程中，当细胞做出命运选择时，分支出现在轨迹中：一个发育谱系沿着一条路径前进，而另一个谱系产生第二条路径。Monocle包含分析这些分支事件的广泛功能。
以Steve Quake's的实验室进行的实验为例，Barbara Treutlein团队从发育中的小鼠肺中捕获了细胞。他们在发育早期捕获细胞，后来当肺中同时含有两种主要类型的上皮细胞（AT1和AT2）时，细胞就要决定变成AT1或AT2。Monocle可以将这个过程重构为一个分支轨迹，允许您对决策点进行非常详细的分析。下图显示了使用部分数据重建的Monocle轨迹。有一个分支，标记为“1”。当细胞从树的左上角经过树枝的早期发育阶段时，哪些基因会发生变化？哪些基因在分枝间有差异表达？为了回答这个问题，Monocle为您提供了一个特殊的统计测试：分支表达式分析建模，或BEAM。BEAM(Branched expression analysis modeling)是一种统计方法，用于寻找以依赖于分支的方式调控的基因。
BEAM进行统计分析
该热图显示的是同一时间点两个谱系的变化。热图的列是伪时间的点，行是基因。从热图中间往右读，是伪时间的一个谱系；往左是另一个谱系。基因是被按照等级聚类的，所以你看到的基因表达模式和谱系的表达模式是非常相似的。
We can plot a couple of these genes, such as Pdpn and Sftpb (both known markers of fate in this system), using the plot_genes_branched_pseudotime function, which works a lot like the plot_genes_in_pseudotime function, except it shows two kinetic trends, one for each lineage, instead of one.
相关资料：
官网
monocle2 拟时间分支点分析结果解读
单细胞分析实录(15): 基于monocle2的拟时序分析
Monocle2包学习笔记（四）：Differential Expression Analysis
