导读   

尽管通过宏基因组学发现的未培养病毒序列数量不断增加，而且它们对健康和疾病具有明显的重要性，但人类肠道病毒种群及其与肠道细菌的相互作用尚未得到很好的了解。一部分原因在于全病毒组数据集的缺乏和当前宏基因组数据中识别病毒序列的方法学限制。本研究结合基于深度学习的宏基因组分箱算法与配对的宏基因组和宏病毒组数据集，开发了Phages from Metagenomics Binning(PHAMB)，这种方法可以直接从大量宏基因组数据中对数千个病毒基因组进行分箱，同时将病毒基因组聚类为准确的分类病毒种群。当应用于人类微生物组计划2(HMP2)数据集时，PHAMB从1024个病毒种群中恢复了6077个高质量基因组，并确定了病毒-微生物宿主相互作用。PHAMB可以广泛地应用于现有和未来的宏基因组，有助于阐明病毒与其他微生物组成分的生态动力学。    

论文ID

原名：Genome binning of viral entities from bulk metagenomics data

译名：从大量宏基因组数据中对病毒进行基因组分箱

期刊：Nature Communications

IF：17.694

发表时间：2022.2

通讯作者：Simon Rasmussen

通讯作者单位：丹麦哥本哈根大学

DOI号：10.1038/s41467-022-28581-5

实验设计

结果

1 从宏基因组学数据中进行病毒种群分箱和组装的方法构架

我们使用VAMB来产生宏基因组bins，它的优点在于既能对微生物基因组进行分类，又能将各样本中的bins归入亚种或同种群。已有研究证明这种方法对于研究细菌和古细菌微生物组很有用，但它在病毒学研究中更有潜力，因为病毒的保守性差得多，更多样，由于没有像在细菌中发现的通用遗传标记，因此更难识别。同种特异性病毒基因组可用于直接识别和溯源同一样本队列中的病毒种群(图1a)。为了开发方法构架，我们采用两个基于Illumina鸟枪法测序的有成对的宏基因组和宏病毒组可用的数据集，其中2010年哥本哈根儿童哮喘前瞻性研究(COPSAC)数据集包括662个配对样本，Diabimmune队列数据集包含112个配对样本。这两个数据集都包含一个完成注释的病毒种类列表，分别为10021和328种病毒，我们将其作为训练和测试我们的方法的黄金标准。与COPSAC相比，Diabimmune宏病毒组的病毒富集程度较低(附图1)，因此，我们利用CheckV的平均氨基酸同源性(AAI)模型将宏病毒组中的基因组分为不同质量等级，即完整、高质量(HQ)、中等质量(MQ)、低质量(LQ)和未确定(ND)，以建立一个可比较的病毒数据集。

图1 从宏基因组数据中进行病毒种群分箱和组装的方法构架。

a 从分箱的宏基因组中探索病毒的工作流程图。首先，在分箱的宏基因组上训练RF模型；使用参考数据库确定细菌bins，使用配对的宏基因组组装病毒序列确定病毒。将病毒和细菌标记的bins用作训练和评估RF模型的输入。来自任何宏基因组(如人类肠道、土壤或海洋)的bins都可以通过RF模型进行注释，以提取病毒bins序列，并使用专用工具(如CheckV)进一步验证HQ病毒。此外，可以在纵向数据集中探索宿主-病毒动态信息，以建立温和噬菌体和病毒对宿主泛基因组的贡献。b 对来自Diabimmune的病毒bins的预测结果计算AUC、F1分数和Matthews相关性。这些性能分数是根据训练过的RF模型的概率分数和各种病毒预测工具的病毒bins分数汇总计算而来的。对于除RF模型外的所有工具，如果所有重叠群的病毒评分(以平均值、中位数或contig长度加权平均值计算)高于阈值，则基因组被标记为病毒。

viralVerify、Seeker、Virsorter2、Virfinder和DeepVirfinder使用的阈值分别为7、0.5、0.9、0.9、0.9。c 在Diabimmune或COPSAC数据集中以三种不同的完整性水平计算的从大量宏基因组中恢复的病毒基因组数量，评估为宏基因组的单个重叠群或病毒bins。用CheckV评估基因组的完整性，MQ≥50%，HQ≥90%，完整基因组定义为基于直接末端重复(DTR)或倒置末端重复的闭环基因组。d 相对于单重叠群评估，使用我们的方法在Diabimmune或COPSAC队列中发现的病毒基因组的百分比增加。e 与(b)类似，对训练过的RF模型和各种病毒预测因子计算预测性能得分，但基于CAMI模拟病毒基因组(包括细菌、病毒和质粒)的预测结果。 

2 病毒分箱技术比单个重叠群层面的方法更有效果 

宏基因组样本的分箱结果可能包含几十万个bins，因此我们首先开发了一个随机森林(RF)模型来区分病毒类基因组bins和细菌类基因组bins。RF模型利用分箱结果的聚类信息，并在样品特定的分箱结果中聚集信息，形成亚种聚类。在这一研究中，我们发现RF模型能够非常有效地区分细菌和病毒种群，在验证集上的曲线下面积(AUC)为0.99，Matthews相关系数(MCC)为0.91(图1b和附表1)。与单重叠群评估方法相比，RF模型的优势在于其他方法的AUC最高为0.86，MCC最高为0.16。

这种性能上的差异可能是由RF模型在bin水平上的评估所解释的，其中一个病毒得分较低的序列不会导致对整个bin的错误预测。例如，与根据CheckV的单重叠群评估相比，我们在Diabimmune和COPSAC数据集中恢复的HQ bins分别增加了200(190%)和771(95%)(图1c，d)。基于CheckV的单重叠群评估，我们发现在COPSAC和Diabmmune数据集中，分别有97.7%和95.3%的HQ contigs被归入HQ bins，这意味着有一小部分的HQ contigs(高达2.3和4.7%)在分箱过程中丢失，造成基因组恢复的净增加降低，但这种影响可以通过单重叠群评估来恢复。最后，我们观察到，当使用两个数据集中的病毒bins时，每个病毒的病毒标志基因数量明显增多(T检验，双侧，t = 16.85，P < 0.0005)，而病毒基因组长度和病毒比例基本相当(附图2)。 

3 病毒分箱技术在模拟病毒组数据中表现出色 

然后，我们使用模拟数据集(包括两个纯病毒数据集和一个包含细菌、质粒和病毒的混合数据集)研究了VAMB的病毒分箱(binning)性能和预测性能。两个纯病毒数据集由80个crAss样病毒和50个从MGV数据库中随机抽取的小基因组(0.5)和不断提高的准确率水平下，复原bins的数量越来越多。对于质粒，这两个工具效率相当，以>0.5的召回率和>0.95的准确率对10/20个质粒进行了分箱(附图3b)。接下来，我们研究了病毒基因组大小和高度相似的病毒如何影响分箱性能。

为此，我们对较小的病毒基因组(0.9)捕获了大多数小基因组病毒，这证实了小病毒易于分箱的特点(附图4c)，表明基因组大小对分选性能的影响较小。最后，为了进一步验证RF模型，我们比较了单个重叠群病毒预测工具在预测bin是病毒还是细菌方面的效果(图1e)。使用混合模拟数据集，与RF模型相比，单重叠群方法鉴别性能明显更低。例如，具有高AUC(高达0.98)的多个单株病毒预测工具显示出较低的MCC得分，这意味着在给定的阈值下预测不是非常准确(图1e和附图5，6)。然后，我们尝试优化每个单重叠群病毒预测工具的决策阈值(附图5，6)，这略微提高了MCC得分。例如，viralVerify在模拟数据上的AUC为0.98，表明它能有效地区分细菌和病毒基因组，但在细菌和病毒得分分布上存在重叠。因此，即使有一个优化的阈值，viralVerify显示的MCC为0.39。相比之下，RF模型的AUC(0.93)和MCC(0.87)均较高。因此，我们发现RF模型，其次是viralVerify，是混合微生物组装数据集中最适合的分箱方法。虽然RF模型将质粒错误地预测为病毒，但我们发现，在下游使用CheckV有助于做出最终的评估，因为质粒bin包含多个细菌起源的基因，通常被归类为 "NA "或被不太精确的HMM模型挑出(附图7)。 

4 宏基因组分箱可以识别宏病毒组无法鉴定出来的病毒基因组

当用我们的VAMB和RF模型方法进行分箱时，我们在COPSAC和Diabimmune数据集中分别获得了4480和916个MQ或HQ代表性分箱的病毒bins。然后，我们将所有VAMB簇视为"病毒种群"，从而得到2428和534个至少具有1个MQ或更高的病毒bins的病毒群。在将从宏基因组数据集获得的病毒群与相应的宏病毒组进行比较后，我们分别在物种(ANI>95)水平和菌株(ANI>97)水平(图2a)上恢复了在宏病毒组中建立的17-36%和9-28%的HQ病毒(对应于527和2676个宏病毒组病毒群)。

在宏基因组中恢复的病毒比例大大高于近期研究中的估计值(8.5-10%)。这很有意思，因为深度测序的宏基因组可能会捕获宏基因组中通常没有发现的多种低丰度病毒。此外，我们发现46-69%的HQ宏基因组病毒种群(对应于Diabimmune中的124个和COPSAC中的839个病毒种群)在宏基因组中没有发现，这表明很大一部分病毒组可能在病毒富集过程中丢失或者在诱导形式中没有体现，因为它们是整合的原噬菌体(图2b)。然而，我们也发现，在宏病毒组中有65-83%的HQ病毒种群未在宏基因组数据中发现(Diabimmune中共有197个，COPSAC中共有2589个)，这表明反过来也是如此。对于在COPSAC大量和宏病毒组中发现的病毒子集，我们估计病毒bins的平均完整性更高(T检验，双侧，T=34.02，CI=24.4；27.4，P=2.2e-16)(图2c)。总体而言，我们发现很大一部分肠道病毒种群可以从宏基因组数据中恢复，并且与宏病毒组数据相比，检索的完整性更高。

图2 宏基因组分箱可以识别宏病毒组无法鉴定出来的病毒基因组。

a COPSAC和Diabimmune宏病毒组中的不同完整度病毒的比例，或全部用CheckV确定，这些病毒在同一队列的大量宏基因组的VAMB bins中确定。如果对齐部分至少为75%，且ANI>90、>95或>97.5，则我们将宏病毒组病毒定义为基于FastANI的VAMB bins的复原病毒。b使用CheckV确定的不同完整性水平上的病毒种群百分比，在宏病毒组(MVX)和宏基因组(MGX)中均可识别，或仅为其中一个数据集所独有。共享种群的最小序列覆盖率为75%，ANI大于95%。

(1) MVX中的MGX：在MGX中发现的病毒种群在MVX中也发现的百分比。

(2) MGX不在MVX中：MGX特有的病毒种群的百分比，即在MVX中没有发现。

(3) MGX中的MVX：在MVX中发现的病毒种群也在MGX中发现的百分比。

(4) MVX不在MGX中：MVX特有的病毒种群的百分比，即在MGX中没有发现。

c 病毒基因组的完整性是基于宏病毒组和宏基因组中发现的n = 2646种病毒而估计的，这些病毒在CheckV数据库中具有相同的最接近的参考序列。d 来自大量宏基因组的病毒bins中与宏病毒组中最接近的病毒参考序列不一致的contig数。在大多数病毒bins中，所有contigs都与最接近的参考序列对齐。ANI，平均核苷酸同源性。 

5 病毒分箱污染低 

最后，我们想研究技术上的“分箱错误”和污染性重叠群的情况，因为这可能会夸大病毒基因组的大小，影响质量评估和下游分析。基于COPSAC数据集中与宏病毒组病毒高度相似的病毒bins(n = 1705)(见方法)，我们发现在91.4%的情况下，每个bin都不包含不相关的重叠群(图2d)。仅考虑多重叠群bins时(n=570)，我们计算出碱基对的平均bin纯度为97.4%，这意味着平均有2.55%的基因组没有与相应的MVX病毒对齐，这表明存在污染，或者说，在大量宏基因组数据集中有更完整的病毒。我们根据模拟数据进一步分析了污染程度，其中87.6%的病毒bins的准确率为1(附图8a)。对于多重叠群bins，我们计算出平均bin纯度为94.5%(中位数100%)，支持了真实数据的结果，即大多数bins的污染程度较低。总之，我们结合分箱和机器学习的方法提高了从宏基因组数据中识别和恢复病毒基因组的能力，并概述了直接从污染程度较低的人类肠道微生物组样本中分箱片段和完整病毒的可能性。

6 重组HMP2 IBD肠道宏基因组队列的病毒组

然后我们将本研究的方法应用于HMP2 IBD队列(包括27名健康对照、65名CD和38名UC患者)。这些样本是以纵向方式收集，每个患者有1-26个样本。重要的是，该队列中没有已有的宏病毒组学数据，使用我们的方法，我们能够鉴定队列中的细菌和病毒种群，并仅使用宏基因组学数据探索它们在IBD中的动态变化。从该队列中，我们恢复了577个完整、6077个HQ、9704个MQ(图3a)和122,107个LQ病毒bins，对应于263个完整、1024个HQ、2238个MQ和44,017个LQ病毒种群。与单重叠群评估相比，我们还发现基因组大小>200 kbp的较大病毒/巨型病毒的基因组完整性有所增加(附图9)。在所有的数据集中，我们观察到54个分箱的假定巨型病毒(补充数据1)。此外，我们还观察到，作为单个重叠群和病毒bins恢复的病毒具有相似的病毒长度分布模式，两者都与CheckV质量等级相关(图3b)。

图3 重组HMP2 IBD肠道宏基因组队列的病毒组。

a HMP2中三种不同完整程度的病毒基因组的数量，以单个重叠群或来自宏基因组的病毒bins进行评估。基因组完整性的评估是用CheckV确定的，这里显示的是中等质量≥50%(MQ)，高质量≥90%(HQ)，完整=基于直接末端重复或倒置末端重复的封闭基因组。b HMP2中以来自宏基因组的单重叠群(n = 215009)或病毒bins(n = 138367)评估的四个不同完整程度的病毒基因组的序列长度分布(kbp)。c 基于同一VAMB簇内各bins之间的成对ANI基因组测定的ANI中值。ANI中值在0-25个bins的小型VAMB簇和300-400个bins的大型VAMB簇中均高于97.5。d 基于TerL的crAss样bins的无根进化树的Cladogram图。五个不同的VAMB簇进行不同颜色着色，说明了高度的单系关系。ANI，平均核苷酸同源性%；DTR，直接末端重复；ITR，倒置末端重复；Kbp，千碱基对。 

7 病毒种群分类高度一致 

然后，我们分析并发现我们的病毒种群的分类一致性很高，从MQ到完整病毒种群的簇内平均核苷酸同源性(ANI)的中位数为97.3-99.3%(附图11)。即使在有超过100个样本特异性病毒bins的病毒簇中，簇内ANI中值也始终很高(中位数=97.1-98.5%)(图3c)。但是，簇间ANI较低，在91.7-92.8%之间，更接近于属水平。因此，我们的方法能够在不同的样本中识别和聚类接近菌株水平的病毒基因组。例如，在HMP2数据集中，我们确定了50个不同的病毒种群，总共有916个MQ或更好的crAss样病毒bins。此时，病毒种群653对应于原型crAssphage噬菌体，占HMP2数据集中发现的916个crAss样基因组中的253个。然后，我们用所有这916个bins来构建一个基于TerL的系统发育树，发现根据其种群的病毒基因组的位置高度一致(图3d和附图12)。病毒种群653形成了一个单系分支，而所有其他crAss样bins中除一个bin外均为单系。因此，将crAss样基因组划分为分箱簇很可能代表了实际的病毒多样性。综上所述，这表明我们的无参分箱产生了分类准确的病毒种群，从而聚集了各样本之间高度相似的病毒基因组。 

8 健康受试者的宏基因组病毒组是个性化且高度稳定的

一些宏病毒组研究报告了人类肠道中存在大量稳定的病毒。我们发现，在HMP2队列中，非IBD受试者的肠道病毒群是高度个性化且稳定的，这体现在非IBD受试者的样本与UC(T检验，双侧P = 0.017，t = -2.47，CI = -0.01;-0.13)和CD受试者(T检验，双侧，P = 0.023，t = -2.3，CI = -0.12;-0.01)相比，Bray-Curtis差异度更低(图4a，b)。此外，Price等人所定义的菌群失调样本可以用主成分分析(PCoA)明确分开，其中病毒组解释了4.2和3.4%的变化(图4c)。这一点通过对病毒(P < 10 - 3，R2 = 1.6%，F = 9.51，permutations = 999)和细菌丰度谱(P < 10-3，R2 = 3.0%，F = 11.97)的PERMANOVA分析得到证实，并表现出影响病毒组和细菌组的失调现象。α-多样性指标支持这一点，因为与UC和CD受试者相比，非IBD受试者的Shannon-多样性(SD)更高(T检验，双侧，P = 0.000155，t = -3.79和P = 7.9e-09，t = -5.81)，而菌群失调影响了每个患者组，导致SD明显降低。

与此相应，病毒丰富度在UC(双侧T检验，P=1.44e-15，t=-8.09，CI=-12.40;-19.80)和CD(双侧T检验，P=