使用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估计的不确定性。
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可以自动识别它们的生物学身份：
- 终末状态 (Terminal states)：这些宏观状态具有很高的自我转移概率，即细胞一旦进入就倾向于停留在那里。CellRank通过一个“稳定性指数”（stability index, SI）来量化这一点，SI高的宏观状态被识别为终末状态。
- 初始状态 (Initial states)：这些宏观状态很难被再次访问。CellRank通过计算马尔可夫链的粗粒度稳态分布（coarse-grained stationary distribution, CGSD）来识别它们。CGSD值最小的宏观状态被认为是初始状态。
- 中间状态 (Intermediate states)：既非初始也非终末的宏观状态被归为中间状态。
计算命运概率
一旦确定了终末状态，作者们就可以计算每个细胞到达每一个终末状态的概率。这在数学上对应于一个“吸收概率”（absorption probabilities）问题。对于任意一个细胞，作者们可以计算从它出发的随机游走（random walk）最终被某个终末状态“吸收”的概率。这个问题可以高效地被求解为一个线性系统，从而为每个细胞得到一个关于所有终末状态的命运概率向量。
2.3 应对RNA Velocity的不确定性
RNA velocity的估计本身受到多种生物和技术噪声的影响，因此是不确定的。CellRank通过将velocity向量$v$视为一个随机变量（服从一个以邻近细胞velocity均值为中心的高斯分布），并将这种不确定性传播到下游计算中，来提高模型的稳健性。作者们提供了一种高效的解析近似方法，通过泰勒展开来计算期望的转移概率，从而避免了耗时的蒙特卡洛采样（Monte Carlo sampling）。
2.4 下游分析
基于计算出的命运概率，CellRank提供了一系列强大的下游分析功能：
- 基因表达趋势可视化：结合任何一种伪时间排序（CellRank可以自动识别起始细胞来锚定伪时间），CellRank可以绘制特定谱系（lineage）的基因表达趋势。与传统方法不同，CellRank通过命运概率对每个细胞进行加权，从而能够精确地描绘出基因在不同分化路径上的特异性表达模式。
- 发现谱系驱动因子：通过计算基因表达与特定命运概率的相关性，可以识别出那些在特定谱系中起关键调控作用的“驱动基因”（driver genes）。
3. 实验分析
作者在一系列真实的生物学数据集上展示了CellRank的强大功能，涵盖了正常发育、重编程和组织再生等多种场景。
实验一：解析胰腺内分泌谱系的形成
作者将CellRank应用于一个E15.5小鼠胰腺发育的数据集。
- 图a：展示了该数据集的UMAP图和scVelo预测的速度流，可以看到从内分泌祖细胞（EPs）到alpha, beta, epsilon, delta四种内分泌细胞的明显分化趋势。
- 图b, c, d：CellRank的分析结果。
- 图b, c：CellRank自动识别出12个宏观状态，并正确地将alpha, beta, epsilon细胞所在的宏观状态识别为终末状态，将Ngn3low EP识别为初始状态。delta细胞由于数量稀少，其宏观状态的稳定性稍低，但仍可被识别。
- 图d展示了每个初始和终末宏观状态中最具代表性的30个细胞。
- 图e：展示了计算出的命运概率图谱。该图谱清晰地显示了不同祖细胞群对不同终末命运的偏好。例如，Ngn3high EP细胞群主要偏向于beta细胞命运，这与已知的生物学事实一致。
- 图f：结合伪时间，CellRank准确地重现了各个谱系关键调控因子（如Arx, Pdx1, Hhex）的动态表达模式。
实验二：揭示delta细胞分化的驱动基因
delta细胞在这个数据集中非常稀少，其分化机理尚不清楚，scVelo也未能有效捕捉其动力学。CellRank通过其全局分析能力，成功揭示了其分化路径。
- 图a, b, c：虽然RNA velocity的信号很弱，但CellRank的命运概率（fate probability）分析明确地指出了一条从Fev+祖细胞群中的一个特定亚群（Fev+ delta）分化到delta细胞的路径。
- 图d：通过关联基因表达与delta细胞的命运概率，CellRank识别出了一个基因表达的级联激活事件。在前50个最相关的基因中，不仅包含了已知的delta细胞标志物Hhex和Cd24a，还发现了一系列新的候选调控基因，如Hadh, Map2k4, Msi1等，为后续研究提供了新的靶点。
实验三：在细胞重编程中预测细胞命运
作者将CellRank应用于一个带有谱系追踪“金标准”的细胞重编程数据集，以验证其预测能力。在这个数据集中，只有约1%的细胞能成功重编程为iEPs，其余则进入“死胡同”（dead-end）状态。
- 图a, b：原始的t-SNE图和RNA velocity流线图未能清晰地揭示通往成功重编程的路径。
- 图c, d：CellRank成功地将细胞 coarse-grained 为5个宏观状态，并准确地将其中两个分别识别为“成功”状态和“死胡同”状态。
- 图e, f：这是最关键的验证。作者比较了CellRank预测的“成功”命运概率与通过CellTag谱系追踪得到的真实细胞命运。
- 图e直观地显示，CellRank的预测（上图）与真实的谱系标签（下图）高度吻合。
- 图f的ROC曲线分析定量地证明了这一点，在不同时间点，分类器的AUC值都非常高（例如，第21天时为0.94），表明CellRank的预测非常准确。
实验四：在肺再生中发现新的去分化轨迹
作者将CellRank应用于一个损伤后的小鼠肺再生数据集，以展示其在方向未知的复杂过程中的发现能力。
- 图a, b, c：CellRank分析了气道上皮细胞的动态，并识别出9个宏观状态。命运概率分析（图c）揭示，杯状细胞（goblet cells）有很高的概率转变为基底细胞（basal cells），这暗示了一条“去分化”（dedifferentiation）的路径。
- 图d, e：通过结合伪时间和命运概率，作者将这条去分化轨迹分为三个阶段：Stage 1（杯状细胞，高表达Bpifb1），Stage 2（中间态，同时表达Bpifb1和基底细胞早期标志物Krt5），以及Stage 3（基底细胞，高表达Krt5和Trp63）。
- 图f, g (实验验证)：为了验证这个新发现的中间态，作者在损伤后第10天的小鼠肺组织中进行了免疫荧光染色。
- 图f的图像清晰地显示，在损伤后的肺组织中，确实存在大量同时表达Bpifb1和Krt5但不表达Trp63的细胞，即Stage 2的中间态细胞。
- 图g的定量分析表明，这种中间态细胞在正常组织（PBS）或再生后期（bleo d22）中几乎不存在，但在损伤后第10天（bleo d10）显著富集。
这个发现表明，杯状细胞在损伤后可能通过去分化来补充干细胞库，CellRank成功地预测了这一先前未知的生物学过程。
与其他方法的性能比较
- 图a, b, c：在胰腺数据上，CellRank在自动识别初始/终末状态、预测命运偏好和揭示基因表达趋势方面，均优于Palantir, FateID, STEMNET和velocyto。
- 图d, e：在计算效率上，CellRank在处理10万个细胞时，仅需几分钟，并且内存消耗远低于Palantir, FateID和velocyto，显示了其卓越的可扩展性。
