0. TL; DR

单细胞测序技术的最新进展为测量单个细胞的基因表达谱和RNA velocity提供了前所未有的机会。然而，由于单细胞基因表达测量的高维、稀疏特性以及非线性的调控关系，对转录动态进行建模在计算上极具挑战性。

本文介绍了一种名为DeepVelo的方法，它是一个基于神经网络的常微分方程（ordinary differential equation）模型，能够通过描述单个细胞内连续时间的基因表达变化来模拟复杂的转录组动态。作者将DeepVelo应用于来自不同测序平台的公开数据集，以：（i）构建不同时间尺度上的转录组动态；（ii）衡量细胞状态的不稳定性；（iii）通过扰动分析识别发育过程中的驱动基因。

与当前最优方法的基准比较表明，DeepVelo能够学习到更准确的速度场表示。此外，作者的扰动研究揭示了单细胞动力系统可能表现出混沌特性。总之，DeepVelo为数据驱动地发现描绘单细胞转录组动态的微分方程提供了可能。

1. 背景介绍

单细胞RNA测序（scRNA-seq）的最新进展为我们在前所未有的分辨率下研究细胞发育开辟了新途径。scRNA-seq的一个主要目标是研究细胞的发育和分化，其中，解析细胞状态转变背后的调控级联对于理解胚胎发育、组织再生和肿瘤发生等关键生物过程至关重要。然而，大多数scRNA-seq实验只能捕捉到基因表达的“快照”，使得连续追踪同一个细胞的转录组谱变得困难。

现有的一些方法，如Monocle和Palantir，可以通过构建“伪时间”（pseudotime）来揭示大致的发育趋势。但这些方法并非为建模单个细胞内连续时间的转录组动态而设计。

最近，RNA velocity的提出通过利用已剪接与未剪接转录本的比例来估计瞬时的转录组变化，为动态建模提供了新思路。然而，RNA velocity只能预测未来数小时内的细胞状态。为了克服这种瞬时性，研究人员开发了基于线性ODE和稀疏近似的方法来预测更长时间尺度上的细胞状态演化。但转录调控是一个精确协调的生物过程，涉及众多调控因子之间的复杂非线性相互作用。因此，这些线性模型可能低估了调控的复杂性，无法捕捉基因表达动态的非线性特征。

受神经ODE和数据驱动的动力系统最新发展的启发，作者们提出了DeepVelo，一个基于神经网络的框架，用于构建scRNA-seq实验背后的动态模型并克服上述挑战。具体来说，DeepVelo训练一个变分自编码器（variational autoencoder, VAE）来从一个测量到的转录组谱中预测基因表达的变化率（即RNA velocity）。然后，作者们将这个VAE嵌入到一个ODE中，以模拟单个细胞内基因表达随时间的连续变化。

这种深度学习架构使得模型能够模拟复杂的非线性基因相互作用、降低基因表达动态的维度、学习对scRNA-seq中的速度场进行去噪。最重要的是，通过将不同发育阶段的细胞片段拼接起来，DeepVelo能够对更远的未来做出准确的细胞状态预测。与现有方法不同，DeepVelo旨在学习能够描述连续时间单细胞基因表达动态的神经微分方程。

2. DeepVelo 方法

DeepVelo框架包含两个主要步骤：推断基因相互作用和预测连续的细胞转变。

2.1 用VAE学习非线性基因相互作用

DeepVelo的核心思想是，将一个细胞的基因表达谱（$\vec{x}$）与其瞬时变化率（即RNA velocity，$\frac{\partial\vec{x}}{\partial t}$）通过一个神经网络连接起来。与假设线性基因相互作用（即$\frac{\partial\vec{x}}{\partial t} = A\vec{x}$）的现有方法不同，作者们训练一个变分自编码器（VAE）$f_A$来捕捉非线性的基因调控关系，并将基因表达状态映射到RNA velocity：

\[\frac{\partial\vec{x}}{\partial t} = f_A(\vec{x})\]

这个训练好的VAE $f_A$因此可以被看作是一个能够从任何（无论是观测到的还是未观测到的）细胞状态直接估计其瞬时基因表达变化率的函数。

VAE由一个编码器和一个解码器组成。编码器将输入（基因表达$\vec{x}$）映射到一个隐层，该隐层参数化一个高斯分布。从这个分布中采样一个隐变量$z_i$，然后解码器再将$z_i$映射回输出空间，以预测RNA velocity $\frac{\partial\vec{x}}{\partial t}$。训练过程通过最小化重构损失（MSE）和KL散度损失来优化。为了防止过拟合并鼓励稀疏表示，还加入了L1正则化。

2.2 用Neural ODE建模连续细胞转变

一旦有了能够预测任意状态下速度的函数$f_A$，作者们就可以通过将其嵌入到一个ODE中来模拟单个细胞随时间的连续基因表达变化。这在数学上是一个初值问题：给定$t=0$时的初始基因表达状态$\vec{x}_0$，可以通过积分来求解任意未来时间$t$的状态$\vec{x}(t)$。

\[\frac{\partial\vec{x}(t)}{\partial t} = f_A(\vec{x}(t)) \\ \vec{x}_{t+1} = \vec{x}_t + f_A(\vec{x}_t)\]

这个积分过程可以使用任何黑箱的ODE求解器（如欧拉法或更高阶的龙格-库塔法）来数值求解。通过迭代这个步骤，DeepVelo可以勾勒出单个细胞随时间的发育轨迹。

2.3 解决漂移效应（Drift Effects）

在控制理论中，仅使用前一个状态和速度来预测下一个状态会导致一种称为“航位推算”（dead reckoning）的现象，即误差会随着每一步的迭代而累积。为了减轻这种漂移效应，作者们采取了两种策略：

去噪VAE（Denoising VAE）：在训练VAE时，对输入的基因表达数据加入少量噪声。这可以增加输入空间的泛化能力，并提高对样本外细胞状态的外推性能。
在数据流形上寻找参考点：在ODE积分的过程中，每隔几个步骤，就将预测出的细胞状态投影回原始数据集，并找到其k-近邻。然后，将这些近邻的平均表达谱作为新的“参考细胞”，并从这个参考点继续进行积分。这确保了预测的轨迹能够紧密地贴合数据流形，进一步减少了自由度，并为动态系统构建了边界条件。

3. 实验分析

作者们将DeepVelo应用于多个来自不同组织、技术平台和发育时间尺度的公开scRNA-seq数据集，以展示其鲁棒性和普适性。

3.1 案例研究一：小鼠胰腺内分泌发生

作者们首先在一个E15.5的小鼠胰腺内分泌发生数据集上进行了概念验证。

图A, B（预测未来细胞状态）：
- 图A展示了数据集中不同细胞类型及其RNA velocity流。
- 图B展示了对一个“样本外”（out-of-sample）细胞的轨迹模拟。这个细胞（红色星号）从一个内分泌祖细胞状态开始，经过一个内分泌前体细胞状态，最终分化为beta细胞。模拟的轨迹（由紫到黄的星号）在UMAP空间中清晰地沿着已知的发育路径移动。值得注意的是，轨迹的步长反映了RNA velocity的大小：速度快的阶段步长大，速度慢的阶段步长小。VAE的隐空间嵌入（右下图）也正确地编码了发育的层级结构。
图D, E, F, G, H（用细胞临界指数CCI表征不稳定性）：作者们定义了一个“细胞临界指数”（Cell Criticality Index, CCI），它衡量了细胞在模拟的发育路径上基因表达分布的累积信息变化（KL散度）。CCI高的细胞状态不稳定，而CCI低的则处于稳定状态。
- 图D显示了CCI在UMAP空间的热图。可以看到，内分泌祖细胞的CCI最高，而导管细胞和已分化的内分泌细胞的CCI最低。
- 图E将CCI与scVelo推断的潜伏时间一起展示，CCI揭示了在潜伏时间轴上的一次剧烈的基因表达转变，这代表了细胞命运决定的“临界点”。
- 图F, G, H显示，与CCI正相关的基因（如Neurog3）是已知的发育驱动因子，而负相关的基因（如ATP1b1）则与维持细胞稳态相关。这证明了CCI作为一个衡量细胞状态不稳定性的指标是有效的。
图I, J, K（In silico扰动研究）：
- 图I显示，通过分析上游的差异表达基因，可以发现与alpha或beta细胞命运相关的候选基因。作者们假设，通过在内分泌前体细胞中“激活”这些候选基因，可以改变最终的细胞命运比例。
- 图J, K的结果表明，当激活与alpha细胞命运相关的“驱动基因”时，最终分化为alpha细胞的比例确实显著增加（从47%增加到61%）。这揭示了单细胞动力学系统可能具有类似于混沌系统的特性，即初始条件的微小扰动可能导致下游命运的巨大变化。

3.2 案例研究二：发育中的小鼠齿状回

作者们将DeepVelo应用于一个发育中的小鼠齿状回数据集，以测试其在不同组织和时间尺度上的适用性。

图A, B（轨迹模拟）：
- 图A展示了该数据集的UMAP图和velocity流。
- 图B模拟了一个样本外Nbl2细胞的轨迹，该轨迹正确地沿着已知的颗粒细胞发育路径前进。
图C, D, E（用CCI识别命运决定点）：结果再次验证了CCI的有效性。CCI在从Nbl1到Nbl2细胞的转变过程中达到峰值，表明这是细胞命运决定的一个关键不稳定点。
图J, K（In silico扰动）：
- 图J识别了与锥体神经元命运相关的上游差异表达基因。
- 图K的扰动模拟显示，激活这些基因可以显著增加最终分化为锥体神经元的比例。

3.3 案例研究三：发育中的小鼠新皮层

作者们进一步在一个跨越多个胚胎发育天数（E14-E17）的小鼠新皮层数据集上评估了DeepVelo。

图A, B（轨迹模拟）：模拟了一个中间祖细胞（IP）的轨迹，该轨迹正确地演变为一个迁移神经元（MN），最终分化为一个皮层投射神经元（CFN）。
图C, D, E（用CCI识别双相动态）：CCI揭示了一个双相的命运决定动态。AP（顶端祖细胞）的CCI较低，表明其主要功能是自我更新。而IP和MN的CCI较高，表明它们分别处于向神经元分化和向特定神经元亚型分化的关键不稳定点。
图F, G, H（In silico扰动）：DeepVelo成功地模拟了MN向ITN和CFN两种神经元亚型的分叉。扰动与ITN命运相关的上游差异表达基因，能够使最终分化为ITN的细胞比例平均增加10%。

3.4 应用：用逆行轨迹重构基因共表达网络

作者们提出了DeepVelo的一个新颖应用：通过模拟“逆行轨迹”（retrograde trajectories）来增强基因共表达网络的推断。

图A展示了其核心思想。从一个分化好的终末细胞状态（如beta细胞）开始，通过在ODE积分中减去速度向量（即时间倒流），可以模拟出导致该终末状态的“历史”轨迹。由于VAE具有去噪功能，这条逆行轨迹上的细胞相比于静态的原始细胞，其基因表达模式更平滑、噪声更低，因此能够揭示出更强的基因共表达关系（图B）。

作者们通过图C GO富集分析来定量评估这一效果。结果显示，在四个不同的数据集上，从逆行轨迹中发现的功能基因模块，其在细胞类型特异性GO词条上的富集程度，比从静态细胞中发现的模块高出至少两个数量级。

作者们将DeepVelo与线性和最先进的向量场学习方法SparseVFC进行了基准比较(图D)。在out-of-sample的速度预测任务中，DeepVelo在所有六个数据集上的预测误差都降低了至少50%，表明它学习到了更准确的速度场表示。