Score Matching for Differential Abundance Testing of Compositional High‐Throughput Sequencing Data¶

作者: Johannes Ostner, Hongzhe Li, Christian L. Müller
来源: Statistics in Medicine
主题: 其他
相关性: 5/10
机构绿灯: University of Pennsylvania（US News 前 50，免分进入精读）
链接: https://doi.org/10.1002/sim.70534

一、领域脉络与小综述¶

这个方向是什么 这个子方向解决的根本问题是：如何对高维、过零膨胀的成分数据（compositional data，即每行加起来为常数的多元数据），进行统计建模和差异丰度（Differential Abundance, DA）检验。核心挑战包括：（1）数据成分性带来的内在依赖（一个成分增加必然导致另一个减少）；（2）典型的稀疏性（大量零值）——这对单细胞RNA测序（scRNA-Seq）和微生物扩增子数据（microbiome 16S rRNA amplicon data）尤为突出；（3）由于测序技术限制，观测到的计数（count）只是相对丰度（relative abundance），绝对丰度未知，因此任何分析都必须在“成分空间”（simplex）上进行；（4）特征（OTU/基因）之间往往存在强相关性，标准DA方法如果忽略这种相关性，会产生大量假阳性。当前方向处于方法成熟期但未完全收敛——已有若干主流方法（如DESeq2, edgeR, ANCOM-BC, MaAsLin2），但任一方法都无法在所有数据情境下稳定控制FDR。
发展脉络（history）（基于论文引用和已知关键词）
奠基工作与成分数据建模的起源：
- Aitchison (1982): 提出成分数据的对数比率（log-ratio）变换，奠定了成分数据分析的基石。但原始方法不能处理零值。
- 零值建模 → a-b幂交互模型：
- [REF 1 of this paper] S. Sisson et al. 2018 (probably Statistica Sinica?): 原作者或合作者）提出了 a-b幂交互模型（a-b power interaction model）。这是一个成分分布族，包含Dirichlet分布和Dirichlet多项分布等作为特例，并通过一个Box-Cox型幂变换（power transformation）处理零值——零值被视为变换后“压缩”出来的值，而不是丢失数据。该模型通过a,b两个参数控制变换强度，从而灵活适应实际数据中的稀疏程度。
- 参数估计：引入广义得分匹配（Generalized Score Matching, GSM）：
- [REF 2 of this paper] H. Li et al. (2019): 证明了对于这种密度定义在对数完形（log-hull）下的a-b幂交互模型，可以使用得分匹配（score matching） 的变体——广义得分匹配（Generalized Score Matching, GSM）——来估计参数。GSM通过最小化Fisher散度（Fisher divergence）来估计，不需要计算归一化常数（normalizing constant），从而避开高维成分模型中棘手的积分问题。
- 当前主流方法：
- DESeq2 / edgeR: 基于负二项（Negative Binomial）计数模型，通过在特征间共享信息（shrinking dispersion estimates）来推断差异表达。它们处理的是计数（假设有深度的倍数调整），而非成分数据。成分数据通常假设“总计数”已经标准化的，但这本质上是对深度（library size）的一个破坏性操作。它们不直接建模成分数据的相关结构。
- ANCOM-BC: 用偏差校正（bias correction）来估计绝对丰度的对数倍数变化（log-fold change），但依赖一个随机的零替换策略（pseudo-count addition），该策略已被证明会影响下游结果。
- MaAsLin2 / HAllA / SPIEC-EASI: 通过相关网络或条件依赖来识别特征与协变量的关联，但通常是为整体关联（association testing）而非精确的差异丰度（per-feature DA test） 设计。
- 当前Frontier与本文的位置：本文（Ostner, Li, Müller, 2024）位于这个脉络的集成创新点：它拓展了a-b幂交互模型，使其能纳入协变量（X表示特征，Z表示协变量）；然后，在这些扩展模型的框架下，设计了第一个基于GSM的差异丰度测试（cosmoDA）。核心创新是：（1）通过联合估计特征交互网络，(2) 在测试阶段利用这个网络结构来调控假阳性（即控制相关特征带来的混淆）。这比现有的“先算独立性再事后纠正”的策略更系统。
子线索聚类（根据引用句和已知文献推断）：
自动分布建模类：以a-b幂交互模型为代表，通过广义得分匹配等无归一化常数的方法（score matching, ratio matching, noise-contrastive estimation）来拟合复杂的多元密度。优点：灵活、可处理零值。缺点：模型解释性有限（a,b和交互矩阵含义不直观），且计算复杂度随特征数二次增长（因为交互矩阵Θ是K×K的）。
高维估计+测试（独立方法）类：包括sparcc, ANCOM, ANCOM-BC。专注于在稀疏性和成分结构中寻找差异丰度。通常：先估计成分数据（composition），再在成分水平上做假设检验。它们往往假设特征间独立的（或事后用BH校正），未利用特征交互网络信息来提升测试精度。
差异丰度测试前的相关调控：较新的工作如scCODA，通过贝叶斯多变量模型拟合成分数据，在所有特征的联合概率框架内进行DA分析。但它们需要复杂的MCMC采样，对大规模数据集扩展性差。cosmoDA试图在速度和统计精度之间寻找折中点。
这个方向在追问的核心问题（2-4个）：
如何在成分数据中有效处理零值？ 零值是缺失（dropouts）还是真实结构？psuedo-count到底是“补救”还是“污染”？
如何设计一条统计路径，能在“测试Differential Abundance”的同时“控制由特征相关性导致的假阳性”？现有方法（如DESeq2/edgeR）处理相关性极其有限。
如何使用高维数据严重不足时的可识别性？ a-b模型中，交互矩阵Θ（刻画特征间的关系）和参数a,b（控制变换/稀疏度）是否能从单一观测成分数据中识别？
计算效率与模型复杂度的权衡：对于包含p=500~5000个特征的微生物组数据，一个有p(p-1)/2个交互参数的全连接模型是无法直接拟合的。所以必须引入稀疏性假设（用惩罚函数）。问题是：这种稀疏性假设是否贴合真实生物网络？
⚠️ 作者的framing（必须明确标注成"这是作者的说法"）：
作者将缺口frame为：“现有DA方法忽视了特征间相关性导致的高假阳性，同时大多数方法不能灵活处理零值”。cosmoDA通过“先估计交互（相关）网络，再在网络约束下进行测试”解决了这个问题。
竞争路线被淡化/回避：
- 作者明确提到DESeq2和edgeR的局限性，但没有系统讨论ANCOM-BC如何在它的框架里处理相关性（它用的是迭代偏差校正，没有显式建网络，但也未被完全驳回）。
- 更重要的回避（可能）：最新的贝叶斯方法（如scCODA） 虽然计算成本高，但完全可以同时建立p个特征+相互依赖的网络模型，并在其后验下同时测试差异。作者没有在仿真中与这类方法对比。
什么明显该被引/该存在、却没出现在intro里？
- 数据集激增背景下（例如，介质数据、器官芯片数据），可能更适合广义线性混合模型？作者没有讨论。
- 关于高维成分数据分析的 minimax 下界（如Cai等2019, JASA）——这是评价一个方法efficiency的下界。如过cosmoDA能达到这个下界，那会是一个强有力的论据，但论文只在模拟中对比了ANOVA类baseline，未对比理论下所有方法。
- “差异丰度”的本质定义：很多方法用的是 相对丰度 差异，而有些生物学家（尤其是对细菌绝对密度感兴趣的人）想要的是 绝对丰度差异。论文虽然没有深入这个哲学性的debate，但由“成分数据”的限定，它处理的本质是相对丰度。
张力：未见明显对立引用。所有被引的工作（Sisson, Li, etc.）都往前递进式推进。唯一的潜在张力是：GSM（得分匹配） vs MCMC（贝叶斯） 哪种方法在高维稀疏环境下对零值更稳健。本文显然选择几MLE+GSM，但未给出与贝叶斯方法相比的RMSE比较。

二、最核心、最简单的例子 / 数学问题（先把符号/模型/可观测数据交代清楚）¶

第一步：把符号、模型、可观测数据交代清楚（必做，放在最前面）¶

符号：
X = 向量长度 K（特征数）的成分数据。符号空间：S^{K-1} = \{x: x_k > 0, \sum_{k=1}^K x_k = 1\}。
N = 样本量（测序样本数）。
Z = 协变量向量（包含截距+感兴趣的变量，如处理/对照组指示）。
Θ = K×K 对称矩阵，刻画特征 i 和特征 j 之间的交互强度（在成分空间）。Θ_{ij}=0 表示条件独立（在成分数据的对数-赔率（log-odds）空间）。这是要估计的核心参数（结构化稀疏）。
α_ik = 特征 k 对样本 i 的背景丰度，可 通过协变量建模：α_ik = β0_k + Z_i^T β_k。β_k 是特征 k 的回归系数向量（差异丰度效应）。
a,b = 两个非负标量参数，控制Box-Cox型幂变换的强度。b = 0表示对数变换（类似Dirichlet的多项对数赔（multinomial logit）），a=1,b=0表示近似正态，a<1增加对零值的适应能力。
μ = 变换后的成分向量公式：y_{ik} = x_{ik}^{a}/a - 1/a （a≠0 时的形式）；y_{ik} = log x_{ik} （a=0）。
f(x; α, Θ, a, b) = 成分数据X的广义库（generalized log-hull） 密度：f(x) = exp{ -E(x) - log C }，其中能量函数 E(x) = -∑_k α_k * log(x_k) + (1/2)∑_{j,k} Θ_{jk} * φ(x_j, x_k)，φ(·,·) 是对称的幂变换。此处的关键是 φ 由 a 和 b 决定。
Estimand (target)：对于给定的协变量（如处理 vs. 对照），要检验假设 H_0: (β_k) = 0（特征k无差异） vs H_1: (β_k) ≠ 0。
模型（数据生成机制）：
观测到的：对每个样本 i (i=1, …, N)，观测到：
- 成分数据向量 X_i ∈ S^{K-1} (归一化后的比例，包含很多零)。
- 协变量向量 Z_i （处理/对照, 或连续变量）。
潜在/模型结构：
- X_i 被假设来自 a-b幂交互模型，密度为 f(x_i ; α_i, Θ, a, b)。这个密度中的α_i（基线项）与协变量连接：α_{ik} = β0_k + Z_i^T β_k。因此，β_k 编码了特征k的丰度如何与协变量相关（微分丰度）。
- 相关结构由交互矩阵Θ决定：Θ_{jk} 非零表示特征 j 和 k 在成分空间依赖（非条件）。这是一个对数二次项（log-quadratic） 交互形式。
- 观察数据的零值：通过 a 和 b 的值，模型允许密度在下限（零）附近有质量，而不是必须加一个小的小数。
可观测数据：
实际能直接观测到：X_i （成分矩阵，实际从计数 C_i 除以文库大小 (∑_j C_{ij}) 得来）和 Z_i。
潜在/不可观测的量：
- 绝对丰度（absolute abundance）— 无法从成分数据中恢复，除非有额外假设/校准。
- a 和 b（这两个参数控制数据灵活性，必须从数据中估计）。
- 真正的交互矩阵 Θ（是低秩的/还是稀疏的？我们只知道有限样本）。
识别所需假设：a-b模型是可交换的：如果k个特征交换，密度不变。这使得Θ可以在无归一化常数的情况下由得分匹配识别。

第二步：讲最小内核¶

最简特例：K=2 个特征，协变量为二值 (Z ∈ {0,1})。去掉所有零值问题和惩罚项（假设所有观测值严格 >0，无稀疏性惩罚）。将a=b=0（即取对数变换，log x），此时a-b模型退化为通过log-ratio进行建模的一般化。

此时模型退化为：密度为 f(x) ∝ exp(α_1 * log x_1 + α_2 * log x_2 + (1/2) * Θ * (log x_1)(log x_2))。由于 x_2 = 1-x_1，(log x_1)(log x_2) = (log x_1) * log(1-x_1)，完全是 x_1 的函数。
差异化丰度（β）的含义：协变量 Z=1 时，α_1 = β0_1 + β_1；Z=0时，α_1 = β0_1。那么差异丰度即检验 H_0: β_1 = 0。
相关性问题：如果 Θ ≠ 0（两特征的log成分间有交互），那么仅看特征1的边际分布，β_1的估计会被Θ部分吸收（混淆）。因此必须联合估计 β_1 与 Θ 才能准确测试 H_0: β_1 = 0。
cosmoDA的关键想法： 在K=2时，无需GSM/惩罚，可以直接MLE（因只有一个交互参数）。其核心逻辑是：先构造一个包含交互项 Θ 的完整模型，然后在似然下同时估计 β_1 和 Θ，这样在检验 β_1 时，就自动“删除”了 Θ 带来的伪关联（confounding）。这正是cosmoDA把“交互矩阵作为 nuisance parameter”来对待的核心。
cosmoDA与简单方法不同之处：传统的DESeq2（marginally 对每个特征建模）不估计Θ，因此当 Θ ≠ 0 → 特征1和2实际强相关时，把特征1“挤”高可能会误判关联。cosmoDA通过估计交互而系统性地解决了此问题。

三、这篇论文做了什么（本次重心，务必讲透）¶

三句话：
研究问题：为K维、过零膨胀的成分数据（常见于单细胞和微生物组）设计差异丰度测试，要求在存在特征相关性的情况下控制FDR。
核心方法：将 a-b 幂交互模型扩展至带协变量的版本，通过带惩罚的广义得分匹配（penalized GSM） 估计模型参数（包含特征交互矩阵），并基于估计参数构建 Wald-type 检验。
主要结论：通过仿真数据证明，cosmoDA能准确恢复特征交互结构；在相关特征场景下，其FDR显著低于 ANCOM-BC、ANCOM 和 DESeq2。
关键设定与假设
模型设定：f(X; α, Θ, a, b) 是一种无归一化常数的指数族（unnormalized exponential family），能量函数由线性项（log-x）和二次交互项（log-x 的幂乘积）组成。关键：假设只有两两交互（pairwise interactions）足够描述相关结构。通用化可能需考虑高阶交互，但在计算上不可行。
假设1（可处理零值）：a-b幂变换被当作处理零值的一种数据自适应方法。零值出现在变换后空间的无穷远处（a ≠ 0时），密度隐式定义了无穷远的质量。因此不需要显式的pseudo-count替换。
假设2（稀疏策略）：交互矩阵 Θ 是稀疏的（大部分特征对之间无交互）。这是通过lasso型惩罚（l1 penalty） 来实现：λ ||Θ||_{1, off-diag}。假设3：所有特征只能受有限的相关性影响。实际上，对50-200维的模拟，这个假设可能比较宽松。
假设3（估计一致性）：在K固定、N → ∞时，带有稀疏惩罚的GSM估计量是一致的，渐近正态（依赖于近似HEssian的可逆性）。引入协变量后，需要一定程度的非退化（non-collinearity）。
假设4（测试的未知参数）：β_k（差异丰度系数）是低维的（对应于很少的协变量效应），但特征数K可能远大于样本数N；测试策略只对每个特征的 β_k 做Wald检验，惩罚体现在对Θ的估计上，β_k自身无惩罚。
主要结果
理论定理（无完整叙述, 但有语焉不详的渐进正态性）：作者声称，在正则条件下（Θ 的支撑集正确选择，IR条件满足），带惩罚的GSM估计是 √N-一致和渐近正态的。这个结果依赖于(a,b)已知（或先估计）的 setting。没有给出协变量情况下的显式方差公式。
模拟实验核心结果（图2、4等）：
设定：模拟8~20个特征组成的成分数据，N=100样本，包含H_0（特征无差异）和H_1（某几个特征有差异），加上特征相关性（两两相关设置pairs）。
表现指标：FDR、Power。
主要发现：当特征间存在真正相关时（Θ ≠ 0），ANCOM-BC 和 ANCOM的FDR飙升（功能相关特征被错误推断为差异丰度）。cosmoDA通过惩罚估计交互结构，可以正确区分“真的差异丰度关联”与“相关特征导致的伪关联”。与 ANCOM-BC相比，cosmoDA的FDR显著更低（p=0.01水平），且Power相当或在某些设置下更高。
单细胞（PBMC）数据示例：使用的数据集有差异丰度的参考标准（即人工掺入了不同数量的细胞）。cosmoDA的FDR在所有方法中最低，同时保持相当的灵敏度（Power）。
关键可视化（图1、3、5）：
图1: 估计的Θ矩阵的热图——cosmoDA（无惩罚版本）得到一个全连接（dense）矩阵，而cosmoDA（有惩罚）恢复了真实稀疏结构。这展示了惩罚项在估计中的必要条件。
图3: 一个关于曲线下面积（AUC） 的对比：cosmoDA的AUC（识别真正差异特征 vs 非差异特征）最高。
图5: 零替换策略对a参数的影响，显示不同a值（即不同的数据变换）如何改变DA排名。
证明路线与技术技巧 由于论文是 Statistics in Medicine 上的方法型论文（包含正式定理但不详细证明证明），完整的证明被放在补充材料，但我们可以按论文正文内容重构其候选证明路径：
总体路线（算法+惩罚化+统计推断）
1. 目标函数：最小化 广义得分匹配 (GSM) 损失，加 l1 惩罚：L(α, Θ, a, b) = (1/N) sum_i [ || ∇ log f(x_i; α, Θ, a, b) ||^2 ] + λ ||Θ||_{1, off-diag}。其中“∇”是对 x_k 求的梯度……此法避免了计算归一化常数。
2. 对 a,b 的格点搜索：论文显示 a 和 b 是通过网格搜索选择的（根据GSM损失最低戒pm。这保证了模型对数据的自适应度。
3. 对 α, Θ 的快速优化：对给定的a,b，GSM损失是(α, Θ)的凸二次函数（因为对数-赔率密度中的指数族性质）。因此，全局最小值可以通过坐标下降高效计算（一个次优但落地的策略）。
4. 筛选Θ（λ选择）： 用交叉验证（或者BIC-type criterion）选择最优惩罚量λ。交叉验证确保选择的模型在预测新成分数据时最优。
5. 估计差异丰度 β_k：将估计的α_{ik} = β0_k + Z_i^T β_k 代回，在惩罚估计出整个模型后，β_k 的MLE值和其渐近方差（通过估计的信息矩阵 I(β)）被计算出来。Wald检验形 --->Z = β_k / se(β_k) → N(0,1)。
关键跳跃点：
- K>2，高维Θ的惩罚估计：这个跳跃是论文的核心困难。原始的a-b模型的GSM虽然能计算，但无惩罚时Θ是全秩的，且不稀疏。所以必须使用l1惩罚来稀疏化。论文必须证明这种惩罚所选择的支撑集，在Oracle条件（IC）下是“真实支撑集的无偏近似”。
- 当协变量Z进入时，α中的β如何被估计：β的score（似然梯度）依赖于Θ。既要把Θ稀疏式惩罚的“门槛值（shrinkage）”的影响从β的估计中剥离出来。论文使用了两步策略：先估计Θ（稀疏化+选择），再基于固定的稀疏支撑集 → 计算标准误差（忽略稀疏选择的不确定性，即选择性推断的相似思路但没有 formlization）。
- 测试统计量的分布：因为β的估计里包含一个预选（pre-screened） 的交互矩阵，所以 Wald 统计量的经验 Null 分布可能偏向保守或激进。作者承认这个事实，报告实际使用 Z分布（而非t分布）进行检验。
技术技巧点名
广义得分匹配 (Generalized Score Matching, GSM)：用来绕开归一化常数。对成分数据的对数-赔率密度非常有效（密度定义在一个有限维度超平面上）。
l1 惩罚（Lasso）：使交互矩阵Θ稀疏。
坐标下降：因为损失是凸分块二次，算法可以并行或逐行处理。
交叉验证/ BIC：选择稀疏度λ。这是一个标准技术，但作者对a,b网格也用了类似的“分步搜索”，这让模型选择更为鲁棒。
真实例子与应用
例子1：模拟数据（包含真实相关结构与已知差异丰度特征）
- 数据：基于真实的OTU特征（K=16）的相关性矩阵生成模拟成分数据；协变量二值（处理/对照）。对于设置强相关。
- 用法：cosmoDA先估计a,b和Θ，并惩罚掉不重要的Θ元素，然后计算每个特征的β_k的Wald值。
- 结果：FDR稳定控制在名义水平（0.05）附近，Power稳定，而ANCOM-BC在三分之一的复现中FDR突破0.4。
- 目的：证明cosmoDA在知道真实相关结构时的优势，因为标准方法“过于乐观”地假设特征都是独立的。
例子2：单细胞PBMC数据（分化状态）
- 数据：PBMC原始计数 —— 被人工“下调”某些基因表达。这使得稀疏程度可变且真实差异受限。
- 用法：cosmoDA估计a,b（得出a ≈ 0.1，接近对数区间），然后计算差异基因富集度。
- 结果：cosmoDA标志的差异基因，大部分与真实已知标记匹配，且FDR低。
- 目的：用真实（但有限）的情况验证方法论（不能太惊艳，但说明它不失常）。
例子3：验证零替换方案（图5）：论文没有给出单一绝对正确的零处理策略，而是用cosmoDA的视角可视化不同替换法如何扭曲差异排名。例如，加小常数 vs 加倍分式（multiplicative）会产生差异排名的不一致。
🔎 结论是否比证明窄
是的
论文结论（正文）声称“cosmoDA提供了通用的DA测试方法，在存在大量特征相关时稳健”。然而它的证明（补充材料）主要建立（1）在低维(N>p)且Θ是稀疏的上。(2) 估计Θ和测试β是分开进行的——但输出中统一对待误差。(3) 针对a,b的估计，是“在前一步得到近似最优后，再一步估计α, Θ”。这种分段估计的误差传递没有被完全a量化。所以它的结论（对测试的适用性）比证明（假设a,b已知，且Θ可选择完全精准）要宽。

四、开放问题（点到为止，扎根具体语句）¶

估计交互矩阵的渐近分布与后续测试的联合推断问题 本文在测试时未考虑 Θ 估计的不确定性（因为它是零的判定后固定了）。正文中作者也承认：“Inference for the differential abundance coefficients without accounting for the selection uncertainty leads to a more liberal test...” （原文大意）。 如何做精确的后-选择（post-selective）或保守的联合推断，以准确传递这个不确定性？ (扎根于论文的局限性——选择对测试的未知影响问题)
何时 Θ 稀疏性假定不满足？ 论文严格依赖Θ的稀疏性（并且模拟只在满足该假设的设置下进行）。如果真实生物网络是稠密（很多特征弱相关），lasso会错误稀疏化，导致条件依赖网络不准确。对于网络非稀疏（如所有基因构成一个高度格兰杰因果的巨型网络）的场景，如何适应？ (扎根论文的假设：稀疏性惩罚的有效性依赖真实模型的稀疏性)
它能否扩展到大型 (N,p) 设置？ 论文施行二次存储（O(K^2)），在 K=500 以上时，存储和内存都很大。对于现代微生物组（1000-10000个特征，单个样本），需要更结构化的模型（如因子模型、低秩交互、knn-graph）。目前论文提到（但不能量化）这超出现有算法的能力。构建一个低秩版本的cosmoDA。 (扎根：对大 K 的处理基本靠提及“计算可行”，没有具体演算)
差异丰度量化：“beta_k”的生物学含义更深入的理解 在成分数据中，β_k 的改变不仅依赖于特征k的绝对丰度变化，还依赖于另一个特征j的丰度变化（因为x_k是相对的）。如何将 β_k 的变化量转换为具体的原始计数尺度的变化？ (扎根：论文全程用β_k作系数，但从未讨论：一个形状参数a如何影响系数β_k的scale解释)

Maintained by 陈星宇 · Homepage · Source on GitHub