通过动态模型推广RNA速度到瞬时细胞状态.
0. TL; DR
RNA velocity框架将其速度计算为观测到的剪接与未剪接mRNA比例与推断出的稳态（steady state）比例的偏差。如果其核心假设——即所有基因共享一个共同的剪接率，并且能够观测到包含稳态水平的完整剪接动态——被违反，那么速度估计就会出现错误。
为了解决这些限制，作者们开发了scVelo，通过求解完整的转录动力学来解决剪接动力学问题。scVelo能够推断出基因特异性的转录、剪接和降解速率，并恢复出底层细胞过程的潜伏时间（latent time）。这个潜伏时间代表了细胞的“内部时钟”，并且仅仅基于其转录动态。此外，scVelo还能识别调控变化的区域，如细胞命运决定的阶段，并在此过程中系统性地检测出潜在的驱动基因。
通过在海马齿状回神经发生和胰腺内分泌细胞生成数据集上的演示，作者们证明了scVelo能够以前所未有的分辨率解析异质亚群的动力学。
1. 背景介绍
单细胞转录组学使得在单细胞分辨率下无偏见地研究细胞分化和谱系选择等生物过程成为可能。由此产生的计算问题被称为轨迹推断（trajectory inference）。然而，单细胞RNA-seq的破坏性是一个核心挑战，它只能揭示细胞状态的静态“快照”。
RNA velocity的概念通过利用新转录的、未剪接的pre-mRNA和成熟的、已剪接的mRNA可以被区分这一事实，为恢复有向的动态信息提供了可能。通过一个简单的反应模型，可以推断出mRNA丰度的变化，即RNA velocity。然而，最初的模型（在此称为“稳态模型”）在推断稳态比例时做出了两个基本假设：
- 在基因层面，能够捕捉到完整的剪接动态，包括转录的诱导、抑制和稳态mRNA水平。
- 在细胞层面，所有基因共享一个共同的剪接率。
这些假设在实际中常常被违反，特别是在一个群体包含具有不同动力学的多个异质亚群时。为了解决这些限制，作者们开发了scVelo，一个基于可能性的动态模型，它通过求解完整的基因特异性转录动力学，从而将RNA velocity的估计推广到瞬时系统和具有异质亚群动力学的系统。
scVelo能够推断基因特异性的转录、剪接和降解速率，以及一个所有基因共享的潜伏时间。这个潜伏时间代表了细胞的内部时钟，准确地描述了细胞在底层生物过程中的位置。与基于相似性的伪时间方法不同，这个潜伏时间仅基于转录动态，并考虑了运动的速度和方向。
2. scVelo 方法
scVelo的核心在于从一个基于常微分方程（ODE）的剪接动力学模型中，通过一个期望最大化（Expectation-Maximization, EM）框架来推断模型的参数。
2.1 剪接动力学的动态模型
与原始的稳态模型一样，scVelo的动态模型也是基于一个描述转录、剪接和降解过程的反应动力学方程组。对于每个基因，其未剪接mRNA丰度$u(t)$和已剪接mRNA丰度$s(t)$的变化由以下ODE描述：
\[\frac{du(t)}{dt} = \alpha_k(t) - \beta u(t) \\ \frac{ds(t)}{dt} = \beta u(t) - \gamma s(t)\]其中$\alpha_k(t)$是转录速率，它依赖于转录状态$k$（例如，“开”或“关”）。$\beta$是剪接速率。$\gamma$是降解速率。RNA velocity被定义为已剪接mRNA丰度的变化率，即$v(t) = \frac{ds(t)}{dt}$。
与稳态模型的关键不同之处在于，scVelo不假设已经观测到了稳态。相反，它直接求解上述ODE系统。这个系统的解析解为：
\[\begin{aligned} u(t) &= u_0 e^{-\beta\tau} + \frac{\alpha(k)}{\beta}(1 - e^{-\beta\tau}) \\ s(t) &= s_0 e^{-\gamma\tau} + \frac{\alpha(k)}{\gamma}(1 - e^{-\gamma\tau}) + \frac{\alpha(k) - \beta u_0}{\gamma - \beta}(e^{-\gamma\tau} - e^{-\beta\tau}) \end{aligned}\]其中，$\tau = t - t_0^{(k)}$是相对于状态切换时间的相对时间，$u_0, s_0$是初始丰度。
2.2 通过期望最大化（EM）推断参数
剪接动力学由两组参数控制：
- 反应速率：$\theta = (\alpha(k), \beta, \gamma)$。
- 细胞特异性潜伏变量：每个细胞$i$的离散转录状态$k_i$和连续时间$t_i$。
这两组参数相互依赖，因此作者们采用期望最大化（EM）算法来迭代地推断它们。
- 初始化：使用稳态模型的估计结果来初始化参数。
- E-步（Expectation）：在给定当前反应速率$\theta$的估计下，为每个细胞的观测值$(u_i, s_i)$分配一个最优的潜伏时间$t_i$和转录状态$k_i$。
- 时间分配：通过最小化观测点到当前模型预测的相位轨迹（phase trajectory）$(û(t), ŝ(t))$的距离来为每个细胞分配一个潜伏时间$t_i$。
- 状态分配：根据观测点与轨迹上不同部分（诱导、抑制、稳态）的距离，为每个细胞分配一个最可能的转录状态$k_i$。
- M-步（Maximization）：在给定潜伏变量（时间$t_i$和状态$k_i$）的分配下，通过最大化数据的总似然（likelihood）来更新反应速率参数$\theta$。
通过迭代这两个步骤，模型最终会收敛到最优的参数估计。
2.3 推断基因共享的潜伏时间
上述模型是针对单个基因的。为了将在不同基因上学习到的动力学联系起来，scVelo引入了一个“通用”的、所有基因共享的潜伏时间概念，它代表了细胞的“内部时钟”。这个通用潜伏时间是通过耦合所有基因的基因特异性潜伏时间来计算的。它不仅使得不同基因的反应速率具有可比性，还提供了一个比伪时间更可靠的、基于真实转录动态的时间度量。
2.4 扩展到随机模型
为了更好地捕捉基因表达的随机性，作者们还将稳态模型扩展到了一个随机框架。通过将转录、剪接和降解视为概率性事件，可以推导出一个包含一阶矩（均值）和二阶矩（协方差）的动力学方程组。这个随机模型不仅利用了剪接与未剪接mRNA的丰度平衡，还利用了它们的协变关系，从而能够更好地捕捉动态方向。
3. 实验分析
作者们在一系列模拟和真实的单细胞数据集上展示了scVelo的动态模型的优越性。
3.1 解析齿状回发育中的异质群体动力学
作者们将scVelo应用于一个包含两个时间点（P12和P35）的小鼠海马齿状回（dentate gyrus, DG）神经发生的数据集。
- 图a（速度流与一致性）：动态模型和稳态模型都捕捉到了从神经母细胞到颗粒细胞的主要分化谱系。然而，在亚群中，只有动态模型正确地识别出OPC（少突胶质前体细胞）分化为OL（少突胶质细胞），并将CR（Cajal-Retzius）细胞识别为终末状态（速度为零）。稳态模型则错误地为CR细胞分配了高速度，并指错了OPC的分化方向。动态模型推断的速度场在邻近细胞间具有更高的一致性（consistency scores更高）。
- 图b（基因层面的解释）：作者们发现稳态模型的错误可追溯到单个基因的错误状态判定。例如，对于基因Fam155a，稳态模型错误地将作为“离群点”的CR细胞分配了高速度。而动态模型则正确地将这些细胞识别为处于稳态。对于驱动神经母细胞分化的关键基因Tmsb10，动态模型显示出更平滑、更一致的速度模式。
scVelo为每个基因计算一个似然值，该值衡量了其表达动态被模型解释得有多好。通过对基因按似然值进行排序，可以系统性地识别出那些表现出显著动态行为的基因，这些基因是驱动生物过程的候选“驱动基因”。
图c高似然的基因（上排）显示出清晰的剪接动态模式，而低似然的基因（下排）的表达则主要由噪声主导。许多被识别出的高似然基因（如Grin2b, Map1b, Dlg2）都是已知的在神经发生中起关键作用的基因。通过只在部分细胞（如某个谱系或某个阶段）中计算似然，还可以识别特定转变阶段的驱动基因。
3.2 描绘内分泌细胞生成中的循环、定型和命运转变
作者们将scVelo应用于小鼠胰腺内分泌细胞生成的数据集。
- 图a, b, c（与稳态模型的比较）：
- 图a显示，动态模型准确地描绘了内分泌祖细胞（EP）的细胞周期、谱系定型（lineage commitment）、细胞周期退出和最终的内分泌分化过程。
- 相比之下，稳态模型（图b）未能捕捉到细胞周期，并且在晚期内分泌细胞中产生了错误的、指向去分化方向的“回流”。
- 图c的基因层面分析表明，这种回流源于稳态模型对α细胞的错误状态判定（例如，对于基因Cpe）。
- 图d, e, f, g（通过潜伏时间解析动态）：
- 图d显示，scVelo推断的潜伏时间能够正确反映不同内分泌细胞产生的真实时间顺序（α细胞较早，β细胞较晚），而基于相似性的伪时间则无法区分。
- 图e通过统计基因转录状态转换点（switching points）的数量，准确地识别出了谱系定型和分叉点的区域。
- 图f展示了沿潜伏时间轴的基因表达“瀑布流”，清晰地揭示了转录激活的时序。
- 图g展示了几个被识别为高似然驱动基因的表达动态，它们分别在细胞周期退出、谱系定型等关键节点被激活。
