《Applied Microbiology》：The Use of Confidence Intervals in Differential Abundance Analysis of Microbiome Data
摘要
差异丰度分析（DAA）是微生物组研究中的一项关键任务，旨在识别能够可靠区分不同组别的微生物特征。研究表明，微生物组系统相对稳定且具有恢复力，但在特定条件下，即使是微小的变化也可能引发菌群失调。微生物数据集的高维度性加剧了检测此类变化的挑战，因为它带来了多重比较问题，需要进行假设过滤。标准的多重比较校正程序（如Benjamini-Hochberg程序）是为真实阳性数量较多的场景设计的，对于真实信号比例可能非常低的微生物组数据往往过于保守。因此，迫切需要针对微生物组数据量身定制的假设过滤方法。组间差异的置信区间（CI）为p值过滤提供了一个强大的替代方案，因为其范围同时传达了效应的显著性、潜在大小和方向，以及估计本身的确定性。微生物数据可以使用负二项（NB）分布进行充分建模，其位置参数可以使用Hodges-Lehmann估计量（HLE）进行稳健估计。利用合成和实验数据，我们证明了基于两样本HLE的CI进行假设过滤是一种比较微生物数据的稳健方法。分析表明，HLE-CI方法在微生物DAA中提供了与使用多重调整方法过滤相同的精度水平，同时实现了显著更高的召回率。本研究的结果表明，基于HLE-CI的过滤可以成为寻找微生物组生物标志物的有效步骤。
1. 引言
1.1. 背景
近年来，人们广泛认识到人类微生物组显著影响宿主的消化、免疫、代谢甚至神经功能，使其平衡对整体健康至关重要。微生物组研究的一个核心任务是差异丰度分析（DAA），旨在识别区分微生物组数据集中不同组别的可靠标记特征。微生物数据集的特点是高维度，通常包含数百甚至数千个特征，例如细菌分类群或基因。生物标志物发现通常涉及测试每个特征与统计关联。由于每个测试都带有非零的假阳性（I类错误）概率，进行数千次测试会产生偶然发现显著结果的高风险。因此，微生物分析本质上提出了一个多重比较问题，凸显了假设过滤的必要性。
鉴于这种背景，微生物数据的高维度性为假阳性创造了充足的机会。然而，与所检查的分类特征总数相比，已识别的微生物标记数量通常少得惊人，并且绝大部分（>80%）的微生物变异始终无法由研究的人口统计学、临床或环境因素解释。这种差异在很大程度上可以通过个体内变异和微生物组系统固有的生物学稳定性来解释。
研究表明，人类微生物组系统通常是稳定和具有恢复力的。例如，纵向研究表明，在重复检查患者时，独特微生物物种的重新识别率很高（超过86%）。类似地，核心微生物组物种在大量人群中持续被识别，并解释了相当大比例的微生物变异。同时，解释的微生物变异性的最大部分归因于个体内差异，这使模式识别复杂化。然而，菌群失调可能由已知的细微微生物失衡引起。微生物组的微小变化可能引发群落范围的结构、代谢活性或两者的变化，最终导致显著的表型效应。
因此，微生物DAA的方法学挑战之一是在具有较大个体内变异、但总体组成似乎稳定的微生物组中检测细微变化。检测微生物组成中这些变化的挑战有三方面。首先，条件致病菌的增加或关键共生菌的减少通常是轻微的。其次，条件致病菌和关键共生菌通常只占总微生物种群的一小部分。第三，它们相对丰度（本身已经很小）的变化也可能很小。在对所有分类特征进行测试并针对多重比较进行调整后，这些细微但具有生物学重要性的变化可能会重新表现为不显著。
这种掩盖的原因在于多重比较校正的性质。尽管这些假设过滤方法旨在控制假阳性（I类错误），但它们固有地降低了统计功效。这对于高维度数据（如微生物数据集）可能尤其成问题，因为大多数特征不是差异丰度的（即“稳定”）。在这种情况下，标准校正变得过于保守，不成比例地惩罚了少数真实信号，导致假阴性（II类错误）。鉴于这一局限性，可以得出结论，微生物DAA需要其他假设过滤方法。
另一种流行的假设排序方法是按效应大小过滤。大多数效应大小度量（例如，来自r和d家族的）提供了关于效应大小的信息，补充了p值所指示的统计显著性。具有合理样本量的大效应大小通常比小的p值更能揭示实际重要性。然而，效应大小度量也带来了自身的挑战。例如，尽管像LEfSe（线性判别分析效应大小）这样的工具在微生物DAA中非常流行，但它们容易产生或夸大假阳性。这种脆弱性出现是因为效应大小，无论是来自LEfSe的LDA估计还是其他常见统计量如Cliff‘s δ，都是单点估计，旨在补充或扩展基础显著性检验。因此，这些点估计仍然高度依赖于样本特性。
解决这个问题的一个众所周知且直接的方法是考虑置信区间（CI）。组间差异的CI可以被视为原始p值和相关点估计统计量的更好替代方案，因为它们提供了位置偏移值的合理范围。该范围同时传达了显著性、效应的潜在大小和方向，以及估计本身的确定性。
一个关键问题是什么是分析微生物数据的最佳CI。由于理想的CI大致描述了两个样本中心趋势之间真实偏移的合理范围，第一步可以是获得微生物数据本身位置参数（中心值）的稳健估计。
微生物组样本典型的过度离散、零膨胀计数数据可以使用负二项（NB）分布进行充分建模。该分布模拟在一系列独立伯努利试验中，在达到指定数量的成功之前的失败次数。NB分布的一个定义特征是其右偏、单峰形状和长右尾，这自然捕获了微生物丰度数据中观察到的过度离散（方差超过均值）。分布的形状，特别是其偏斜程度和峰的位置，由其参数控制。随着所需成功次数的增加，分布变得不那么偏斜，其形状接近对称的钟形曲线。因此，具有显著正偏斜的NB分布最适合模拟低丰度、稀疏的分类群，而接近对称钟形的分布更好地描述了常见的、丰富的分类群。悬而未决的问题是：哪种集中趋势的度量最能准确代表NB分布微生物数据的中心？
均值是最常用的集中趋势度量。在NB分布的背景下，它直接对控制分布扩散和位置的参数敏感。它也与微生物数据的零膨胀性质很好地吻合，使其成为反映分类群丰度和普遍性的理论上有用的度量。然而，微生物数据集经常包含异常值，这在实践中使均值成为一个不可靠的统计量，容易受到极端值的过度影响。相关的稳健度量，如修剪平均值和缩尾平均值，对极端值提供了更大的抵抗力。然而，这些方法引入了一个新的挑战：确定一个最优且生物学上合理的修剪或缩尾阈值。此外，即使应用了这些调整，所得估计的可靠性也无法保证。中位数定义为排序向量的中间值，比均值更不易受异常值影响，为样本比较提供了更大的稳健性。然而，它可能是NB分布真实中心参数最不准确的理论估计，因为它没有考虑右偏斜。此外，在零膨胀微生物数据的背景下，中位数往往低估了真实的集中趋势，并且对于低普遍性的样本很容易产生零估计。
传统度量的一个稳健替代品可以是伪中位数，通过Hodges–Lehmann估计量（HLE）计算。该估计量对NB分布的右偏斜敏感，同时不会受到极端异常值的过度影响。在单样本情况下，HLE定义为所有成对平均值（也称为Walsh平均值）的中位数。在对称分布中，HLE等于中位数。然而，对于右偏斜的NB数据，HLE位于中位数和均值之间，有效地平衡了中位数的稳健性和均值的敏感性。这种平衡是通过所有成对平均值（类似于均值）合并有关数据分布的信息，同时通过取它们的中位数来保持对极端值的抵抗力来实现的。因此，HLE提供了一个有利的折衷方案。两样本HLE定义为两个样本中观测值之间所有成对差异的中位数。该估计量的CI可以通过反转Mann–Whitney U计数统计量导出，从而产生两个NB样本位置参数之间差异的非参数CI。
此外，与主要产生p值和相应检验统计量的常规测试不同，两样本HLE估计直接产生组间中心趋势偏移（即两个样本中所有观察到的成对差异的中位数）的直观且可解释的度量。因此，单样本HLE（伪中位数）和用于位置间差异的两样本HLE可能代表一种强大的、假设宽松且稳健的微生物数据分析工具，值得进一步研究。
1.2. 研究目的与目标
总之，在微生物DAA中，II类错误（假阴性）可能与I类错误（假阳性）同样关键。标准的多重测试校正倾向于在真实信号较少的数据集（如许多微生物组研究典型的数据集）中过度惩罚比较，通常达到禁止任何标记分类群检测的程度。
本研究的目的是证明基于HLE-CI的假设过滤在缓解使用多重测试校正进行假设过滤期间标记分类群发现过度惩罚问题方面的效用。
为此，该研究解决了两个具体目标。首先，我们使用合成数据调查标记发现错误率，比较几种方法：来自流行统计检验的原始和调整后p值，以及基于效应大小的过滤。其次，我们在真实世界数据集上对HLE-CI过滤性能进行了比较分析。
1.3. 结果总结
本研究的结果表明，在具有稀疏真实信号的数据集中，传统的多重比较调整会无意中过度惩罚和掩盖真实的生物信号。相比之下，基于HLE-CI的过滤在保持优异或相当的精度的同时，提供了显著更高的召回率。因此，我们表明HLE-CI实现了显著更高的F1分数，这一优势即使在样本量减少的情况下仍然存在。本研究的结果说明，基于HLE-CI的假设过滤可以成为微生物组数据中生物标志物发现的一种有益且强大的方法。
2. 材料与方法
2.1. 合成数据集
模拟数据通过使用Python v3.12中的NumPy v2.0.1包从NB分布中随机抽样生成。每个模拟的分布参数是随机抽取的：参数'n'（成功次数）从1到10的离散均匀区间中抽样，成功概率'p'从0到1的连续均匀区间中抽样。为确保可重复性，每次迭代的随机种子是固定的。每个模拟实验包含100,000次抽取，以最小化随机效应。
2.2. 实验数据集
炎症性肠病（IBD）数据集包括IBD患者及其健康亲属（HR）的16s微生物组分析数据。患者在哈萨克斯坦阿拉木图的临床地点招募。所有参与者都提供了书面知情同意书。总共招募了270名18岁以上的个体，形成了135个匹配对（135名IBD患者和135名健康亲属）。该数据集的收集和使用得到了纳扎尔巴耶夫大学“阿斯塔纳国家实验室”PI伦理委员会的批准（协议号#05-2022，2022年10月21日）。
生物信息学与实验数据的统计分析
使用LotuS2生物信息学管道处理测序数据。使用Lambda分类器和SILVA v138.1参考数据库对ASV进行 taxonomic 注释。在多样性分析之前，数据进行了稀释。使用Python v3.12中的NumPy v2.0.1, SciPy v1.15.1, scikit-bio v0.6.3, scikit-learn v1.6.1, Matplotlib v3.10.0和seaborn v0.13.2包进行生物多样性分析、特征重要性分析和可视化。仅考虑在任何组中50%样本中普遍存在的分类群进行分析。使用Jaccard和Bray–Curtis度量评估Beta多样性（组间多样性）。使用999次排列的ANOSIM检验评估分组显著性。在ASV水平上使用Observed、Pielou和Faith指数评估Alpha多样性（样本内多样性）。
2.3. 统计方法实现
研究中使用的所有标准统计检验，包括独立t检验、Mann–Whitney U检验、Brunner-Munzel检验和Mood’s中位数检验，均使用SciPy v1.15.1进行。通过Benjamini–Hochberg程序和Holm–Bonferroni方法进行的多重比较调整使用statsmodels v0.14.4实现。Cliff’s delta（δ）效应大小直接从Mann–Whitney U统计量计算。单样本和两样本Hodges–Lehmann估计量使用NumPy v2.0.1和SciPy v1.15.1包实现。为了进行比较分析，使用了R v4.4.2中ANCOMBC v2.8.1包的ANCOM-BC2算法。重现本工作中提出的分析所需的所有数据和脚本可在<
https://github.com/VeaLi/hle-ci-for-microbial-daa >获取（访问于2025年12月10日）。
3. 结果
3.1. 负二项数据中的经验假阳性率
我们进行了一项模拟实验，以估计几种统计方法的假阳性率（FPR）：独立t检验（TT）；带有Welch校正的t检验（TTW）；应用于带有伪计数的对数归一化数据的t检验（PL-TT和PL-TTw）；Mann–Whitney U检验（MWU）；Brunner-Munzel检验（BMT），一种不假设相同形状的MWU替代方法；Mood’s中位数检验（MMT），一种比MWU更稳健的替代方法；Cliff’s Delta（δ），一种源自MWU统计量的效应大小度量，代表不重叠程度；以及基于两样本Hodges–Lehmann估计量（HLE）的检验（HLE-CI，使用非重叠CI进行显著性估计）。
我们从具有不同位置参数的NB分布中生成了100,000个随机样本。对于每个样本，我们将数据随机分成相等的两半。由于这些半部分是同一大样本的随机子集，它们之间不应存在真实差异。然后，我们应用每种统计检验来比较这两半，以确定是否有任何方法错误地检测到显著差异。此过程允许我们估计每种方法应用于NB分布数据时的经验假阳性率。该分析的结果如表1所示。
在显著性水平α = 0.05下，所有基于t检验的方法（TT, TTW, PL-TT, PL-TTW）、MWU和BMT在比较来自同一NB分布的 subsamples 时，显示出比HLE-CI方法更高的假阳性率（FPR）。具体来说，这些检验对于大小为15的 subsamples 产生假阳性的可能性大约是HLE-CI的2.5倍，对于大小为25的是3.2倍，对于大小为50的是4.5倍，对于大小为100的是5.6倍。MWU的表现仅略优于t检验，而BMT的表现略差。MMT显示出比其他检验低1.5倍的FPR，但仍劣于HLE-CI。值得注意的是，虽然MMT在样本量为15时产生的FPR与HLE-CI相似，但对于大小为25的样本，产生假阳性的可能性是HLE-CI的1.9倍，对于大小为50的是2.9倍，对于大小为100的是4.4倍。
随着样本量的增加，t检验、MWU、BMT和MMT的FPR也增加。这种模式与其增长的统计功效一致，这使得能够检测到越来越微不足道、生物学上无关的差异。相比之下，HLE-CI的FPR随着样本量的增大而减小，表明中位数差异估计的准确性相应提高。在α = 0.01和α = 0.001时观察到类似但更明显的趋势。
无条件按Cliff’s δ效应大小（小和中等效应）过滤导致最高的FPR，尽管该速率随着样本量的增加而减小。将Cliff’s δ与MWU p值耦合限制了这些虚假发现；然而，所得性能并不比单独使用MWU而不配合相应效应大小有实质性改善。当使用大效应大小阈值进行过滤时，基于Cliff’s δ的过滤性能与HLE-CI方法相当。
该实验的结果表明：（1）HLE-CI方法产生的结果与MWU或相关检验不相似或相同（与t检验及其互补CI的关系相反）；（2）HLE-CI在从NB分布抽取的样本中产生最低的FPR；（3）HLE-CI即使在小样本量（n = 15–25）下也表现出最优的性能。
3.2. 负二项数据中的经验假阴性率
我们进行了进一步的模拟实验，以比较表现最佳的方法：MWU和MMT（带和不带多重比较校正（FDR-BH，表示为q））、基于Cliff’s δ的过滤和HLE-CI方法在旨在检测少量真实阳性标记的场景中的性能。与第一个实验类似，我们从NB分布中生成了1000对样本来模拟两个实验组的分类特征丰度。为了创建一个已知、小部分真实阳性的数据集，我们通过在一个组中将10个随机选择特征的值移动固定量（15%、25%或50%）来引入显著差异。这产生了一个包含10个真实标记的1000个特征的最终数据集。我们将模拟重复了100次，并计算了F1分数，该分数总结了精确度和召回率之间的平衡。F1分数是精确度和召回率的调和平均值，计算公式如方程（1）所示。结果如表2、表3和表4所示。
在15%的标记偏移下，具有95%和99% CI的HLE检验在所有数据集大小（表2）中表现出最佳性能，在所有α水平上 consistently 优于MWU和MMT。使用多重比较调整导致次优性能，在标记中心差异为15%时，有效地禁止了在具有少量真实阳性的数据集中进行标记检测。只有当每组的样本量增加到100时，按校正后p值过滤才显示出任何检测能力，产生比按原始MWU p值过滤更好的结果。然而，最高的F1分数仍然由HLE-CI过滤实现。
对于最大样本量（n = 100），按原始p值、校正后p值和HLE-CI结果过滤的最大F1分数分别为32.65%、18.2%和50.0%。因此，HLE检验相对于校正后p值提供了31.8%的F1分数相对改进。使用条件性Cliff’s δ实现的最高F1分数为19.65%，发生在中等效应大小阈值处。此外，即使样本量较小（n = 15–25），基于HLE的过滤也表现最佳。最后，MMT的表现并不优于MWU，这表明MMT较低的经验假阳性率（在实验1中观察到；表1）是其较低统计功效的结果，而不是更高的准确性。
类似地，在25%的标记偏移下，具有95%和99% CI的HLE过滤在所有数据集大小（表3）中表现出最佳性能，在所有α水平上 consistently 优于MWU和MMT。使用多重比较调整再次导致次优性能。只有当每组的样本量增加到50时，按校正后p值过滤才显示出任何检测能力。按原始MWU p值过滤通常由于更高的召回率而比按校正后p值过滤表现更好。最高的F1分数仍然由HLE-CI过滤实现。
对于最大样本量（n = 100），按原始p值、校正后p值和HLE检验过滤的最大F1分数分别为57.1%、46.2%和70.6%。因此，HLE检验相对于校正后p值提供了24.4%的F1分数相对改进。使用条件性Cliff’s δ实现的最高F1分数为46.2%，发生在中等效应大小阈值处。对于较小的样本量（n = 15–25），基于HLE的过滤再次在方法中表现出最佳性能。MMT的表现并不优于MWU。
最后，在50%的标记偏移下，具有99% CI的HLE检验在所有数据集大小（表4）中表现出最佳性能，在所有α水平上 consistently 优于MWU和MMT。使用多重比较调整导致较小样本量（n = 15）的次优性能，过度惩罚了标记检测。按校正后MWU p值过滤通常比按原始p值过滤表现更好（通过提高精度）。最高的F1分数仍然由HLE-CI过滤实现。
对于最大样本量（n = 100），按原始p值、校正后p值和HLE检验过滤的最大F1分数分别为84.2%、84.2%和87.3%。因此，HLE检验在检测大尺度变化（即50%偏移）的任务中与校正后p值表现相似。使用条件性Cliff’s δ实现的最高F1分数为82.4%，发生在中等效应大小阈值处。再次，基于HLE的过滤对于较小样本量（n = 15–25）显示出比任何其他方法更好的结果。
该实验的结果表明：（1）多重测试校正倾向于在真实阳性比例较低的数据集（如微生物组研究中常见的数据集）中过度惩罚比较；（2）HLE-CI方法在所有测试的样本量（n = 15–100）和效应大小（15–50%标记偏移）上 consistently 实现最佳或至少相当的F1分数。（3）基于HLE-CI的过滤在95–99%置信水平下表现最优，并且（4）可以作为微生物DAA中假设过滤的强大工具，与多重比较校正表现相当。
3.3. 人类肠道微生物组数据中假设过滤的精度和灵敏度
在上一节中，我们证明了在NB数据中可能偶然发现虚假的显著关联，并且多重比较校正会过度惩罚在真实阳性比例较低的数据集中的发现。在这里，我们评估了我们模拟中表现最佳的假设过滤方法HLE-CI在真实世界炎症性肠病（IBD）数据上的性能，比较IBD（n = 135）患者与其健康亲属（HR, n = 135）的16S微生物组谱。我们将这种基于HLE-CI的方法与ANCOM-BC2（无校正、带默认Holm–Bonferroni方法和带Benjamini–Hochberg程序）进行比较，ANCOM-BC2是一种专为微生物组数据设计的最先进方法。ANCOM-BC2基于具有测序深度校正的回归框架，并利用多重比较校正进行假设过滤。关于ANCOM-BC2的更多细节可以在[5]中找到。DAA仅考虑在完整数据集中普遍性至少为50%的特征（数据集中所有2272个已分类分类特征中的343个）。
3.3.1. 初步生物多样性分析
IBD是一个重大的全球负担，影响数百万人，并与生活质量下降和结直肠癌等并发症风险增加相关。如图1所示，IBD患者的肠道微生物组结构与健康对照组有显著差异。图1A显示了在ASV水平上基于Jaccard和Bray–Curtis距离的IBD和HR组的Beta多样性分析。IBD和HR样本形成显著分离的簇（ANOSIM, R = 0.13, p < 0.001 和 R = 0.11, p < 0.001），表明不同的微生物组谱，并提示IBD患者中存在菌群失调转变。IBD也以慢性肠道炎症和独特分类群的急剧减少为特征。在这个IBD患者数据集中，观察到ASV水平的Alpha多样性减少了2倍（图1B, Observed, p = 0.001 和 Pielou, p = 0.001）。无需任何特征选择的留一法交叉验证也显示了在属水平上的稳健可分离性（图1C, AUC = 0.84）。
总体而言，多样性分析表明IBD患者和HR对照的微生物谱存在显著差异，提示存在可靠的分类标记。这一点以及IBD数据研究深入的性质，使得该数据集以及一般的IBD数据成为微生物组DAA方法基准测试的热门选择。
3.3.2. 差异分析
使用HLE-CI过滤对稀释后的、TSS归一化的相对丰度数据进行DAA，并使用ANCOM-BC2对未稀释的、未归一化的计数进行DAA。ANCOM-BC2使用其默认设置运行。ANCOM-BC2的作者推荐使用Holm–Bonferroni方法而不是Benjamini–Hochberg程序进行多重比较校正。在我们的实验中，我们应用了推荐的Holm–Bonferroni方法和流行的Benjamini–Hoc h b e r g
