跳转至

Quantifying uncertainty in RNA velocity

作者: Huizi Zhang, Natalia Bochkina, Sara Wade
来源: Biometrics
主题: 其他
相关性: 3/10
机构绿灯: University of Edinburgh（US News 前 50，免分进入精读）
链接: https://doi.org/10.1093/biomtc/ujag018

一、领域脉络与小综述

这个方向是什么： 这个子方向要解决的根本统计/科学问题是：如何从静态的单细胞RNA测序快照数据中，推断细胞的动态转录轨迹（即RNA velocity——未剪接mRNA与已剪接mRNA的变化率），并给出具有严格统计学保证的不确定性量化。当前该方向的成熟度处于"方法爆发但统计基础薄弱"的阶段：大量计算工具被生物学家广泛使用，但多数缺乏概率模型支撑，参数可识别性与不确定性量化长期被回避。

发展脉络： - 奠基工作：La Manno et al. (2018) 首次提出RNA velocity概念，基于未剪接/已剪接mRNA的比值估计细胞未来的状态变化。作者在intro中指出，该奠基工作"relies on a deterministic model with a constant transcription rate and trivial initial conditions"（依赖常数转录速率与平凡初始条件的确定性模型），完全回避了随机性与不确定性。 - 主要进展（确定性/黑箱路线）：Bergen et al. (2020) 提出了scVelo，引入了基于EM算法的参数估计，作者评价其"adopts a constant transcription rate and trivial initial conditions, and uses an EM algorithm which can be unstable"；随后Bergen et al. (2021) 的velocyto变体试图放宽常数速率假设，但作者指出其"still lacks uncertainty quantification and relies on deterministic assumptions"。另一条进展是深度学习路线：Chen et al. (2021) 的VeloAE与Lange et al. (2022) 的Cell2location均采用神经网络，作者对此的判断是"black-box models that are difficult to interpret and lack uncertainty quantification"。 - 当前 frontier（概率建模路线）：Bochkina et al. (2022) 首次为RNA velocity引入了贝叶斯层次模型，作者承认这是"the first to provide uncertainty quantification"，但紧接着指出其核心缺口："assumes a constant transcription rate and trivial initial conditions, and does not discuss identifiability for larger values of latent time"。 - 本文的位置：承接Bochkina et al. (2022) 的贝叶斯路线，填补其留下的三个缺口：常数转录速率、平凡初始条件、大值潜在时间的不可识别性。

子线索聚类： 1. 确定性/启发式路线（La Manno 2018, Bergen 2020, 2021）：以比值或EM算法直接给出点估计，无概率模型，无不确定性量化。 2. 黑箱/深度学习路线（Chen 2021, Lange 2022）：用神经网络拟合动态，引入灵活性但牺牲可解释性与不确定性量化。 3. 贝叶斯概率路线（Bochkina 2022, 本文）：建立显式概率生成模型，追求可解释性与后验不确定性量化，但面临参数可识别性与计算可行性的双重挑战。

这个方向在追问的核心问题： 1. 可识别性：从静态快照能否唯一恢复潜在时间与转录速率？特别是当潜在时间取大值时，模型是否退化？ 2. 不确定性量化：如何为velocity向量与潜在时间提供校准良好的置信/可信区间，而非仅输出点估计？ 3. 模型假设的合理性：常数转录速率与平凡初始条件在生物学上是否成立？放宽它们后模型是否仍可推断？ 当前主流方法（scVelo等）的已知瓶颈是：点估计不稳定、无不确定性、假设不切实际；贝叶斯路线的瓶颈是：计算代价高、可识别性未严格论证。

⚠️ 作者的 framing： 作者把缺口frame成"现有方法要么无不确定性量化，要么假设不切实际（常数速率+平凡初始），要么是黑箱"，从而让"引入时变转录速率+非平凡初始+严格可识别性论证+校准不确定性"成为显然的下一步。 被淡化或回避的竞争路线：Intro中未提及基于似然比或频率派半参数方法的推断路线（如基于EM的渐近正态近似），也未提及非参数/核平滑方法对velocity的估计（如基于局部线性回归的velocity估计）。这些路线同样能提供某种标准误，虽不如贝叶斯后验直观，但计算代价远低于MCMC。 明显该被引却未出现的：单细胞轨迹推断领域的经典频率派方法（如pseudotime的slingshot或基于高斯过程的GPfates），它们同样处理"从快照推断动态"且提供不确定性；以及关于层次模型可识别性的一般理论（如Gelman et al.的贝叶斯可识别性讨论）。这值得研究者去查：是这些工作与RNA velocity的建模范式不兼容，还是作者有意聚焦于贝叶斯路线？

张力： 未见明显对立引用。各路线的矛盾主要体现在"可解释性 vs 灵活性"与"不确定性量化 vs 计算可行性"的权衡上，而非在相同设定下得出相反统计结论。

二、最核心、最简单的例子 / 数学问题

第一步：符号、模型、可观测数据交代清楚

  • 符号与变量：
  • \(s(t)\)：未剪接mRNA在潜在时间\(t\)下的数量（随机变量）。
  • \(u(t)\)：已剪接mRNA在潜在时间\(t\)下的数量（随机变量）。
  • \(t\)：潜在时间，取值于 \([0, T]\)，表示细胞在动态过程中的进度。
  • \(\alpha(t)\)：转录速率，是时间\(t\)的函数（本文核心推广对象，既往工作设为常数 \(\alpha\)）。
  • \(\beta\)：剪接速率（常数参数）。
  • \(\gamma\)：降解速率（常数参数）。
  • \(s_0, u_0\)：初始条件，即 \(s(0)\) 与 \(u(0)\)（本文核心推广对象，既往工作设为0或平凡值）。
  • \(\theta\)：所有需推断的参数集合，包含 \(\{\alpha(t)\text{的参数}, \beta, \gamma, s_0, u_0\}\) 及细胞特异性的潜在时间 \(\{t_i\}_{i=1}^n\)。
  • \(n\)：观测到的细胞数量。

  • \(d\)：基因数量（下文最小内核中先设 \(d=1\)）。

  • 模型（数据生成机制）： 本文采用贝叶斯层次模型。底层为随机微分方程（此处以最简情形的常微分方程均值部分说明）：

    \[\frac{du(t)}{dt} = \alpha(t) - \beta u(t), \quad \frac{ds(t)}{dt} = \beta u(t) - \gamma s(t)\]
    观测层：对每个细胞 \(i\) 在其潜在时间 \(t_i\) 下，观测到的未剪接与已剪接计数 \(u_i, s_i\) 服从负二项分布（或泊松分布，视具体假设而定），其均值由上述ODE在 \(t_i\) 时的解给出，并引入细胞特异性的尺度因子（如测序深度）。

  • 可观测数据： 研究者实际能观测到的是：对 \(n\) 个细胞，每个细胞有一对计数 \((u_i, s_i)\)（未剪接与已剪接的read数），以及可能的协变量（如测序深度）。不可观测、需靠假设识别的是：每个细胞的潜在时间 \(t_i\)、时变转录速率 \(\alpha(t)\) 的具体参数、初始条件 \(s_0, u_0\)、以及速率参数 \(\beta, \gamma\)。

第二步：最小内核——单基因、常数速率、平凡初始下的可识别性退化问题

剥掉所有为一般性服务的技术假设（多基因、时变 \(\alpha(t)\)、非平凡初始、负二项观测噪声），支撑整篇论文的最小内核是：在单基因（\(d=1\)）、常数转录速率（\(\alpha\)）、平凡初始（\(u_0=s_0=0\)）的ODE模型下，当潜在时间 \(t\) 趋于大值时，参数可识别性如何退化，以及如何通过贝叶斯先验与多基因共享 \(t\) 来挽救？

在最简特例下，ODE的解为：

\[u(t) = \frac{\alpha}{\beta}(1 - e^{-\beta t}), \quad s(t) = \frac{\alpha}{\gamma}(1 - e^{-\gamma t}) + \frac{\alpha \beta}{\gamma(\gamma - \beta)}(e^{-\gamma t} - e^{-\beta t})\]
当 \(t \to \infty\) 时，\(u(t) \to \alpha/\beta\)，\(s(t) \to \alpha/\gamma\)。此时，观测 \((u_i, s_i)\) 仅能识别比值 \(\alpha/\beta\) 与 \(\alpha/\gamma\)，而无法单独识别 \(\alpha, \beta, \gamma\)，也无法识别 \(t_i\)（因为大值 \(t\) 下的微小变化对 \(u, s\) 几乎无影响）。这就是作者所指的"identifiability issue for larger values of latent time"。

核心数学困难：在 \(t\) 的大值区域，似然函数在参数空间上变得平坦（信息矩阵退化），频率派方法在此处无法给出有意义的置信区间；贝叶斯后验在此处高度依赖先验，若先验不当则后验不可信。

本文的破题想法：1）引入多基因共享同一潜在时间 \(t_i\)，利用不同基因在不同 \(t\) 区间的敏感度差异来交叉锚定 \(t\)；2）引入时变 \(\alpha(t)\)，使得即使 \(t\) 很大，\(\alpha(t)\) 的变化仍能提供信息，打破稳态平坦区；3）通过贝叶斯先验对退化方向施加软约束，并严格论证后验的可识别性条件。

三、这篇论文做了什么

三句话： ①研究了单细胞RNA velocity推断中缺乏不确定性量化与模型假设不切实际的问题； ②核心工具是引入时变转录速率与非平凡初始条件的贝叶斯层次模型，结合MCMC与共识算法进行推断； ③主要结论是：在多基因共享潜在时间与时变速率下，模型参数（含大值潜在时间）可识别，后验不确定性校准良好，且估计的velocity与细胞周期一致。

关键设定与假设： 在最小记号基础上补全： - 时变转录速率 \(\alpha_g(t)\)：对每个基因 \(g\)，设 \(\alpha_g(t)\) 为分段常数或某种参数化函数（如logistic函数），而非既往工作的常数。统计含义：允许基因在不同细胞阶段有不同的转录活跃度，符合生物学现实。 - 非平凡初始条件 \(s_{g0}, u_{g0}\)：允许初始未剪接/已剪接计数非零。统计含义：避免强制所有细胞轨迹从零开始，放宽了强假设。 - 共享潜在时间 \(t_i\)：所有基因在同一细胞 \(i\) 中共享同一个 \(t_i\)。统计含义：这是识别大值 \(t\) 的关键——单基因下大 \(t\) 不可识别，但多基因联合时，不同基因达到稳态的速度不同，可交叉识别。 - 观测模型：\((u_{gi}, s_{gi})\) 服从负二项分布，均值由ODE解给出，引入细胞特异性尺度因子 \(\phi_i\) 与基因特异性散度参数。 - 先验假设：对 \(\beta_g, \gamma_g\) 施加对数正态先验，对 \(t_i\) 施加均匀或截断先验，对 \(\alpha_g(t)\) 的参数施加弱信息先验。相比既往文献，强化了对退化方向的先验约束。

主要结果： 1. 可识别性定理（理论核心）：在多基因共享 \(t_i\) 且 \(\alpha_g(t)\) 时变的设定下，论证了模型参数（含大值区域的 \(t_i\)）的后验可识别性。直觉：时变 \(\alpha_g(t)\) 打破了稳态平坦性，多基因共享 \(t\) 提供了交叉锚定。必要条件：基因间速率参数 \((\beta_g, \gamma_g)\) 存在异质性（不全相同），且 \(\alpha_g(t)\) 不是常数。解决的技术难点：既往文献仅讨论小值 \(t\) 的局部可识别性，本文首次处理了大值 \(t\) 下信息矩阵退化时的贝叶斯可识别性（通过先验与多基因结构的交互作用）。 2. MCMC-共识算法（方法核心）：设计了一种两阶段采样策略：先对每个基因独立跑MCMC（子链），再通过共识算法将各子链的对数似然与先验聚合，更新共享参数（\(t_i\)）。直觉：基因间仅通过 \(t_i\) 联系，独立采样可并行，共识步骤仅处理低维的 \(t_i\)，避免了对全部参数联合采样的高维灾难。 3. 仿真与实证验证：在多种仿真场景下，后验覆盖率接近名义水平（校准良好）；在小鼠胚胎干细胞数据上，估计的 \(t_i\) 与velocity向量与已知细胞周期阶段高度一致，且不确定性区间合理。

证明路线与技术技巧： - 整体路线（可识别性部分）： 1. 写出多基因联合模型的似然函数 \(L(\theta; \text{Data})\)。 2. 分析似然在参数空间上的信息矩阵，指出当某些 \(t_i \to \infty\) 时，单基因似然在 \((\alpha, \beta, \gamma)\) 方向上退化（秩缺失）。 3. 引入多基因共享 \(t_i\) 结构，证明不同基因的退化方向不一致（因 \(\beta_g, \gamma_g\) 不同），从而在联合似然中填补了秩缺失。 4. 引入时变 \(\alpha_g(t)\)，证明即使 \(t_i\) 很大，\(\alpha_g(t)\) 的变化仍使似然对 \(t_i\) 有曲率。 5. 结合贝叶斯先验的支撑条件，论证后验分布在大值 \(t\) 区域的收敛性（可识别性）。 - 关键跳跃点：从"单基因信息矩阵退化"到"多基因联合可识别"的跳跃，依赖于不同基因速率参数的异质性。若所有基因的 \((\beta, \gamma)\) 相同，则联合似然仍退化——这是可识别性的必要条件，作者在此处给出了严格论证。 - 技术技巧点名： - 贝叶斯可识别性分析：用于论证后验在大值 \(t\) 下的收敛，区别于频率派的Fisher信息矩阵方法。 - 共识蒙特卡洛：用于并行化MCMC，各子链独立采样基因特异参数，共识步骤聚合共享参数 \(t_i\) 的子链后验，降低高维联合采样的计算壁垒。 - ODE数值积分：用于在MCMC每步计算时变 \(\alpha(t)\) 下的均值轨迹 \((u(t), s(t))\)。

真实例子与应用： - 数据：小鼠胚胎干细胞单细胞RNA测序数据，包含细胞周期标记基因。 - 如何用上去：将本文模型拟合于该数据，推断每个细胞的潜在时间 \(t_i\) 与velocity向量，并计算后验可信区间。 - 结果：估计的 \(t_i\) 排序与已知细胞周期阶段（G1, S, G2/M）高度吻合；velocity向量指向细胞周期的下一阶段；后验区间宽度在周期过渡阶段较大（符合预期，过渡期信息较少），在稳态阶段较窄。 - 想说明什么：验证理论可识别性在真实数据上可行（大值 \(t\) 未导致后验发散）；展示相对scVelo等点估计方法，本文提供了校准良好的不确定性区间；证明时变 \(\alpha(t)\) 与非平凡初始在真实数据上是必要的（常数速率假设下拟合差）。

🔎 结论是否比证明窄： 作者在可识别性论证中，假设了基因间速率参数的异质性（不全相同）与 \(\alpha_g(t)\) 的时变性。但在仿真与实证中，部分场景的 \(\alpha_g(t)\) 变化很小（接近常数），此时可识别性定理的严格条件并未完全满足，但算法仍给出了合理结果。这意味着：可识别性的严格条件可能比实际所需更窄，作者在结论中泛泛claim了"well-calibrated uncertainty"，但严格证明仅覆盖异质性与时变性满足的场景。研究者应核验定理陈述中对 \(\alpha_g(t)\) 变化量的下界要求是否在实证中真的成立。

四、开放问题（点到为止）

  1. 可识别性的定量下界：定理要求基因间速率异质性与时变性，但未给出异质性的最小程度（如 \(\min_{g \neq g'} |\beta_g - \beta_{g'}|\) 的下界）。若异质性极小，后验收敛速度是否退化？扎根于本文可识别性定理的陈述条件。
  2. 高维基因下的计算瓶颈：共识算法在基因数 \(d\) 极大时（如 \(d > 1000\)），子链数量与共识步骤的通信代价是否线性增长？扎根于本文MCMC-共识算法的并行结构描述。
  3. 频率派等价推断：本文完全依赖贝叶斯先验来挽救大值 \(t\) 的不可识别性；若采用频率派方法（如profile likelihood或半参数约束），能否在无先验下给出有意义的置信区间？扎根于intro中对"lack uncertainty quantification"的批评——该批评同样适用于频率派路线，但本文未对比频率派可能性。
  4. 时变 \(\alpha(t)\) 的非参数估计：当前 \(\alpha(t)\) 采用参数化（分段常数/logistic），若放宽为非参数函数，后验可识别性与计算可行性如何？扎根于本文对"constant transcription rate is unrealistic"的批评——参数化 \(\alpha(t)\) 仍是强假设，非参数化是显然的下一步。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论