Inference after latent variable estimation for single-cell RNA sequencing data¶
作者: Anna Neufeld, Lucy L Gao, Joshua Popp, Alexis Battle, Daniela Witten
来源: Biostatistics
主题: 数理统计 / 假设检验
相关性: 6/10
链接: 期刊页 · arXiv
一、领域脉络与小综述¶
这个方向是什么 这个子方向要解决的根本统计问题是“数据复用导致的假阳性膨胀”（即选择性推断 / 后选择推断问题）。在单细胞 RNA 测序数据分析中，研究者通常先用全量数据估计表征细胞状态的隐变量（如细胞类型聚类或连续分化轨迹伪时间），随后再用同一批数据逐基因检验该隐变量与基因表达的关联。由于隐变量是从同一批数据中“挖出来”的，第二步的 \(p\)-值不再满足名义 Type 1 误差控制。当前该方向的成熟度处于方法框架刚建立、特定应用场景正在被攻坚的阶段：选择性推断在回归与聚类中已有较成熟的条件推断框架，但在单细胞计数数据这种具有强 Poisson 结构、且隐变量估计步骤极其多样的场景下，如何给出一个通用、不依赖特定算法的推断框架，仍是未决问题。
发展脉络 - 奠基工作（选择性推断的条件化框架）：Lee et al. (2013) 与 Fithian et al. (2014) 提出了 selective Type 1 error 的概念与条件推断框架，通过在给定选择事件的条件下计算 \(p\)-值，恢复了频率学派的误差控制。这为所有后选择推断问题定下了基本原则。 - 主要进展（聚类与降维中的后选择推断）：Gao et al. (2020) 与 Chen & Witten (2022) 将选择性推断应用到层次聚类与 K-means，通过精确计算选择事件的几何约束，给出了有限样本下精确的 \(p\)-值。Chung & Storey (2013, 2020) 则从另一条路线出发，提出 jackstraw 方法，通过置换残差来检验主成分或聚类成员的显著性，试图绕过精确的条件化计算。 - 当前 frontier（单细胞数据中的数据复用问题）：Lähnemann et al. (2020) 明确将“数据复用导致的假阳性”列为单细胞数据科学的十一大挑战之一；Deconinck et al. (2021) 指出轨迹推断中的循环性导致“人为偏低的 \(p\)-值”；Zhang et al. (2019) 针对单细胞聚类后的差异表达提出了一个特定的修正框架。 - 本文的位置：本文避开需要针对每个算法定制选择事件的条件推断路线，也避开需要特定模型假设的 jackstraw 路线，转而利用单细胞数据的 Poisson 分布结构，提出 count splitting——在泊松假设下，将每个基因的计数随机拆分为两份独立子集，一份用于隐变量估计，一份用于推断，从而将后选择推断问题转化为经典的样本拆分问题。
子线索聚类 1. 条件选择性推断路线：以 Lee et al. (2013), Fithian et al. (2014), Gao et al. (2020), Chen & Witten (2022) 为代表。这一簇通过精确刻画“选择事件”（如 lasso 选中某变量、K-means 划分某簇），在给定选择事件的条件下推导检验统计量的分布。瓶颈：必须针对每个具体的隐变量估计算法写出选择事件的解析表达，对于 Monocle3、Seurat 等极其复杂的黑箱算法几乎不可行。 2. 残差置换 / Jackstraw 路线：以 Chung & Storey (2013, 2020) 为代表。通过在原数据上估计隐变量后，对残差进行置换以构造零分布。瓶颈：置换破坏了数据的原始依赖结构，且在隐变量估计步骤非线性、非参数时，残差的定义与置换的合理性难以保证。 3. 数据拆分 / Binomial thinning 路线：以 Gerard (2020) 为代表，提出用二项拆分从真实数据生成合成数据以评估方法。本文的 count splitting 是这一思路在推断问题上的直接演进：Gerard (2020) 用 binomial thinning 做方法评估，本文用 Poisson 下的 count splitting 做严格推断。
这个方向在追问的核心问题 1. 如何在不依赖特定隐变量估计算法解析形式的前提下，实现后选择推断的 Type 1 误差控制？（当前主流条件推断路线要求算法可解析，这是核心瓶颈。） 2. 单细胞计数数据的分布结构能否被利用来构造通用的数据拆分机制？（Poisson 假设是否足够？偏离时后果多严重？） 3. 在保证 Type 1 误差的前提下，如何尽可能保留检验功效？（样本拆分天然损失功效，如何最小化这种损失？）
⚠️ 作者的 framing（这是作者的说法） - 作者将缺口 frame 为：现有选择性推断路线（条件化）要求针对每个算法定制，而单细胞领域的隐变量估计步骤（如 Monocle3 的降维+构图+排序）过于复杂，无法写出选择事件；同时，传统的样本拆分（将细胞随机分成两半）在隐变量是细胞级别时不可行（拆分后两半细胞的隐变量不再一致）。因此，利用 Poisson 结构的 count splitting 是“显然的下一步”。 - 被淡化或回避的竞争路线：作者对 jackstraw 路线仅在引用中提及，未深入比较其与 count splitting 在功效与适用范围上的差异；对 Zhang et al. (2019) 的单细胞专用修正框架也未做正面对比。 - 明显该被引却未出现的：在讨论 Poisson 假设的合理性时，作者引用了 Sarkar & Stephens (2021) 支持“Poisson 测量模型足够”的观点，但未引用任何对 UMI 数据负二项或零膨胀模型提供强证据的工作（如 Van den Berge et al. (2019) 的 tradeSeq 明确使用负二项分布）。这是一个值得研究者去查的缺口：Poisson 假设的稳健性到底有多大？
张力 未见明显对立引用。但存在一个隐含张力：Sarkar & Stephens (2021) 主张 UMI 数据的 Poisson 测量模型足够且应优先使用，而大量差异表达分析实践（如 tradeSeq）使用负二项分布以捕捉过度散布。本文的 count splitting 严格依赖 Poisson 假设，这两条路线的假设基础存在张力。
二、最核心、最简单的例子 / 数学问题¶
第一步：符号、模型、可观测数据交代清楚
- \(n\)：细胞数（样本量）。
- \(p\)：基因数（特征数）。
- \(X_{ij}\)：可观测的随机变量，表示第 \(i\) 个细胞、第 \(j\) 个基因的 UMI 计数（read count）。这是一个 \(n \times p\) 的计数矩阵。
- \(\gamma_i\)：第 \(i\) 个细胞的 size factor（已知或已准确估计的测序深度偏移参数），反映细胞间总计数量的差异。
- \(\Lambda_{ij}\)：不可观测的潜在量（latent true expression rate），表示第 \(i\) 个细胞第 \(j\) 个基因的真实表达速率。它包含了我们想要推断的隐变量效应。
- \(Z_i\)：不可观测的潜在量（latent variable），表示第 \(i\) 个细胞的隐状态（如离散的细胞类型标签，或连续的伪时间值）。这是第一步要估计的对象。
- 模型（数据生成机制）：
\[X_{ij} \sim \text{Poisson}(\gamma_i \Lambda_{ij})\]且所有 \(X_{ij}\) 相互独立。这里 \(\gamma_i\) 是已知偏移，\(\Lambda_{ij}\) 是真实信号。在第二步检验中，我们假设 \(\Lambda_{ij}\) 依赖于 \(Z_i\)，例如 \(\Lambda_{ij} = \exp(\beta_{0j} + \beta_{1j} Z_i)\)。
- 可观测数据：研究者实际能观测到的是整个计数矩阵 \(X \in \mathbb{N}^{n \times p}\) 以及 size factor \(\gamma_1, \dots, \gamma_n\)。不可观测的是 \(Z_i\)（隐变量）和 \(\Lambda_{ij}\)（真实表达速率），只能靠假设与算法去估计。
第二步：最小内核
整篇论文的证明与方法本质上是泊松分布可加性这一性质的推广与包装。最简特例是：单个基因（\(p=1\)），单个细胞（\(n=1\)），无 size factor（\(\gamma=1\)）。
在这个特例下： - 可观测数据：\(X \sim \text{Poisson}(\Lambda)\)。 - 我们想做的：先从 \(X\) 估计隐变量 \(Z\)（在此特例中，\(Z\) 的估计可能就是 \(X\) 本身或某个函数 \(\hat{Z} = f(X)\)），然后检验 \(X\) 是否与 \(\hat{Z}\) 关联。 - 核心数学困难：\(\hat{Z}\) 是 \(X\) 的函数，用 \(X\) 去检验 \(\hat{Z}\) 必然导致循环依赖，Type 1 误差膨胀。传统的样本拆分是把 \(n\) 个细胞分成两半，但这里 \(n=1\)，无法拆分细胞。 - 本文的破局点（最小内核）：利用 Poisson 分布的性质。若 \(X \sim \text{Poisson}(\Lambda)\)，我们可以生成一个随机拆分比例 \(\delta \in (0,1)\)，令：
三、这篇论文做了什么¶
三句话 ① 研究了单细胞 RNA-seq 数据中，先估计隐变量（细胞类型/伪时间）再检验基因与隐变量关联时，数据复用导致 Type 1 误差膨胀的问题。 ② 核心工具是 count splitting：在 Poisson 计数假设下，将每个细胞的每个基因的计数按比例 \(\delta\) 拆分为两份独立子集，一份用于隐变量估计，另一份用于推断。 ③ 主要结论是：在 Poisson 假设下，count splitting 对任意隐变量估计算法和任意推断方法均保证 Type 1 误差控制，且在模拟与真实数据中展示了可接受的检验功效。
关键设定与假设 在第二节最小记号的基础上补全： - 假设 1（Poisson 测量模型）：\(X_{ij} \sim \text{Poisson}(\gamma_i \Lambda_{ij})\)，且所有 \(X_{ij}\) 跨细胞跨基因相互独立。这是最核心假设。统计含义：UMI 计数的噪声纯粹由泊松抽样产生，无过度散布或零膨胀。相比已有文献（如 Wang et al. 2017 的 DESCEND 允许过度散布，Van den Berge et al. 2019 的 tradeSeq 使用负二项），本文强化了分布假设以换取拆分的精确独立性。 - 假设 2（Size factor 已知或准确估计）：\(\gamma_i\) 不依赖于 \(X_{ij}\)。若 \(\gamma_i\) 需从数据估计，本文要求估计 \(\gamma_i\) 的步骤不使用基因 \(j\) 的数据（例如用全基因集的总量估计），以保证拆分后的独立性不被破坏。 - 定义（Count splitting 操作）：对每个 \(X_{ij}\)，独立生成 \(X_{ij}^{(1)} \sim \text{Binomial}(X_{ij}, \delta)\)，\(X_{ij}^{(2)} = X_{ij} - X_{ij}^{(1)}\)。其中 \(\delta \in (0,1)\) 是预设的拆分比例。 - 推论（拆分后的分布）：\(X_{ij}^{(1)} \sim \text{Poisson}(\delta \gamma_i \Lambda_{ij})\)，\(X_{ij}^{(2)} \sim \text{Poisson}((1-\delta)\gamma_i \Lambda_{ij})\)，且 \(X^{(1)}\) 矩阵与 \(X^{(2)}\) 矩阵完全独立。
主要结果 1. 定理（Type 1 误差控制）：在假设 1 与 2 下，若使用 \(X^{(1)}\) 估计隐变量 \(\hat{Z} = \hat{L}(X^{(1)})\)，使用 \(X^{(2)}\) 进行任意检验（如 GLM 似然比检验）检验 \(\hat{Z}\) 与基因 \(j\) 的关联，则在 \(\hat{Z}\) 的条件下，该检验的 Type 1 误差严格等于名义水平 \(\alpha\)。直觉：因为 \(X^{(2)}\) 与 \(\hat{Z}\)（仅依赖 \(X^{(1)}\)）独立，问题退化回经典的“固定协变量下的检验”。必要条件：Poisson 假设与 size factor 的可分离性。 2. 功效的渐近行为（模拟展示）：由于 \(X^{(2)}\) 的期望只有原数据的 \((1-\delta)\) 倍，检验功效必然下降。作者通过模拟表明，当 \(\delta=0.5\) 时，功效损失在中等样本量下可接受；且 \(\delta\) 的选择在 Type 1 误差控制上无差别，只在功效上有权衡。 3. 对聚类与轨迹推断的适用性：作者分别展示了在 K-means 聚类（离散 \(Z\)）和 Monocle3 轨迹推断（连续 \(Z\)）下，count splitting 均能恢复 Type 1 误差，而未拆分数据的检验 \(p\)-值严重偏向 0。
证明路线与技术技巧 - 整体路线： 1. 从模型出发：写出 \(X_{ij} \sim \text{Poisson}(\gamma_i \Lambda_{ij})\) 的独立结构。 2. 应用泊松-二项拆分：对每个 \(X_{ij}\) 施加 Binomial 拆分，利用泊松分布的卷积/分解性质，得出 \(X^{(1)}\) 与 \(X^{(2)}\) 独立且分别服从 \(\text{Poisson}(\delta \gamma_i \Lambda_{ij})\) 与 \(\text{Poisson}((1-\delta)\gamma_i \Lambda_{ij})\)。 3. 建立条件推断框架：将 \(\hat{Z} = \hat{L}(X^{(1)})\) 视为基于 \(X^{(1)}\) 的任意函数。由于 \(X^{(2)}\) 与 \(X^{(1)}\) 独立，\(X^{(2)}\) 在给定 \(\hat{Z}\) 下的条件分布，等于 \(X^{(2)}\) 在给定真实参数 \((\gamma_i, \Lambda_{ij})\) 下的边际分布（因为 \(\hat{Z}\) 不提供关于 \(X^{(2)}\) 的额外信息）。 4. 回归经典检验：在 \(X^{(2)}\) 上应用标准 GLM 检验，由于 \(X^{(2)}\) 仍服从 Poisson 分布（仅 size factor 变为 \((1-\delta)\gamma_i\)），经典检验的零分布成立，Type 1 误差控制得证。 - 关键跳跃点：从“拆分后的边际分布是 Poisson”到“拆分后的两部分联合独立”。这是整个证明最吃劲的地方。若数据不是 Poisson（如负二项），边际拆分仍可做（Binomial thinning），但联合独立性立刻崩溃，后续所有条件推断失效。 - 技术技巧点名： - Poisson-Binomial 分解：核心工具，用在每个 \(X_{ij}\) 上，保证拆分的独立性与分布保持。 - 条件化：将 \(\hat{Z}\) 条件化，把随机协变量问题转化为固定协变量问题。这比 Lee et al. (2013) 的选择性推断条件化更粗暴但也更通用——不需要写出 \(\hat{Z}\) 的选择事件，只需利用 \(X^{(1)}\) 与 \(X^{(2)}\) 的独立性。
真实例子与应用 - 数据：分化干细胞数据（pluripotent stem cells differentiating to cardiomyocytes），使用 Monocle3 估计连续伪时间轨迹。 - 如何用上去：对原始计数矩阵应用 count splitting（\(\delta=0.5\)），用 \(X^{(1)}\) 运行 Monocle3 的降维、构图、排序全流程得到 \(\hat{Z}\)（伪时间），用 \(X^{(2)}\) 对每个基因拟合泊松 GLM（以 \(\hat{Z}\) 为协变量，\(\gamma_i\) 为 offset），检验 \(\hat{Z}\) 的系数是否为 0。 - 结果：未拆分数据下，大量基因的 \(p\)-值极度偏向 0（假阳性膨胀）；count splitting 下，\(p\)-值在零假设下均匀分布，Type 1 误差恢复；同时，已知与心肌分化相关的关键基因（如 MYH6, MYH7）仍被检测出显著性（功效保留）。 - 想说明什么：展示 count splitting 在极其复杂的黑箱算法（Monocle3 包含 PCA、UMAP、图学习、伪时间排序等多步非线性操作）下依然有效，这是条件选择性推断路线无法做到的。
🔎 结论是否比证明窄 - 作者在正文严格证明了 Poisson 假设下的 Type 1 误差控制，但在讨论部分泛泛 claim：“count splitting 可能对偏离 Poisson 的数据也有一定稳健性”。这个 claim 无理论支撑，仅在模拟中展示了轻微过度散布下 Type 1 误差仍可控。这是一个典型的“条件 X 下严格证明，却被泛泛 conjecture”的地方，扎根在原文对模拟结果的解读中。
四、开放问题（点到为止）¶
- Poisson 假设偏离时的 Type 1 误差界：当 \(X_{ij}\) 服从负二项或零膨胀泊松时，count splitting 的联合独立性崩溃，Type 1 误差膨胀的定量界是什么？扎根在本文讨论部分对“robustness to model violations”的 conjecture，以及 Sarkar & Stephens (2021) 与 Van den Berge et al. (2019) 之间对分布假设的张力。
- 过度散布数据的拆分机制：能否构造一种类似 count splitting 的拆分，在负二项分布下保持两部分独立？扎根在本文 Section 3 明确指出的 limitation：“Our framework relies on the Poisson assumption”。
- 最优 \(\delta\) 的选择与功效-误差权衡的渐近理论：本文仅通过模拟展示 \(\delta\) 的影响，缺乏渐近功效的解析表达。扎根在本文对 \(\delta\) 选择的纯经验讨论。
- 多重检验下的 FDR 控制：本文仅保证单检验的 Type 1 误差，对 \(p\) 个基因同时检验时的 FDR 控制机制未触及。扎根在 Lähnemann et al. (2020) 提出的“grand challenge”中对多重比较的强调。
提醒：要确认某条是不是真 gap，去读同子领域近期约 5 篇的 intro——都指向 Poisson 假设的局限 = 共识（真 gap），互相打架 = 机会。
Maintained by 陈星宇 · Homepage · Source on GitHub