CardamomOT: a mechanistic optimal transport-based framework for gene regulatory network inference, trajectory reconstruction and generative modeling

本文提出了 CardamomOT，一种基于机理最优传输的新框架，通过整合精确时间标签和先验知识，在单细胞转录组时间序列分析中实现了基因调控网络推断、未观测蛋白轨迹重建及扰动响应生成的统一与高精度建模。

要点

这是一篇关于生物信息学新工具 CardamomOT 的论文介绍。为了让你轻松理解，我们可以把细胞内的基因调控想象成一场宏大的交响乐演出，而 CardamomOT 就是一位超级侦探兼指挥家。

🎻 核心故事：看不见的指挥家与断断续续的乐谱

1. 面临的难题：只有乐谱，没有指挥
想象一下，细胞里的基因（DNA）就像乐谱，蛋白质（Proteins）是指挥家，而 mRNA 是乐手们演奏的声音。

• 现实困境：科学家通常只能听到乐手们（mRNA）发出的声音（通过测序技术），而且只能听到瞬间的快照（比如每隔一小时拍一张照片）。
• 致命伤：
  • 我们看不见指挥家（蛋白质），但指挥家才是决定音乐走向的关键。
  • 乐手们的演奏非常随机且嘈杂（转录爆发），光听声音很难推断出指挥家到底在打什么拍子。
  • 以前的方法（如旧版 CARDAMOM）就像是在猜：“既然乐手 A 在 1 点响了，乐手 B 在 2 点响了，那它们之间肯定有关系。”但这往往猜不准，因为忽略了指挥家（蛋白质）的缓慢变化和真实的时间节奏。

2. 新工具登场：CardamomOT
这篇论文提出的 CardamomOT 就像是一个拥有“透视眼”和“时间机器”的超级侦探。它不仅能听到声音，还能反推出看不见的指挥家（蛋白质）在每一刻的动作，并还原出整场演出的完整乐谱（基因调控网络）。

🕵️‍♂️ 它是如何工作的？（三个魔法步骤）

CardamomOT 的工作流程可以比作修复一部被撕碎且模糊的侦探电影：

• 第一步：整理线索（预处理）
  它先仔细分析每一张“快照”（单细胞数据），把混乱的乐手声音分类。它把细胞分成不同的“状态组”（比如：准备睡觉的、正在跳舞的、正在逃跑的）。这就像先把观众按情绪分组。

• 第二步：侦探推理与时间旅行（核心循环）
  这是最精彩的部分，它在一个循环中不断自我修正：

  • A. 预测指挥家的动作（轨迹重建）：
    它利用“最优传输”（Optimal Transport，一种数学方法，可以想象成最省力的搬运路径）来猜测：如果细胞 A 在 1 点是这个状态，细胞 B 在 2 点是那个状态，那么中间的指挥家（蛋白质） 是怎么移动的？它不再假设指挥家是乱跑的，而是根据生物学规律（比如蛋白质降解速度）来推算最合理的移动路线。
  • B. 修正乐谱（基因网络推断）：
    一旦猜出了指挥家的动作，它就能反过来推断：是谁在指挥谁？是“基因 X"在指挥“基因 Y"吗？它根据刚才猜出的蛋白质轨迹，重新绘制基因之间的控制关系图（GRN）。
  • C. 循环验证：
    它用新画出的关系图，再去修正对指挥家动作的猜测；再用新的猜测去修正关系图。就像侦探不断比对线索，直到逻辑完美闭环。

• 第三步：生成“平行宇宙”（生成模型）
  一旦它完全搞懂了这套规则，它就不再只是观察者，而变成了创造者。它可以模拟出“如果没发生这件事会怎样”或者“如果强行让某个基因过表达会怎样”的平行宇宙。

🌟 它的超能力（主要贡献）

1. 透视隐形指挥家：
   以前的方法只能看 mRNA（声音），CardamomOT 能重建蛋白质（指挥家）的轨迹。就像它不仅能听到音乐，还能在脑海里画出指挥家挥舞指挥棒的完整视频。这让推断出的基因关系更准确。

2. 利用真实时间：
   以前的方法只能知道“先发生 A，后发生 B"，但不知道具体过了多久。CardamomOT 能利用精确的时间标签，结合已知的生物学知识（比如蛋白质降解有多快），让推理更精准，不需要那么多人为设定的“猜测参数”。

3. 预测未来（生成式模型）：
   这是最酷的一点。训练好之后，它就像一个数字孪生体。

   • 例子：科学家想测试“如果我把基因 X 关掉，细胞会变成什么样？”
   • 以前必须真的去实验室做实验（杀细胞、重做）。
   • 现在，CardamomOT 可以在电脑里模拟出关掉基因 X 后的结果。论文中，它成功预测了某些基因过表达会提高干细胞重编程的效率，这与后来的实验结果完全一致！

🎯 总结：为什么这很重要？

如果把细胞分化（比如干细胞变成皮肤细胞）比作从山顶滑向山谷：

• 旧方法：只能看到山顶和山谷的几张照片，猜中间的路径，容易猜错，因为不知道重力（蛋白质动力学）是怎么作用的。
• CardamomOT：不仅猜出了路径，还还原了重力场，甚至能告诉你：“如果我把这块石头（某个基因）移走，滑行的路线会怎么变？”

一句话总结：
CardamomOT 是一个基于物理规律的 AI 侦探，它通过结合数学（最优传输）和生物学原理，从混乱的单细胞数据中，还原了看不见的蛋白质舞蹈，画出了精准的基因控制地图，并能预测未来的细胞命运。这大大加速了我们对生命过程的理解和新药研发的过程。

技术摘要

这篇论文介绍了一种名为 CardamomOT 的新框架，旨在通过单细胞 RNA 测序（scRNA-seq）时间序列数据，联合推断基因调控网络（GRN）、重建未观测的蛋白质轨迹，并构建生成式模型。该方法基于**机制性最优传输（Mechanistic Optimal Transport, OT）**理论，解决了现有方法在因果推断、动力学建模和生成能力方面的局限性。

以下是该论文的详细技术总结：

1. 研究背景与问题 (Problem)

• 核心挑战：从实验数据中推断控制细胞分化等过程的基因调控网络（GRN）面临巨大困难。主要障碍在于无法在单细胞水平直接测量蛋白质动力学，导致难以建立调节因子活性与靶标响应之间的因果关系。
• 现有方法的局限：
  • 数据特性：标准 scRNA-seq 是破坏性的，只能提供不同时间点的独立快照，缺乏连续追踪。
  • 转录组噪声：mRNA 合成具有高度随机性（转录爆发），且 mRNA 与蛋白质水平之间存在复杂的非线性关系，仅基于 mRNA 的统计模型往往无法准确反映真实的生物学动力学。
  • 现有 OT 方法的不足：基于最优传输（OT）的轨迹推断方法通常假设细胞遵循布朗运动或线性过程（如 Ornstein-Uhlenbeck 过程），忽略了基因调控的机制性约束（Mechanistic constraints）。此外，之前的方法（如作者提出的 CARDAMOM）依赖准稳态假设（quasi-stationary approximation），无法利用精确的时间标签，且需要大量超参数。

2. 方法论 (Methodology)

CardamomOT 在一个统一的机制性最优传输框架下，通过迭代过程联合优化 GRN 和蛋白质轨迹。

2.1 生物学模型基础

• 混合两态模型：采用基于生物物理原理的基因表达模型。基因在“非激活”和“激活”状态间随机切换，激活时产生 mRNA 爆发，进而翻译为蛋白质。
• 动力学方程：
  • mRNA：受爆发频率 $k_{on}$ 和降解率 $d_0$ 控制，条件分布近似为负二项分布（NB）。
  • 蛋白质：在 mRNA 快速降解的极限下，蛋白质动力学由确定性微分方程描述：
    $$\frac{dP}{dt} = v_1(P) = d_1 \left( \frac{k_{on}^\theta(P)}{k_1} - P(t) \right)$$
    其中 $k_{on}^\theta(P)$ 是依赖于蛋白质浓度 $P$ 和 GRN 参数 $\theta$ 的爆发频率函数（通常使用 Sigmoid 函数建模）。
• 隐马尔可夫视角：将细胞状态离散化为“势阱”（basins），细胞在势阱内遵循确定性动力学，在势阱间发生随机跃迁。

2.2 CardamomOT 核心算法流程

算法采用类似期望最大化（EM）的迭代循环，包含三个主要步骤：

1. 预处理与初始化：

   • 对每个基因和时间点拟合负二项混合模型，推断机制参数（爆发频率、大小等）。
   • 将细胞初始化为离散的“势阱”标签（basin labels）。

2. 迭代推断（核心步骤）：

   • 步骤 1：基于机制性 OT 的蛋白质轨迹重建
     • 给定当前的 GRN ($\theta$) 和势阱标签，计算相邻时间点细胞间的传输成本。
     • 创新点：成本函数不是简单的欧氏距离，而是基于机制模型。通过积分确定性 ODE（允许在区间内发生一次模式切换），计算从 $t_j$ 的细胞 $c$ 到 $t_{j+1}$ 的候选细胞 $c'$ 的蛋白质轨迹成本。
     • 求解熵正则化的最优传输问题，得到细胞间的耦合（coupling），从而推断出每个细胞的完整蛋白质轨迹 $\{P_c(t)\}$。
   • 步骤 2：GRN 更新
     • 利用重建的蛋白质轨迹，通过回归最小化预测的爆发率与观测到的势阱模式之间的损失，更新 GRN 参数 $\theta$。
     • 引入弹性网络惩罚（Elastic Net），可结合先验知识。
   • 步骤 3：势阱标签细化
     • 根据更新后的 GRN，重新分配细胞的势阱标签，平衡 NB 似然度与机制模型的一致性。

3. 后处理：

   • 利用推断出的 GRN 和轨迹，通过 NeuralODE 方法重新校准动力学参数（如蛋白质降解率 $d_1$），确保模拟数据能复现观测到的随机性和时间尺度。

3. 关键贡献 (Key Contributions)

1. 联合推断框架：首次在一个统一框架中显式地联合推断 GRN 结构和未观测的蛋白质轨迹，放弃了 CARDAMOM 中的准稳态假设。
2. 机制性最优传输：将 OT 的成本函数从几何距离替换为基于生物物理机制（GRN 驱动的 ODE）的代价，使轨迹推断更符合生物学原理，而非简单的布朗运动假设。
3. 利用精确时间标签与先验知识：能够直接利用实验的时间标签和文献中的蛋白质降解率先验，显著提高了鲁棒性，并减少了超参数依赖。
4. 生成式建模能力：校准后的模型是一个生成式模型，不仅能重建数据，还能在 in silico 环境中模拟未见过的遗传扰动（如基因敲除或过表达）对细胞命运的影响。

4. 实验结果 (Results)

研究在合成数据集和三个真实的实验数据集上进行了验证：

• 数据集：
  1. 小鼠胚胎干细胞（mESC）分化（Semrau et al.）。
  2. 交感肾上腺分化（Kameneva et al.，伪时间排序）。
  3. 成纤维细胞重编程为 iPSC（Schiebinger et al.）。
• GRN 推断性能：
  • 在多个基准网络拓扑上，CardamomOT 的 AUPR（精确率 - 召回率曲线下面积）显著优于现有方法（如 CARDAMOM, Reference Fitting, GENIE3, SINCERITIES）。
  • 能够准确推断相互作用的符号（激活/抑制）和方向。
  • 对蛋白质降解率的输入误差具有鲁棒性（即使误差达 $\pm 50\%$ 仍表现优异）。
• 轨迹与速度场重建：
  • 重建的蛋白质速度场与真实机制模型的高度一致（余弦相似度 > 0.8），远超基于 mRNA 的线性模型。
  • 成功恢复了隐藏的蛋白质轨迹和细胞状态转换。
• 生成与扰动预测：
  • 模型能高质量地再生实验数据的分布、相关性结构和细胞类型比例。
  • 扰动实验：成功预测了 Obox6 和 Zfp42 过表达会提高 iPSC 重编程效率（与实验验证一致），且无需预先知道这些因子的功能。还预测了 Dnmt3a 过表达会抑制多能性，揭示了其表观遗传调控作用。

5. 意义与影响 (Significance)

• 理论突破：将最优传输理论与机制性生物物理模型深度融合，解决了传统 OT 方法缺乏生物学机理约束的问题。
• 工具价值：提供了一个从原始 scRNA-seq 计数矩阵到可解释的 GRN、蛋白质动力学及生成式模拟的完整流水线。
• 应用前景：
  • 数字孪生：校准后的模型可作为生物系统的“数字孪生”，用于预测药物或基因编辑的效果，减少湿实验成本。
  • 因果推断：通过显式建模蛋白质动力学，增强了从相关性数据中推断因果调控关系的能力。
  • 可扩展性：虽然计算成本随基因数量平方增长，但在 100 个基因规模下仍具有可行性，且未来可通过引入先验知识（如 ChIP-seq 数据）来约束参数空间，解决更大规模网络的识别性问题。

综上所述，CardamomOT 代表了单细胞时间序列分析领域的一个重要进展，它通过引入机制性约束和最优传输理论，显著提升了基因调控网络推断的准确性、可解释性和预测能力。