基于16S rDNA-seq计数数据的微生物互作网络推断方法综合基准评估：从共现网络到生态关联

《PLOS Computational Biology》：From 16S rDNA-seq count data to microbial interactions: A comprehensive benchmark of co-occurrence network inference methods

　　

摘要专业翻译 微生物群落内生态相互作用的研究对于深入理解群落的组装与进化至关重要，能够实现对群落功能的定向改造。目前已发展出多种推断方法以获取生态关联。本研究聚焦于共现计算方法，该类方法通过测量不同采样群体间物种计数的统计依赖性，从宏基因组丰度数据中推断关联。基于培养实验的局限性阻碍了评估这些方法有效性的途径，导致其在宏基因组学研究中的应用十分有限。为全面评估网络推断方法的性能，研究人员在广泛的实验场景下对18种不同的计算方法进行了基准测试。研究人员使用由已知底层相互作用网络生成的合成16S rDNA-seq计数数据，从而实现系统且受控的评估。该基准测试涵盖了物种数量与样本数量、样本丰富度、测序深度以及相互作用网络结构等不同场景。研究人员从推断网络的准确性（使用二元分类理论中的多种指标以及网络科学中的指标）和计算负担两方面评估了各方法的性能。研究发现，各方法仅在受限且往往不切实际的应用领域中表现出有效结果，而在更真实的模拟参数下其有效性显著下降。研究结果凸显了基于共现方法的局限性，以及微生物生态关联推断领域进一步发展的必要性。

论文解读文章

一、研究背景与意义

微生物群落广泛存在于土壤、植物、海洋及人体等各类生态系统中，对宿主健康或生态系统平衡具有重要影响。人类共生微生物群落的失调与心血管疾病、炎症性肠病、阴道病、HIV、癌症及糖尿病等多种疾病密切相关。随着下一代测序技术（Next Generation Sequencing, NGS）特别是16S rRNA基因靶向测序（16S rDNA-seq）的发展，微生物群落研究已突破传统培养实验的局限。在众多分析方法中，从物种共现数据推断微生物互作网络引起了广泛关注，该方法仅依赖横断面微生物丰度谱即可生成生态学假说，具有数据获取成本低、技术成熟等优势。

然而，从16S rDNA-seq数据推断微生物互作面临诸多挑战，包括维度灾难、成分性（compositionality）、批次效应以及缺失数据/零膨胀等问题。已有多种基于共现的推断方法被开发，从简单的Pearson相关、Spearman相关等基础度量，到SparCC、CCLasso等成分感知方法，再到SPIEC-EASI、SPRING、gCoda、FlashWeave等条件依赖方法，以及HARMONIES、COZINE、SECOM等针对零膨胀和测序偏差建模的新方法。尽管方法众多，但共现方法在宏基因组研究中的应用仍然有限，其结果常受到质疑。根本原因在于缺乏大规模生物学金标准来验证推断结果，而基于培养实验的验证途径又受到限制。因此，利用已知底层互作网络的模拟数据进行系统评估成为必要选择。

本研究的主要贡献包括：纳入了迄今最大规模的方法集合（共18种推断技术）；涵盖广泛的模拟场景，包括样本量、网络结构、物种丰富度、文库大小、互作类型和样本多样性等多个因素；使用多种评估指标，包括二元分类指标（Precision、Recall、MCC）以及网络科学指标（图拓扑指数和中心性度量）；采用新开发的N2SIMBA模拟器，该模拟器不依赖成对互作模型，减少了生成数据中的不切实际的过度简化；并首次评估了计算负担这一此前研究未涉及的方面。

二、关键技术与方法

研究人员使用N2SIMBA模拟器生成合成16S rDNA-seq数据，该模拟器基于社区模拟器（Community Simulator）实现微生物群落消费者-资源模型（Microbial Consumer-Resource Model, MiCRM），通过非线性方程组模拟群落动态。模拟流程包括：根据用户定义的互作网络生成C和D矩阵编码物种-资源关系；从真实测序数据预设中随机选取初始丰度；模拟至稳态后通过扰动资源产生不同β多样性的样本；最后采用多元超几何模型（Multivariate Hypergeometric Model）采样生成16S rDNA-seq计数矩阵。研究共生成189种参数配置、9450个模拟计数表。测试方法分为无条件依赖方法和条件依赖方法两大类，共18种推断技术，包括5种基础度量（Pearson相关、Spearman相关、Bray-Curtis相异度、偏Pearson相关、偏Spearman相关）和13种专用推断工具。评估指标涵盖边层面（MCC等）、节点层面（介数中心性、特征向量中心性、紧密中心性）和图层面（度分布KS统计量、密度、直径、半径、平均距离、平均聚类系数）三类指标，并测量了各方法的运行时间和峰值内存。

三、研究结果

图生成模型和样本多样性对性能产生一致且显著的影响

在默认模拟参数下，研究人员考察了3种图生成模型（ER、BA、WS）和3个β多样性水平（低、中、高）共9种场景。结果显示，方法整体表现中等，中位MCC值低于0.5，仅在WS图模型下部分方法取得较高值（SPIEC-EASI（MB）MCC≈0.7、SPIEC-EASI（Glasso）MCC≈0.68、SPRING MCC≈0.65）。主要限制因素在于方法通常推断出较小的边子集，导致大量假阴性，而已调用的边通常是正确的。条件依赖方法整体表现最佳，SPRING几乎始终取得最高MCC值，CCLasso是唯一能与条件依赖方法表现相当的非条件方法。WS图模型对方法最有利，其次是BA和ER，这可能源于WS模型的规则拓扑结构使度分布近乎均匀。β多样性水平越高，方法性能越好，表明在资源丰富的环境（如人类肠道）中推断互作比资源贫乏环境（如土壤）更困难。

增加样本量可改善网络推断性能

研究人员模拟了50个物种的群落，样本量从50到1000变化。结果表明，方法性能随样本量增加而提升。当n=50=p时所有方法表现不佳（MCC<0.2）；当n=1000时部分方法可达MCC≈0.6，WS模型下MCC≈0.7；大多数方法在n>300后进入饱和状态。值得注意的是，偏Pearson相关在n>500时成为所有方法中表现最优者，而Pearson和Spearman相关的结果与其他方法相当但略逊。在节点层面，样本量增加提升了推断图与真实图的介数中心性R2值，其中偏Pearson相关（R2=0.61）、CCLasso（R2=0.59）、Bray-Curtis相异度（R2=0.46）和SPRING（R2=0.43）提升最为显著。在图层面，KS统计量随样本量增加而降低，SPRING取得最佳结果（KS≈0.2）。在计算负担方面，n=500时gCoda和FlashWeave（两种实现）运行时间最短（约2秒），HARMONIES最慢（约100分钟）；内存方面CCLasso最优（197MB），HARMONIES和SPIEC-EASI（GL）最差（分别为1.8GB和1.5GB）。

物种数量增加导致性能下降

研究人员模拟了物种数p=25至p=300的场景。结果显示，随着物种数增加，所有方法的MCC值普遍下降，这反映了高维数据带来的推断困难。在p=25时部分方法表现较好，但当p达到300时，多数方法的MCC值降至0.2以下。条件依赖方法在高维场景下相对更具优势，但整体性能仍然有限。

网络密度对性能的影响

研究人员考察了网络密度从2%到15%的变化。结果表明，低密度网络（2%）下方法表现最佳，这与多数推断方法假设网络稀疏的前提一致。随着密度增加至15%，所有方法的性能均显著下降，说明网络密度是影响推断效果的关键因素之一。

负互作比例的影响

研究人员变化了竞争性互作比例（负互作比例）从0%到100%。结果显示，负互作比例对方法性能的影响相对较小，但总体趋势是负互作比例适中时性能略优，极端比例（0%或100%）时性能略有下降。

文库大小的影响

研究人员通过缩放因子K（从100%降至10%）模拟不同文库大小。结果表明，文库大小减小导致计数矩阵更稀疏、成分效应更强，所有方法的性能均随之下降。条件依赖方法对文库大小的变化相对更敏感。

计算负担评估

研究人员测量了各方法在代表性场景下的运行时间和峰值内存。结果显示，基础方法和Proportionality的计算成本随样本量增加而显著上升；SECOM（NL）和HARMONIES在样本量增加时运行时间分别增加5340%和4500%。在默认参数下，FlashWeave（两种实现）和gCoda运行最快，HARMONIES最慢；内存方面CCLasso最优，HARMONIES和SPIEC-EASI（GL）最差。

四、讨论与结论

研究结果表明，共现网络推断方法仅在受限且往往不切实际的应用领域中表现有效，而在更真实的模拟参数下其有效性显著降低。条件依赖方法在边层面推断上整体优于无条件依赖方法，SPRING和CCLasso是综合表现最优的工具。然而，所有方法都存在大量假阴性的问题，即推断出的边通常正确但覆盖范围有限。样本量是影响性能的最重要可控因素，n>300后性能趋于饱和。网络拓扑结构（WS优于BA优于ER）和β多样性水平对性能有显著影响。计算负担方面，不同方法差异巨大，实际应用中需根据数据规模和资源条件选择合适的方法。

研究结论强调，基于共现的方法在微生物生态关联推断中存在固有局限，需要进一步发展更有效的推断方法。研究团队公开了合成数据集、结果文件和分析代码，以及包含所有测试方法的Singularity容器，以确保研究结果的可重复性，并为研究人员选择合适方法提供参考依据。