使用SDEvelo对具有细胞特异性潜伏期的转录动力学进行多元随机建模.
0. TL; DR
RNA velocity的出现为单细胞RNA测序（scRNA-seq）研究带来了范式转变，使得重构和预测细胞分化及状态转变的有向轨迹成为可能。然而，大多数现有的动态建模方法使用常微分方程（ODE）对单个基因进行建模，而没有采用多变量方法。这种建模策略未能充分捕捉由跨多个基因的、细胞特异性的潜伏时间（latent time）所控制的、内在随机的转录动态，可能导致错误的结果。
为了解决这些问题，作者们提出了SDEvelo，一个通过多变量随机微分方程（multivariate stochastic differential equations, SDE）来建模未剪接和已剪接RNA动态，从而推断RNA velocity的生成式方法。SDEvelo的独特之处在于，它在估计跨基因的细胞特异性潜伏时间的同时，显式地建模了转录动态中固有的不确定性。
通过在模拟数据以及四个真实的scRNA-seq和空间转录组学数据集上的验证，作者们证明了SDEvelo能够模拟成熟状态细胞的随机动态模式，并能准确检测癌变过程。此外，估计出的基因共享的潜伏时间有助于许多下游的生物学发现分析。作者们还证明了SDEvelo在计算上是可扩展的，并且同时适用于scRNA-seq和基于测序的空间转录组学数据。
1. 背景介绍
单细胞RNA测序（scRNA-seq）技术彻底改变了我们对分子和细胞过程的理解，为在单细胞分辨率下探索基因表达谱、识别稀有细胞类型、追踪细胞谱系等提供了前所未有的机会。轨迹推断（trajectory inference）方法旨在利用在发育过程的不同阶段对细胞群体进行的scRNA-seq测量，来重构单个细胞随时间的发育路径。然而，一个主要的挑战是这些技术的破坏性，它们只能捕捉到转录组的静态“快照”。因此，传统的轨迹推断方法通常需要足够的细胞样本覆盖所有中间状态，并且需要先验知识来确定分化方向。
为了从描述性的轨迹模型转向预测性的动态模型，RNA velocity应运而生。它通过对mRNA的转录、剪接和降解动力学进行建模，来恢复动态信息。一个正的velocity表示某个基因的剪接后转录本正在被上调，负的则表示下调。然而，现有的RNA velocity方法存在一些不容忽视的局限性：
- 忽略基因间的协同性：大多数方法对每个基因独立建模，忽略了在同一个细胞中，多个基因的动态是由一个共同的细胞特异性潜伏时间（cell-specific latent time）所协同调控的。这可能导致在稳态细胞中出现错误的命运方向预测。
- 忽略内在随机性：大多数方法基于确定性的ODE模型，虽然考虑了测量误差等外在噪声，但没有考虑转录过程本身固有的随机性和不确定性。
- 对不同技术的适用性：将现有的RNA velocity方法应用于结构复杂、测序深度较低的空间转录组学数据仍然是一个挑战。
为了解决这些挑战，作者们提出了SDEvelo，一种基于SDE的RNA velocity推断方法。SDEvelo旨在通过一个统一的框架，实现以下目标：
- 协同建模：通过一个多变量SDE框架，协同估计跨多个基因的RNA velocity，并由一个共享的细胞特异性潜伏时间所控制。
- 捕捉随机性：明确地对转录动力学中的内在随机过程进行建模。
- 推断潜伏时间：推断出能够揭示细胞分化进程的细胞特异性潜伏时间。
- 广泛适用性：能够同时应用于scRNA-seq和空间转录组学数据。
SDEvelo的核心思想是，用一个多变量随机微分方程（multivariate SDE）系统来取代传统RNA velocity方法中对单个基因使用的确定性ODE，并引入一个非线性的、可微的转录速率函数和一个细胞特异性的潜伏时间。整个模型通过一个生成式对抗学习框架进行训练。
对于一个细胞中的$p$个基因，作者们用$p$维向量$U(t)$和$S(t)$分别表示它们在细胞特异性潜伏时间$t$时的未剪接和已剪接mRNA丰度。它们的动态由以下多变量SDE系统描述：
\[dU(t) = (\alpha(t) - \beta \odot U(t))dt + \sigma_1 dB(t) \\ dS(t) = (\beta \odot U(t) - \gamma \odot S(t))dt + \sigma_2 dB(t)\]其中$\alpha(t)$, $\beta$, $\gamma$分别是$p$维的转录速率、剪接速率和降解速率向量。$\odot$表示元素级别的乘法。$t$是一个对该细胞内所有基因都共享的潜伏时间。$B(t)$是维纳过程（Wiener process），代表了内在的随机波动。$\sigma_1$和$\sigma_2$是控制噪声水平的参数。
与传统的ODE模型相比，这个SDE模型有三个关键创新：
- 多变量框架：它将所有$p$个基因的动态耦合在一起，通过一个共享的潜伏时间$t$来协同建模，而不是独立处理每个基因。
- 随机性建模：通过引入随机项$\sigma dB(t)$，模型能够显式地捕捉转录过程中的内在随机性。
- 非线性转录速率：作者们不再假设转录速率$\alpha(t)$是分段常数，而是使用一个非线性且可微的sigmoid函数来建模，这更符合基因调控的生物学现实：
其中，对于第$i$个基因，$a_i$是“转换时间”点，$c_i$是最大转录速率，$b$是控制转换陡峭度的超参数。
由于上述多变量SDE系统没有解析解，传统的基于EM算法的参数估计方法不再适用。因此，SDEvelo采用了一个生成式对抗学习的框架来训练模型。
生成器（Generator）即作者们的SDE模型。在训练的每一步，它使用当前的参数$\Theta = {a, c, \beta, \gamma, \sigma_1, \sigma_2}$，并通过欧拉-马力山方法（Euler-Maruyama method）对SDE进行数值积分，来生成一批模拟的$(U, S)$数据。判别器（Discriminator）是通过一个可微的损失函数——最大均值差异（Maximum Mean Discrepancy, MMD）——来充当“判别器”的角色。MMD通过一个核函数（如高斯核）来衡量真实数据分布与生成数据分布之间的距离。
模型的目标是最小化MMD损失$L_{MMD^2}(\Theta)$。通过随机梯度下降（stochastic gradient descent）来更新参数$\Theta$，使得生成的数据分布越来越接近真实的数据分布。
一旦模型训练完成，获得了最优的参数$\Theta$，SDEvelo就可以为每个真实的细胞推断其潜伏时间。这是通过一个类似最优传输的思想来实现的：
1.使用训练好的SDE模型生成大量的模拟轨迹$G(\Theta)$，每个模拟点$g_i(\Theta)$都带有一个已知的潜伏时间$t_i$。
2.寻找一个从模拟数据点到真实观测数据点$R$的一一映射（bijection）$\phi^*$，使得总的“运输成本”（欧几里得距离的平方和）最小,这个问题可以通过Sinkhorn-Knopp算法高效求解。
\[\phi^* = \arg\min_{\phi: G(\Theta) \to R} \sum_{i=1}^{N_g} |g_i(\Theta) - \phi(g_i(\Theta))|^2\]3.一旦找到了最优的映射 \(\phi^*\)，对于每个真实细胞$r_i$，就可以找到它在模拟数据中对应的点 \(g_i(\Theta) = (\phi^*)^{-1}(r_i)\)，并将其潜伏时间$t_i$赋给该真实细胞。
通过这种方式，SDEvelo为每个细胞估计了一个跨多个基因共享的、连续的潜伏时间。
3. 实验分析
作者们通过一系列模拟实验和在四个真实数据集上的应用，全面地评估了SDEvelo的性能。
3.1 模拟数据验证
作者们在两种设置下生成了模拟数据：基于确定性ODE（外加噪声）和基于SDE。在两种设置下，SDEvelo都表现出色。
- 图e的KDE散点图（左）显示，SDEvelo估计的潜伏时间与真实值高度相关（Pearson相关系数接近1），而scVelo (dyn)的估计则存在明显偏差。SDEvelo的速度流线图也准确地重现了真实的动态方向。
- 图f的定量比较显示，在估计稳态比率（steady-state ratio，上图）和潜伏时间（下图）方面，SDEvelo的误差最小，相关性最高，并且在不同基因数量下都保持了稳健的性能。
3.2 对成熟状态细胞群体的负向控制
一个理想的RNA velocity方法应该能够在没有动态变化的成熟细胞群体中检测到随机、无方向的速度场。
作者们将SDEvelo和其他方法应用于一个PBMC（外周血单个核细胞）数据集，这些细胞大多处于成熟的稳态。
- 图b显示，SDEvelo正确地检测到了随机、无方向的速度场。而scVelo (dyn), VeloVI, UniTVelo等方法却错误地推断出了强烈的、但显然是人为的定向模式。
- 图c：在一个模拟的无动态数据集上，SDEvelo同样成功地识别出随机的速度场，而其他方法则推断出错误的方向。
- 图d的潜伏时间热图显示，SDEvelo将大部分细胞的潜伏时间估计为接近终点（接近1），这与它们是成熟细胞的事实相符。而scVelo (dyn)等方法则错误地将大部分细胞置于早期阶段。
- 图e, f进一步分析了几个被scVelo错误识别为“驱动基因”的基因。SDEvelo的分析表明，这些基因的动态实际上是随机的，从而避免了错误的生物学推断。
3.3 SDEvelo正确推断肝细胞癌变的空间动态
作者们将SDEvelo应用于一个包含肿瘤组织（HCC1-2）和癌旁组织（HCC3-4）的肝细胞癌（HCC）空间转录组学数据集。
- 图a：第二行和第四行显示了SDEvelo和scVelo (dyn)估计的潜伏时间热图。SDEvelo的潜伏时间模式与病理学家手动标注的基质区（stroma，早期）和肿瘤/正常上皮区（TNE，晚期）高度一致。而scVelo等方法则无法得到清晰或正确的模式。
- 图b的速度流线图显示，SDEvelo正确地推断出从基质区到TNE区的癌变方向。
- 图c, d的定量评估（CBDir, AUC, Rank Correlation）表明，SDEvelo在速度方向的正确性和潜伏时间的一致性上，均显著优于所有其他比较方法。
- 图e, f (下游分析)：通过关联基因表达与SDEvelo推断的潜伏时间，作者们识别出了一系列与癌变过程强相关的驱动基因，其中许多是载脂蛋白（APO）基因。通过构建这些基因的蛋白-蛋白相互作用（PPI）网络，并与TCGA数据库进行比较，作者们发现这些基因在肝癌（LIHC）中显著富集，证明了SDEvelo推断结果的生物学意义。
3.4 SDEvelo预测小鼠胚胎重编程的细胞命运
作者们将SDEvelo应用于一个从小鼠胚胎成纤维细胞（MEF）到诱导内胚层祖细胞（iEPs）的重编程数据集，该数据集带有真实的实验时间记录。
- 图a, b显示，SDEvelo推断的速度流和潜伏时间热图与真实的重编程天数高度一致，而scVelo (dyn)等方法则出现了方向错误。
- 图c, d的定量比较证实，SDEvelo估计的潜伏时间与真实实验天数的皮尔逊相关系数最高，AUC也最优，甚至优于经典的伪时间方法Monocle3。
- 图e (命运轨迹识别)：SDEvelo成功地区分了两条不同的轨迹：一条通往成功重编程的状态（高表达Apoa1），另一条则通往“死胡同”状态（高表达Col1a2）。通过与CellRank结合，模型进一步识别出了这两个终末状态，并发现了驱动这两条不同路径的候选基因。
3.5 SDEvelo准确识别小鼠红细胞生成中的转变
作者们在一个小鼠红细胞生成数据上验证了SDEvelo。
- 图a, b：SDEvelo的速度流线图（图a）和潜伏时间热图（图b）都清晰、平滑地重现了从造血祖细胞（BP1/2）到成熟红细胞（Ery1-3）的分化过程，而scVelo的结果则显得混乱且方向错误。
- 图c：展示了SDEvelo作为一个生成模型的能力，它能为存在转录爆发或方向逆转的复杂基因动态生成灵活且符合生物学现实的轨迹。
- 图e-h (下游分析)：通过对与潜伏时间相关的基因进行富集分析（图f）和细胞间通讯分析（图g, h），SDEvelo揭示了在红细胞分化过程中起关键作用的生物学通路（如细胞增殖、红细胞分化）和配体-受体相互作用（如Tgfb1-Itgb1）。
