PRESCIENT用生成模型预测细胞干预下的命运轨迹.

paper：Generative modeling of single-cell time series with PRESCIENT enables prediction of cell trajectories with interventions

0. TL; DR

现有的单细胞计算方法大多无法模拟细胞在真实物理时间下的随机演化过程，也无法预测基因干预如何改变细胞的分化轨迹。本文介绍了一种名为PRESCIENT（Potential eneRgy undErlying Single Cell gradIENTs）的生成式建模框架。该框架能够从时间序列单细胞RNA测序（scRNA-seq）数据中学习一个潜在的分化“能量景观”。

通过在一个真实的谱系追踪数据集上进行验证，作者证明了在考虑细胞增殖的情况下，PRESCIENT能够准确预测造血过程中祖细胞的命运偏好，其性能超越了现有最优方法。更重要的是，PRESCIENT能够模拟基因被扰动后细胞的动态轨迹。作者在造血系统和胰腺β细胞分化模型中，成功地通过PRESCIENT的in silico（计算机模拟）干预实验，重现了已知细胞命运调控因子的预期效果。PRESCIENT支持在不同时间点、从不同起始细胞群对多个基因进行复杂的组合扰动，为大规模的in silico扰动实验提供了可能。

1. 背景介绍

理解细胞分化的“能量景观”对于揭示细胞如何走向不同的终末状态至关重要，并且有助于在体外精确地操控细胞命运。scRNA-seq技术通过在分化过程的多个时间点对单个细胞进行采样，为研究这一景观提供了可能。然而，这些研究只能提供分化过程的“快照”，无法直接观察到不同时间点细胞之间的谱系关系。

为了解决这个问题，领域内发展了多种计算方法，但它们各有局限：

伪时间（Pseudo-temporal）推断：这是最主流的方法，它将细胞沿着一个代表“分化时间”的一维坐标进行排序。然而，这个“时间”是任意的，并非真实的物理时间，因此无法模拟真实的动态过程。
细胞命运预测：一些方法专注于预测细胞的最终命运。例如，Waddington-OT将此问题构建为一个最优传输问题来预测细胞间的概率耦合；FateID则通过迭代构建分类器来预测命运。但这些方法本质上是对观测结果的“总结”，而不是对底层过程的“建模”。
过程建模方法：少数方法尝试将分化作为一个动态过程来建模。例如，Population Balance Analysis (PBA)使用偏微分方程描述分化，但受计算限制只能得到非参数解；pseudodynamics则只在一维细胞状态空间中模拟扩散过程。

这些现有方法普遍缺乏一个关键能力：预测基因干预（interventions）对细胞分化轨迹的影响。为此，作者引入了PRESCIENT，一个生成式建模框架。PRESCIENT将细胞分化建模为一个在高维基因表达空间中的扩散过程，该过程由一个潜在的能量函数驱动。通过学习这个能量景观，PRESCIENT不仅能描述细胞的自然分化过程，还能模拟在基因表达被计算扰动后，细胞的命运会如何改变，这是现有方法无法实现的。

2. PRESCIENT 方法

PRESCIENT的核心思想是将细胞分化过程建模为一个在基因表达“能量景观”上的扩散过程。这个景观由一个潜在的能量函数（potential function）所定义，而作者的目标就是从时间序列的scRNA-seq数据中学习这个函数。

2.1 核心方程：随机微分方程

作者将细胞状态的演化过程用一个随机微分方程（stochastic differential equation, SDE）来描述：

\[dX(t) = \mu(X(t))dt + \sqrt{2\sigma^2}dW(t)\]

其中$X(t)$表示一个细胞在时间$t$的$k$维状态（例如，基因表达的PCA投影）。$\mu(X(t))$是一个漂移项（drift term），代表在当前状态$X(t)$下作用于细胞的“力”。$W(t)$是单位布朗运动，代表了细胞演化过程中的随机性或“噪声”。$\sigma^2$：是扩散系数，控制噪声的强度。

一个关键的设定是，漂移项$\mu(x)$被定义为能量函数$\Psi(x)$的负梯度：

\[\mu(x) = -\nabla\Psi(x)\]

这个设定非常直观：细胞总是倾向于从高能量区域移动到低能量区域，就像一个球会从山顶滚到山谷一样。因此，能量函数$\Psi(x)$就对应了沃丁顿（Waddington）提出的表观遗传景观的高度。

通过对上述SDE进行一阶时间离散化，就可以模拟细胞的轨迹：

\[X(t + \Delta t) = X(t) + \mu(X(t))\Delta t + \sqrt{2\sigma^2\Delta t}Z(t)\]

其中$Z(t)$是独立的标准高斯噪声。

2.2 模型学习：基于Wasserstein距离的优化目标

PRESCIENT的目标是仅根据在不同时间点观测到的细胞群体“快照”来推断未知的能量函数$\Psi(x)$。为此，作者构建了一个优化问题，其目标函数如下：

\[\min_{\Psi \in K} \sum_{i=1}^{n} W_2^2(\hat{\rho}(t_i;x), \rho_{\Psi}(t_i;x)) + \tau \sum_{j=1}^{m_n} \frac{\Psi(x_j)}{\sigma^2}\]

这个目标函数包含两个部分：

Wasserstein距离：$W_2^2(\hat{\rho}, \rho_{\Psi})$是观测到的细胞群体分布$\hat{\rho}$与模型在能量函数$\Psi$下预测的分布$\rho_{\Psi}$之间的Wasserstein距离。最小化这个距离，就是让模型模拟出的细胞群体分布尽可能地接近真实观测。利用GeomLoss库高效地计算Wasserstein距离及其梯度，并使用GPU进行加速，使得模型可以处理大规模数据集。
熵正则化项：第二项是一个正则化项，用于惩罚过于复杂的能量函数，防止模型过拟合，其中$\tau$是正则化强度参数。

为了能够处理高维度的scRNA-seq数据并学习复杂的能量景观，作者使用一个softplus激活函数的双层神经网络来参数化能量函数$\Psi(x)$。

2.3 考虑细胞增殖

细胞在分化过程中会增殖和凋亡，这会影响后代细胞的数量。PRESCIENT通过在计算Wasserstein距离时为每个细胞赋予权重来考虑这一因素。具体来说，在从一个时间点映射到下一个时间点时，每个源细胞$x_i$的权重$\alpha_i$被设定为其预期的后代数量。这样，增殖能力更强的细胞在优化过程中会获得更大的“影响力”。

预期后代数量可以通过两种方式获得：

基于谱系追踪数据：如果实验中使用了谱系追踪技术，可以直接计算每个克隆在不同时间点的细胞数量来得到真实的增殖率。
基于基因表达估计：在没有谱系追踪数据时，可以使用KEGG等数据库中与细胞周期和凋亡相关的基因集的表达水平来估计增殖率。

3. 实验分析

作者通过一系列实验，全面验证了PRESCIENT的性能及其在in silico扰动模拟中的独特能力。

实验一：在造血谱系追踪数据集上的验证

作者使用了一个带有DNA条形码的小鼠造血系统谱系追踪数据集来验证PRESCIENT。

图a (Held-out时间点恢复)：模型在第2天和第6天的数据上训练，并被要求预测第4天（未参与训练）的细胞分布。结果显示，PRESCIENT模拟的第4天细胞群与真实观测的Wasserstein距离，显著低于其他基线方法（如WOT的线性插值），证明了其生成能力的准确性。
图b (增殖率估计)：作者发现，基于基因表达（KEGG基因集）估计的细胞增殖率与谱系追踪得到的真实增殖率之间存在显著的正相关（$r=0.207$），这为在没有谱系数据时使用估计值提供了依据。
图c, d (能量景观可视化)：与不考虑增殖的模型相比，考虑了细胞增殖的模型学习到的能量景观（potential）和漂移场（drift）在早期时间点附近有明显不同，这表明增殖信息对于正确刻画分化初始阶段至关重要。
图f, g (细胞命运预测)：作者评估了模型预测祖细胞最终分化为中性粒细胞（neutrophil）或单核细胞（monocyte）的偏好（clonal fate bias）的能力。
- 当不考虑细胞增殖时，PRESCIENT的性能（Pearson相关系数$r \approx 0.196$, AUROC $\approx 0.601$）与其他方法（如Waddington-OT和FateID）相当，表现平平。
- 当引入从谱系追踪数据中得到的真实增殖率后，PRESCIENT的性能大幅提升（$r \approx 0.347$, AUROC $\approx 0.692$）。
- 当使用基于基因表达估计的增殖率时，性能甚至略有超过（$r \approx 0.399$, AUROC $\approx 0.725$），这极大地增强了模型在缺乏谱系数据时的实用性。
图h (泛化能力)：作者发现，即使在训练时排除掉测试集中的细胞，模型的预测性能依然保持稳定，并且由于训练数据更多，模型能够对更大比例的测试细胞做出有效预测。这证明了PRESCIENT作为生成模型具有良好的泛化能力，能够预测未见过细胞的命运。

实验二：造血系统中转录因子的in silico扰动模拟

作者展示了PRESCIENT最强大的功能：模拟基因扰动。

图a, b：以一个例子直观展示了扰动的效果。在模拟中上调已知的促进中性粒细胞分化的转录因子（TF）后，最终的细胞群体中中性粒细胞的比例显著增加。
图c：对一系列已知的调控中性粒细胞或单核细胞分化的TF进行单基因扰动。结果显示，上调中性粒细胞相关TF会增加中性粒细胞比例，而上调单核细胞相关TF则会增加单核细胞比例，与已知的生物学功能完全一致。
图d, e (剂量效应)：扰动的影响是剂量依赖的。无论是上调还是下调，扰动的幅度越大，对最终细胞命运比例的影响也越大。
图f, g (组合效应)：相比于单基因扰动，同时扰动一组功能相关的TF（ensemble perturbations）会产生更强烈、更显著的细胞命运转变。
图h (对照实验)：作为阴性对照，同时扰动一组随机挑选的、与分化无关的非TF基因，最终的细胞命运比例没有发生显著变化。这证明了模型是稳健的，其预测的扰动效应是特异性的。

实验三：胰腺β细胞分化中的in silico扰动模拟

作者将PRESCIENT应用于另一个经典分化系统——体外诱导胰腺β细胞分化，进一步展示了其模拟复杂扰动的能力。

图c (内分泌/外分泌轴)：作者模拟了调控内分泌/外分泌谱系分化的关键TF。结果显示，在第0天过表达促进内分泌的TF（NEUROG3和NKX6.1）会显著增加最终内分泌细胞的比例；而过表达促进外分泌的TF（PTF1A和HES1）则效果相反。
图d (扰动时间的重要性)：内分泌/外分泌的分化发生在分化早期。PRESCIENT的模拟结果完美重现了这一“时间窗口”效应：在早期（如第0天）扰动NEUROG3/NKX6.1效果显著，而在晚期（如第4天）进行同样的扰动，其效果则大大减弱。

图a (大规模TF筛选)：作者进行了一次大规模的in silico筛选，模拟了超过200个TF的过表达效应。该筛选成功地“找回”了许多已知的调控α细胞和β细胞命运的关键TF（如ARX, PAX4, NKX2.2等），并为研究较少的肠嗜铬细胞（SC-EC）的命运调控提出了新的候选基因，展示了PRESCIENT作为假设生成工具的巨大潜力。
图b, c (扰动细胞状态的重要性)：除了时间，扰动的起始细胞状态也至关重要。
- 图b：对于一些在内分泌分化早期起作用的TF（如NKX2.2），在早期祖细胞中进行扰动能显著增加β细胞比例，但在更成熟的NEUROG3+细胞中扰动则效果甚微。
- 图b：而对于另一些在分化晚期依然活跃的TF（如PDX1, MAFA），无论是在早期祖细胞还是晚期NEUROG3+细胞中进行扰动，都能促进β细胞的生成。
- 图c：对于α细胞，几乎所有的有效扰动都发生在最早的祖细胞阶段，这表明α细胞的命运决定可能发生在分化过程的极早期。