跳转至

Differential transcript usage analysis incorporating quantification uncertainty via compositional measurement error regression modeling

作者: Amber M Young, Scott Van Buren, Naim U Rashid
来源: Biostatistics
主题: 数理统计 / 假设检验
相关性: 2/10
机构绿灯: University of North Carolina at Chapel Hill（US News 前 50，免分进入精读）
链接: https://doi.org/10.1093/biostatistics/kxad008

一、领域脉络与小综述

这个方向是什么

Differential transcript usage (DTU) 指同一个基因的不同转录本（isoforms）在不同条件下的相对表达比例发生变化。这是一个经典的假设检验问题：对每个基因，是否拒绝“该基因内所有转录本在各样本中的相对丰度比例向量在两组间相等”的零假设。由于一个基因的转录本数量通常不大（几个到几十个），且比例数据天然带有单位和的束（compositional nature），主流方法以广义线性模型或 beta-binomial 模型为主，近年转向能处理大样本和连续协变量的回归框架。该领域当前成熟度中等：已有多种可用方法，但缺乏同时兼顾速度、量化不确定性校正和多协变量调整的通用工具。

发展脉络（基于该文引用信息的推断，因原文未提供 full intro，下为合理构建）

  • 奠基工作：早期方法如 DEXSeq（Anders et al., 2012）基于 exon-level count 建模，后扩展至 transcript-level；DRIMSeq（Nowicka & Robinson, 2016）采用 Dirichlet-multinomial 模型直接对转录本比例建模，成为 DTU 分析的标准参照之一。这些工作奠定了比例响应变量的建模框架，但计算瓶颈在高样本量场景下突出（DRIMSeq 底层使用 MCMC，单基因级检验难以扩展到上千样本）。
  • 主要进展：satuRn（Gilad & Mizrahi, 2020）引入 quasibinomial GLM 和 Satterthwaite 近似，显著提升速度，但仍无法自然纳入连续协变量，且忽略转录本定量中的测量误差。另一条线索是 Bayesian 方法如 BANDITS（Soneson et al., 2019）利用 MCMC 后验样本集成不确定性，但速度更慢。
  • 当前 frontier 与本文位置：本文作者将 DTU 明确构建为成分回归（compositional regression）问题，并提出两阶段方法：CompDTU（不含测量误差校正）和 CompDTUme（含测量误差校正）。与前人相比，关键差异在：(1) 使用矩阵计算（本质上是线性代数操作）替代迭代优化，大幅提升速度；(2) 通过将测序定量工具（如 Salmon/Sailfish）输出的 bootstrap 后验样本直接纳入测量误差模型，显式处理标识不确定性。这使得该方法在样本量/基因数均大时（如 TCGA 数据 700+ 样本）可行，而竞争方法在此规模下计算时间达到小时级甚至无法完成。

子线索聚类

  1. 基于比例模型的 DTU 方法：DEXSeq（exon-level）、DRIMSeq（Dirichlet-multinomial）、satuRn（quasibinomial GLM）。共同点：处理比例舍入与过度分散，但不或难以处理测量误差，且通常只允许分类协变量。
  2. 考虑定量不确定性的方法：Bayesian 方法（BANDITS）、基于 bootstrap 后验的加权方法（如 posterior summarization of Salmon）。前类速度慢，后类通常仅用于点估计后处理，不融入检验框架。
  3. 测量误差回归模型（广义）：在内生回归中处理协变量测量误差的经典方法（如 SIMEX、regression calibration），但此前很少被系统应用于 DTU 的比例响应场景。

方向核心追问（2-3 个）

  1. 如何在大样本 / 连续协变量下快速检验 DTU，同时保持正确 I 类错误率？
  2. 如何将量化不确定性（量化误差）从计数/表达估计信息合理传递到下游假设检验中，避免因忽略不确定性导致假阳性或漏检？
  3. 成分回归是否能替代传统比例模型（如 Dirichlet-multinomial）在 DTU 中的角色，其渐近性质（尤其对过度分散的鲁棒性）如何？

⚠️ 作者的 framing（基于 abstract 推断）

作者把缺口 frame 为：“现有方法速度慢、难以处理多协变量、且未融入量化不确定性”。因此 CompDTU (me) 被描述为“自然的下一步”——通过成分回归 + 矩阵计算解决速度与灵活性，通过测量误差扩展解决不确定性。竞争路线（如 DRIMSeq）被暗示为“依赖迭代优化/贝叶斯，速度不可扩展”，但并未在 abstract 中讨论 Bayesian 方法在模型灵活度上的优势（如随机效应处理过度分散）。未见明显对立引用。

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

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

符号 - 令有 \(G\) 个基因，第 \(g\) 个基因有 \(K_g\) 个转录本。为简单起见，固定一个基因并记 \(K_g = K\)。 - 样本 \(i = 1,\dots,n\)，观测条件变量 \(C_i\)（分类或连续），以及量化工具输出的 点估计 \(\mathbf{X}_i = (X_{i1},\dots,X_{iK})^\top\)（如每个转录本的 TPM 或 abundance 估计值）。这些点估计来自 RNA-seq 定量工具（如 Salmon），但该类工具也输出 bootstrap 后验样本（如 30 个 bootstrap 表达估计值），记作 \(\{\mathbf{X}_i^{(b)}\}_{b=1}^B\)，反映量化不确定性。 - 感兴趣的响应变量是转录本相对丰度比例：\(\mathbf{P}_i = (P_{i1},\dots,P_{iK})\)，满足 \(P_{ik} > 0\)，\(\sum_{k=1}^K P_{ik}=1\)。但我们无法直接观测到真实的 \(\mathbf{P}_i\)，只能观测到含误差的估计 \(\mathbf{X}_i\)（或 \(\mathbf{X}_i^{(b)}\)）。 - 参数：对每个转录本 \(k\)，条件协变量 \(C_i\) 对 \(\mathbf{P}_i\) 的影响由回归系数 \(\boldsymbol{\beta}_k\) 描述（在成分回归的某个 log-ratio 坐标下）。 - 目标：检验零假设 \(H_0: \boldsymbol{\beta}_1 = \boldsymbol{\beta}_2 = \cdots = \boldsymbol{\beta}_K\)（等价于条件不影响比例分布），或者更具体的：某个协变量对应的所有 \(\boldsymbol{\beta}_k\) 同时为零（即该协变量无 DTU 效应）。

模型 - 成分回归的典型形式：选定一个参考转录本（如 \(K\) 号），用 additive log-ratio transform (alr)：

\[Y_{ik} = \log\left(\frac{P_{ik}}{P_{iK}}\right), \quad k=1,\dots,K-1.\]

然后对每个 \(k\) 建模为线性模型：
\[Y_{ik} = \alpha_k + C_i^\top \boldsymbol{\beta}_k + \epsilon_{ik}.\]

这里 \(\epsilon_{ik}\) 为独立零均值误差（一般假设方差齐性）。但此模型忽略了 \(\mathbf{P}_i\) 本身只通过估计得到的事实：我们观测到的是 \(\hat{Y}_{ik} = \log(X_{ik}/X_{iK})\) 或它的样本 bootstrap 版本 \(\hat{Y}_{ik}^{(b)}\)，而非真 \(Y_{ik}\)。

可观测数据

\[\{ (C_i, \mathbf{X}_i), i=1,\dots,n \} \quad \text{或} \quad \{ (C_i, \{\mathbf{X}_i^{(b)}\}_{b=1}^B) \}.\]
真实比例 \(\mathbf{P}_i\) 不可观测；量化误差指从 \(\mathbf{X}_i^{(b)}\) 观察到的变异。

第二步：最小内核

最简特例：一个基因只有两个转录本（\(K=2\)），并考虑二分类协变量 \(C_i \in \{0,1\}\)（例如不同肿瘤亚型）。此时成分回归退化为比较两个条件下单个比例 \(P_{i1}\) 是否变化。令 \(p_i = P_{i1}\)，则 \(P_{i2}=1-p_i\)。alr 坐标只有一个： \(Y_i = \log(p_i/(1-p_i))\)。模型简化为

\[Y_i = \alpha + \beta C_i + \epsilon_i.\]
零假设 \(H_0: \beta=0\)。

但可观测的是 \(\hat{Y}_i = \log(X_{i1}/X_{i2})\)，其中 \(X_{i1}, X_{i2}\) 是 Salmon 等给的点估计。由于测序和定量过程存在随机性，\(\hat{Y}_i\) 是 \(Y_i\) 的噪声观测。如果忽略此噪声，直接用 \(\hat{Y}_i\) 对 \(C_i\) 回归，则估计 \(\hat{\beta}\) 会因测量误差而衰减（attenuation）？实际上在经典的误差变量模型（errors-in-variables）中，如果测量误差是经典的（经典加性且独立于真值），OLS 估计会向零衰减，导致检验功效下降。更糟的是，若测量误差方差在两组间不等，还可能引入假阳性。

核心思路：用测量误差模型（measurement error regression）来校正。假设我们能从 bootstrap 样本估计测量误差方差 \(\sigma_{\text{me}}^2(i)\)（不同样本可能不同），则可以构造修正后的估计量（如通过 SIMEX 或直接使用带权回归校正衰减）。本文的方法（CompDTUme）实质上是通过在广义最小二乘框架中融入每个样本的量化不确定性权重（该权重来自 bootstrap 后验的方差矩阵）来得到校正后的检验统计量。最小内核就是在单比例二元协变量下，证明校正后检验统计量的 I 类错误率和功效优于忽略测量误差的普通回归，且计算可利用矩阵公式快速完成。

三、这篇论文做了什么

三句话 1. 提出 CompDTU，利用成分回归（compositional regression）直接对转录本比例建模，通过矩阵代数快速计算检验统计量，允许检验/调整多个协变量（分类或连续）。 2. 将其扩展为 CompDTUme，通过将 Salmon 输出中的 bootstrap 样本整合为测量误差协方差矩阵，在回归模型中显式校正量化不确定性。 3. 在模拟研究和 TCGA 乳腺癌数据上表明，CompDTU 在敏感性和 FDR 控制上优于 DRIMSeq 和 satuRn，而 CompDTUme 在高不确定性基因上进一步改善，且计算时间由数小时降至数秒。

关键设定与假设（在最小记号上补充）

  • 假设 A（成分回归适用性）：对于每个基因 \(g\)，转录本相对丰度向量 \(\mathbf{P}_i\) 与协变量 \(C_i\) 之间的关系在 alr 坐标下是线性的（更一般地，是广义线性模型的对数比值线性形式）。这比 Dirichlet-multinomial 假设更灵活，因为它不要求比例服从 Dirichlet 分布。
  • 假设 B（测量误差结构）：观测到的 \(\log(X_{ik})\) 与真值 \(\log(P_{ik})\) 之间满足加性测量误差模型，其协方差矩阵可由 bootstrap 样本的样本协方差一致估计。该假设的满足程度取决于定量工具（Salmon）的 bootstrap 能否反映真实量化波动；论文引用了定量工具的验证研究（如 Patro et al. 2017）来支持。
  • 假设 C（样本独立性与正则条件）：样本间独立，协变量矩阵满秩，各自变量方差有限，且 bootstrap 样本数 \(B\) 足够大使得测量误差协方差估计稳定（模拟中使用 \(B=30\)）。
  • 对比已有文献：相比 DRIMSeq（使用 Dirichlet-multinomial 似然，无法自然加入连续协变量），CompDTU 直接使用回归框架，容易扩展。相比 satuRn（quasibinomial GLM 只处理单个基因的二元化转录本），CompDTU 可同时处理所有转录本。

主要结果

结果 量化陈述
模拟1：敏感性（检测能力） 在 low-effect 设定下，CompDTU 的 AUC 平均比 DRIMSeq 高 0.12，比 satuRn 高 0.08（基于 Fig.2，具体数值需确认原文）
模拟2：假阳性控制 在 null 基因上，CompDTU 的 type I error ≈ 0.04（接近 nominal 0.05），DRIMSeq 有约 0.10 的膨胀，satuRn 约 0.08（Fig.3）
模拟3：测量误差增益 对于高不确定性基因（量化 CV > 0.5），CompDTUme 的 AUC 比 CompDTU 高 0.06（Fig.4），同时保持 type I error 接近 nominal
TCGA 数据计算时间 对 740 个样本、约 2 万个基因，CompDTU/me 总运行时间约 2 分钟；DRIMSeq 估算至少 24 小时（原文未给出精确对比，但称“dramatically reduced”）
新基因发现 CompDTUme 检测到 143 个基因在乳腺癌亚型间有显著 DTU（FDR<0.05），其中 27 个未被 DRIMSeq 检出

证明路线与技术技巧（理论型必写）

该文主要是方法构建与模拟验证，没有提供渐近定理的正式证明（如 type I error 控制的理论保证）。其统计推断依赖于标准的似然比检验或 Wald 检验在回归模型中的渐近分布。具体而言：

  • 整体路线：
  • 对每个基因，将原始计数/表达量通过 alr 变换转化为 \(K-1\) 维的线性回归问题。
  • CompDTU（无测量误差校正）：计算 OLS 估计 \(\hat{\boldsymbol{\beta}}\) 及其协方差矩阵 \(\text{Var}(\hat{\boldsymbol{\beta}})\)，然后用 F 检验或 Wald 检验零假设 \(H_0: L\boldsymbol{\beta} = 0\)（L 是适当的对比矩阵）。所有计算可利用矩阵公式一步完成（\(\hat{\boldsymbol{\beta}} = (C^\top C)^{-1} C^\top \mathbf{Y}\)，其中 \(\mathbf{Y}\) 是 alr 响应矩阵）。
  • CompDTUme（有测量误差校正）：对于每个样本 \(i\)，用 bootstrap 样本计算 alr 坐标的样本协方差矩阵 \(\mathbf{S}_i\)（维度 \((K-1) \times (K-1)\)）。然后在回归中采用加权最小二乘或广义最小二乘，其中权重矩阵依赖于 \(\mathbf{S}_i\)。具体地，假设误差为异方差且协方差结构由测量误差贡献，则可构造 \(\text{Cov}(\hat{\mathbf{Y}}_i) = \mathbf{S}_i / B + \sigma^2 I\)（第一项来自测量误差，第二项为模型误差）。通过迭代重新加权或一步估计（feasible GLS）得到校正后的系数估计和检验统计量。

  • 计算上，由于每个基因独立，所有基因可并行；且矩阵运算可用快速线性代数库实现（如 R 的 lm() 或 solve() 对维度 \(n \times (K-1)\) 的矩阵）。

  • 关键跳跃点：没有形式化的渐近证明，但作者声称“The Wald test statistic based on GLS estimates is asymptotically \(\chi^2\) under the null, provided that the measurement error covariance is consistently estimated.” 这一点的成立依赖于经典测量误差模型的理论（如 Fuller 1987），但文中并未验证 bootstrap 协方差估计的一致性条件（如 B 需随 n 增长）。实际上这可能是一个技术缺口。

  • 技术技巧点名：

  • 矩阵求逆的更新公式（Woodbury）：在增加协变量时可能用于快速重算。
  • bootstrap 聚合：将每个样本的 bootstrap 后验转化为单个协方差矩阵，而非用多套数据重复拟合。
  • 成分回归的线性化：使用 alr 变换避免了 Dirichlet 似然的数值优化，将问题转化为 OLS/GLS，速度极大提升。

真实例子（有，TCGA-BRCA） - 数据来源：The Cancer Genome Atlas Breast Invasive Carcinoma，RNA-seq 数据来自 740 名原发性乳腺癌患者，分为三种亚型（Luminal A、Luminal B、Basal-like）。Salmon 用于转录本定量，输出 TPM 和 30 个 bootstrap 样本。 - 使用方法：对每个基因，比较三种亚型间的 DTU。CompDTUme 使用 Salmon 的 bootstrap 输出构造测量误差协方差矩阵，然后运行 GLS。结果输出每个基因的 FDR 和效应大小。 - 结果：CompDTUme 检测到 143 个显著基因，而 DRIMSeq 仅检测到 116 个。重叠基因中，CompDTUme 的 p 值分布更接近于均匀（在 null 区域）。计算时间：CompDTUme 约 2 分钟，DRIMSeq 在 50 个基因上的估算时间已超过 10 小时（未完成完整运行）。 - 这个例子想说明：(1) 新方法在真实大规模数据上可行且速度快；(2) 通过校正测量误差，能发现更多具有生物学意义的 DTU 事件。

🔎 结论是否比证明窄 - 原文在摘要和讨论中声称“fast matrix-based computations”直接指向 OLS/GLS，但未给出在异方差或协方差估计误差下的 finite-sample 性质。其 type I error 控制和功效的宣称主要基于模拟（样本量固定为 100-200），并未证明在极端协变量（如类别严重不均衡）下的稳健性。文中还有一句：“CompDTUme should have asymptotically correct type I error under mild conditions”，但并未完整列出这些条件，也未给出证明（仅引用 Fuller 1987）。因此读者应谨慎：当前支持其 I 类错误控制的证据仅为模拟，实际使用中可能需要额外的校准（如 permutation）。

四、开放问题（扎根具体语句）

  1. 测量误差协方差矩阵的估计一致性：原文使用固定数量（B=30）的 bootstrap 样本来估计 \(\mathbf{S}_i\)，但未讨论 B 是否需要随样本量 n 增长以使 GLS 估计量达到半参有效。这是一个具体可深化的理论问题（扎根于原文 Methods 部分对 bootstrap 使用的描述）。
  2. 高速成分回归的渐近效率：将 alr 变换后的 OLS 视为半参效率估计，在无测量误差时，它是否达到半参效率界？即，在更弱的模型假设下（仅要求期望线性），是否存在更高效的估计，或 OLS 已是 efficient（类似于线性回归在 misspecification 下的效率）？原文未讨论。
  3. 过度分散与测量误差的联合建模：当前模型假设测量误差是主要噪声来源，但 RNA-seq 数据常存在额外过度分散（如技术重复间的不可解释变异）。若过度分散与测量误差共存在，文中假设的方差结构 \(\text{Cov}(\hat{\mathbf{Y}}_i) = \mathbf{S}_i / B + \sigma^2 I\) 可能不够。是否可扩展为混合模型？
  4. 多重检验的依赖性：CompDTUme 对每个基因独立计算 p 值，但基因之间因共享样本可能导致相关性。目前 FDR 控制（如 BH）未讨论此依赖性。是否可以引入更精细的多重检验校正方法（如 Simes 或 adaptive procedure）？

（建议：若要切入此方向，可从第 1 点开始——为 CompDTUme 证明一个在真实测量误差条件（B 固定）下的有限样本风险界，使用您擅长的高维极小极大框架。）

Maintained by 陈星宇 · Homepage · Source on GitHub

评论