使用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方法存在一些不容忽视的局限性：

  1. 忽略基因间的协同性：大多数方法对每个基因独立建模，忽略了在同一个细胞中，多个基因的动态是由一个共同的细胞特异性潜伏时间（cell-specific latent time）所协同调控的。这可能导致在稳态细胞中出现错误的命运方向预测。
  2. 忽略内在随机性：大多数方法基于确定性的ODE模型，虽然考虑了测量误差等外在噪声，但没有考虑转录过程本身固有的随机性和不确定性。
  3. 对不同技术的适用性：将现有的RNA velocity方法应用于结构复杂、测序深度较低的空间转录组学数据仍然是一个挑战。

为了解决这些挑战，作者们提出了SDEvelo，一种基于SDE的RNA velocity推断方法。SDEvelo旨在通过一个统一的框架，实现以下目标：

2. SDEvelo 框架

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模型有三个关键创新：

  1. 多变量框架：它将所有$p$个基因的动态耦合在一起，通过一个共享的潜伏时间$t$来协同建模，而不是独立处理每个基因。
  2. 随机性建模：通过引入随机项$\sigma dB(t)$，模型能够显式地捕捉转录过程中的内在随机性。
  3. 非线性转录速率：作者们不再假设转录速率$\alpha(t)$是分段常数，而是使用一个非线性且可微的sigmoid函数来建模，这更符合基因调控的生物学现实：
\[\alpha(t)_i = \frac{c_i}{1 + \exp\{b(t-a_i)\}}\]

其中，对于第$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都表现出色。

3.2 对成熟状态细胞群体的负向控制

一个理想的RNA velocity方法应该能够在没有动态变化的成熟细胞群体中检测到随机、无方向的速度场。

作者们将SDEvelo和其他方法应用于一个PBMC（外周血单个核细胞）数据集，这些细胞大多处于成熟的稳态。

3.3 SDEvelo正确推断肝细胞癌变的空间动态

作者们将SDEvelo应用于一个包含肿瘤组织（HCC1-2）和癌旁组织（HCC3-4）的肝细胞癌（HCC）空间转录组学数据集。

3.4 SDEvelo预测小鼠胚胎重编程的细胞命运

作者们将SDEvelo应用于一个从小鼠胚胎成纤维细胞（MEF）到诱导内胚层祖细胞（iEPs）的重编程数据集，该数据集带有真实的实验时间记录。

3.5 SDEvelo准确识别小鼠红细胞生成中的转变

作者们在一个小鼠红细胞生成数据上验证了SDEvelo。