从单细胞转录组数据中重构生长与动态轨迹.

0. TL; DR

在处理时间序列单细胞转录组（scRNA-seq）数据时，一个核心挑战是如何连接不同时间点的“快照”，因为测序过程会破坏细胞，导致细胞谱系和生长信息丢失。这篇论文介绍了一种名为TIGON（Trajectory Inference with Growth via Optimal transport and Neural network）的新算法。

TIGON基于动态非平衡最优传输（dynamic, unbalanced optimal transport）理论，能够同时推断细胞基因表达的velocity（速度，即状态转变）和growth（生长，即增殖与凋亡）。为了解决高维度的最优传输问题，作者引入了一种基于Wasserstein–Fisher–Rao (WFR)距离的无量纲公式，并利用深度学习方法（特别是神经常微分方程 neural ODEs）进行求解。

通过在模拟数据和三个真实的scRNA-seq数据集上的测试，TIGON不仅准确地重构了细胞轨迹和种群生长，还能推断未测量时间点的基因表达、时序性的基因调控网络（GRN）以及细胞间的通讯动态，展示了其强大的功能和广泛的应用前景。

1. 背景介绍

时间序列单细胞RNA测序（scRNA-seq）技术为我们观察细胞系统的动态过程提供了前所未有的机会。然而，这项技术本身具有破坏性，细胞在测序过程中会被杀死。因此，我们只能获得一系列离散时间点的细胞“快照”，而不同快照之间的细胞谱系关系或单个细胞的动态轨迹是缺失的。

为了解决这个问题，领域内已经发展出多种计算方法：

然而，现有方法大多忽略了一个关键的生物学过程：细胞种群的生长，即细胞的增殖和凋亡。细胞数量的变化会显著影响我们对细胞轨迹的推断。如果不考虑生长，可能会导致对动态过程的理解不完整甚至产生错误。一些前沿工作（如Waddington-OT和PRESCIENT）尝试通过KEGG或GO等数据库中的生长标志基因来近似估计生长率，但这种方法高度依赖于先验知识库的选择，可能引入偏见。

因此，当前迫切需要一种能够将基因表达速度和细胞种群生长同时纳入统一框架的计算工具。为此，作者提出了TIGON，一个基于深度学习的动态非平衡最优传输模型，旨在从时间序列scRNA-seq数据中同时、无偏地推断细胞velocity和growth。

2. TIGON 模型

TIGON的核心思想是将细胞群体在基因表达空间中的分布描述为一个随时间变化的密度函数 $ρ(x, t)$，其中 $x$ 是 $d$ 维基因表达状态， $t$ 是时间。这个密度函数的变化由一个包含平流项和生长项的偏微分方程（PDE）控制。

2.1 包含生长的连续性方程

细胞密度的动态演化遵循以下方程：

\[\partial_t\rho (x, t) + \nabla \cdot (\mathbf{v} (x, t)\rho (x, t)) = g (x, t)\rho (x, t)\]

其中：

这个方程直观地表明，一个区域内细胞密度的变化，来源于两个方面：细胞“流”进或“流”出该区域（由$\mathbf{v}$决定），以及该区域内细胞自身的增殖或死亡（由$g$决定）。

2.2 优化目标：Wasserstein–Fisher–Rao (WFR) 距离

为了求解上述方程中的未知函数 $\mathbf{v}$ 和 $g$，作者引入了动态非平衡最优传输理论，并使用Wasserstein–Fisher–Rao (WFR)距离作为优化目标。其目标函数（或称为cost）定义为：

\[W_{0,T} = \int_{0}^{T}\int_{\mathbb{R}^d} (|\mathbf{v} (x, t)|^2 + \alpha|g (x, t)|^2)\rho (x, t)dx dt\]

这个cost函数由两部分组成：

  1. $|\mathbf{v}(x, t)|^2$ 项：对应于Wasserstein度量，代表了细胞状态转变的“动能”消耗。
  2. $|g(x, t)|^2$ 项：对应于Fisher–Rao度量，代表了细胞种群生长或消亡的“能量”消耗。
  3. $\alpha$ 是一个超参数，用于平衡velocity和growth在总cost中的权重。

通过最小化这个WFR cost，算法可以找到在满足数据约束（即在给定的时间点上，计算出的密度与观测到的细胞分布相匹配）的同时，最为“经济”的velocity场$\mathbf{v}$和growth场$g$。

2.3 深度学习求解：无量纲公式与神经ODE

直接求解上述高维PDE和优化问题在计算上是极其困难的，存在“维度灾难”问题。TIGON采用了一种基于深度学习的无量纲（mesh-free）方法来解决这一难题。

作者使用两个全连接神经网络来近似未知的velocity和growth函数。通过数学推导，作者证明了可以将高维PDE问题转化为沿着单个细胞轨迹的常微分方程（ODE）系统。这使得计算可以从整个高维空间上的积分，简化为对一系列采样点进行轨迹追踪。

整个系统被构建为一个neural ODE模型。模型的训练过程如下：

2.4 下游分析

在训练完成后，TIGON可以提供丰富的下游分析：

3. 实验分析

作者通过一个模拟数据集和三个真实的scRNA-seq数据集，全面评估了TIGON的性能。

实验一：三基因模拟模型基准测试

作者构建了一个包含三个基因的随机模型，模拟了两种细胞动态：一部分细胞处于静态（quiescent），另一部分细胞从状态A转变到状态B，并在转变过程中增殖。

实验二：与谱系追踪实验结果的对齐

作者将TIGON应用于一个带有谱系追踪（lineage tracing）条形码的小鼠造血干细胞分化数据集。该数据集提供了细胞分化到中性粒细胞（Neu）和单核细胞（Mo）的真实谱系信息。

实验三：上皮-间质转化(EMT)过程的动态重构

作者分析了一个诱导A549癌细胞发生EMT的时间序列scRNA-seq数据集。

与其他方法的比较。

实验四：iPSC定向分化中的分叉识别

作者最后分析了一个人诱导多能干细胞（iPSC）向心肌细胞分化的qPCR数据集，该过程包含一个向中胚层（M）和内胚层（En）的分叉。