使用CelRank进行单细胞命运映射.

0. TL; DR

计算轨迹推断使得从单细胞RNA测序（scRNA-seq）实验中重建细胞状态动态成为可能。然而，传统的轨迹推断方法需要预先知道生物过程的方向，这在很大程度上限制了其应用，使其主要适用于正常发育等方向明确的场景。

本文介绍了一种名为CellRank的方法，它能够在方向未知的多种场景下（如再生、重编程和疾病）绘制单细胞的命运图谱。CellRank的核心思想是，将基于相似性的轨迹推断的稳健性与RNA velocity提供的方向性信息相结合。它将细胞状态的转变建模为一个马尔可夫链，不仅考虑了细胞命运决定的渐进和随机性，还考虑了RNA velocity向量本身的不确定性。

作者在一系列生物学问题中展示了CellRank的强大能力：在胰腺发育数据中，CellRank自动检测到初始、中间和终末细胞群，准确预测了细胞的命运潜能，并可视化了沿特定谱系连续的基因表达趋势。在带有谱系追踪信息的细胞重编程数据中，CellRank预测的命运概率与实验的“金标准”结果高度一致。在损伤后的肺再生数据中，CellRank预测了一条新的去分化轨迹，并发现了一个先前未知的中间细胞状态，这一发现得到了后续实验的证实。

1. 背景介绍

细胞在发育、再生、重编程和癌症等众多生物过程中会经历状态的转变，而且这种转变通常是高度异步的。scRNA-seq技术成功地捕捉了由此产生的细胞异质性，但由于每个细胞只能被测量一次，细胞间的谱系关系信息丢失了。

为了解决这个问题，研究人员发展了多种方法。实验上，可以使用谱系追踪（lineage tracing）或代谢标记（metabolic labeling）来追踪细胞的演化，但这些方法大多局限于体外应用。计算上，伪时间（pseudotime）轨迹推断方法被广泛使用，它基于“发育上相关的细胞在基因表达上更相似”这一假设来排序细胞。

然而，传统的计算轨迹推断方法通常需要先验知识来确定细胞状态变化的方向，例如需要手动指定一个“起始细胞”，这限制了它们在那些方向不明确的生物过程（如再生或疾病）中的应用。

近年来，RNA velocity的出现部分解决了方向性问题。它通过分析mRNA的剪接与未剪接比例，为每个细胞推断出一个指向其未来状态的“速度向量”。尽管这一方法已被推广到更复杂的场景，但velocity的估计本身充满噪声，而且对高维速度向量的解释大多局限于低维投影，这使得我们难以获得定量的、长期的细胞命运概率。

因此，当前领域迫切需要一种方法，既能利用RNA velocity提供的方向性，又能像轨迹推断方法一样稳健地捕捉全局结构，从而在无需先验知识的情况下，推断出细胞的初始、终末状态以及它们之间的概率性转变关系。

为了应对这一挑战，作者提出了CellRank。CellRank巧妙地将细胞间的相似性与RNA velocity的方向性信息结合起来，构建了一个马尔可夫模型，用以学习细胞状态的有向、概率性转变轨迹。CellRank不仅能够自动识别初始和终末细胞群，还能计算每个细胞到达各个终点的“命运概率”，并考虑了细胞命运决定的随机性和velocity估计的不确定性。

2. CellRank 算法

CellRank算法的核心目标是模拟细胞状态的动态，并自动识别系统的初始、中间和终末状态，最终为每个细胞计算一个全局的命运图谱，即到达每个终末状态的概率。

2.1 结合相似性与方向性的马尔可夫链

CellRank将细胞状态的转变建模为一个马尔可夫链（Markov chain）。在这个链中，每个观测到的细胞都是一个“状态”，而状态之间的“转移概率”则由细胞间的相似性和RNA velocity共同决定。

构建基于相似性的基础图谱

与伪时间方法类似，CellRank首先基于细胞在高维基因表达空间的相似性（通常在PCA空间中计算），构建一个K-近邻图（K nearest neighbor, KNN graph）。这个图定义了细胞状态流形的拓扑结构，即哪些细胞之间可能存在直接的转变关系。

利用RNA Velocity赋予方向性

接下来，CellRank利用RNA velocity为这个无向的KNN图赋予“方向”。对于一个细胞$i$和它的一个邻居$j$，计算细胞$i$的velocity向量$v_i$与它们之间的状态变化向量$\delta_{i,j} = x_j - x_i$之间的一致性（通过皮尔逊相关系数，即向量夹角的余弦值）。如果邻居$j$恰好在$v_i$所指的方向上，那么从$i$到$j$的转移概率就高；反之则低。这个概率可以通过一个softmax函数来计算：

\[p_{i,k} = \frac{e^{\sigma \cdot c_{ik}}}{\sum_l e^{\sigma \cdot c_{il}}}\]

其中，$c_{ik}$是$v_i$和$\delta_{i,k}$的相关性，$\sigma$是一个缩放常数。

构建最终的转移概率矩阵

为了增加模型的稳健性，CellRank将上述基于velocity的转移矩阵$P_v$与一个纯粹基于基因表达相似性（即KNN图的连接性）的转移矩阵$P_s$进行加权平均，得到最终的转移矩阵$P$：

\[P = (1-\lambda)P_v + \lambda P_s\]

其中，$\lambda$是一个权重参数（默认为0.2），用于平衡方向性信息和相似性信息。这个矩阵$P$就是作者们描述细胞动态的马尔可夫链的核心。

2.2 从马尔可夫链中提取生物学洞见

直接分析这个巨大的、充满噪声的转移矩阵$P$是很困难的。因此，CellRank采用了一种“粗粒化”（coarse-graining）的策略，将成千上万的单细胞状态总结为少数几个关键的“宏观状态”（macrostates）。

识别宏观状态

CellRank使用一种名为“广义佩龙聚类分析”（Generalized Perron Cluster Cluster Analysis, GPCCA）的方法来识别宏观状态。这些宏观状态是细胞在表型流形上的一些区域，一旦细胞进入这些区域，就很难再离开。GPCCA通过对转移矩阵$P$进行实舒尔分解（real Schur decomposition），找到其不变子空间，从而将单个细胞“软分配”（soft assignment）到不同的宏观状态。

自动识别初始和终末状态

在得到宏观状态及其之间的转移概率后，CellRank可以自动识别它们的生物学身份：

计算命运概率

一旦确定了终末状态，作者们就可以计算每个细胞到达每一个终末状态的概率。这在数学上对应于一个“吸收概率”（absorption probabilities）问题。对于任意一个细胞，作者们可以计算从它出发的随机游走（random walk）最终被某个终末状态“吸收”的概率。这个问题可以高效地被求解为一个线性系统，从而为每个细胞得到一个关于所有终末状态的命运概率向量。

2.3 应对RNA Velocity的不确定性

RNA velocity的估计本身受到多种生物和技术噪声的影响，因此是不确定的。CellRank通过将velocity向量$v$视为一个随机变量（服从一个以邻近细胞velocity均值为中心的高斯分布），并将这种不确定性传播到下游计算中，来提高模型的稳健性。作者们提供了一种高效的解析近似方法，通过泰勒展开来计算期望的转移概率，从而避免了耗时的蒙特卡洛采样（Monte Carlo sampling）。

2.4 下游分析

基于计算出的命运概率，CellRank提供了一系列强大的下游分析功能：

3. 实验分析

作者在一系列真实的生物学数据集上展示了CellRank的强大功能，涵盖了正常发育、重编程和组织再生等多种场景。

实验一：解析胰腺内分泌谱系的形成

作者将CellRank应用于一个E15.5小鼠胰腺发育的数据集。

实验二：揭示delta细胞分化的驱动基因

delta细胞在这个数据集中非常稀少，其分化机理尚不清楚，scVelo也未能有效捕捉其动力学。CellRank通过其全局分析能力，成功揭示了其分化路径。

实验三：在细胞重编程中预测细胞命运

作者将CellRank应用于一个带有谱系追踪“金标准”的细胞重编程数据集，以验证其预测能力。在这个数据集中，只有约1%的细胞能成功重编程为iEPs，其余则进入“死胡同”（dead-end）状态。

实验四：在肺再生中发现新的去分化轨迹

作者将CellRank应用于一个损伤后的小鼠肺再生数据集，以展示其在方向未知的复杂过程中的发现能力。

这个发现表明，杯状细胞在损伤后可能通过去分化来补充干细胞库，CellRank成功地预测了这一先前未知的生物学过程。

与其他方法的性能比较