2023-04-27 14:00

衡量你的噬菌体:噬菌体鉴定工具在宏基因组测序数据的基准

摘要

背景

在宏基因组数据集中预测噬菌体序列已经成为一个相当感兴趣的话题,导致许多新的生物信息学工具的发展。对十种最先进的噬菌体鉴定工具进行了比较分析,以告知它们在微生物组研究中的使用情况。

方法

由RefSeq完整基因组生成的人工contigs代表噬菌体、质粒和染色体,以及先前测序的包含四种噬菌体的模拟群落,用于评估工具的精度、召回率和F1分数。我们还生成了一个随机洗牌序列的数据集,以量化假阳性呼叫。此外,还使用了一组先前模拟的病毒组来评估每种工具输出中的多样性偏差。

结果

在RefSeq人工配置数据集中,VIBRANT和VirSorter2获得了最高的F1分数(0.93),其他几个工具也表现良好。Kraken2在模拟社区基准中F1得分最高(0.86),比排名第二的DeepVirFinder高出0.3,这主要是因为它的精度很高(0.96)。一般来说,基于k-mer的工具比参考相似性工具和基于基因的方法表现更好。一些工具,最著名的是PPR-meta,在随机洗牌的序列中发现了大量的假阳性。当分析每个工具从病毒组预测的基因组多样性时,大多数工具产生的病毒基因组集具有与原始种群相似的α和β多样性模式,但Seeker是一个明显的例外。

结论

本研究提供了用于评估噬菌体检测工具性能的关键指标,为进一步比较其他病毒发现工具提供了框架,并讨论了使用这些工具的最佳策略。我们强调,在宏基因组数据集中选择噬菌体鉴定工具及其参数可能会对结果产生偏差,并为不同的用例场景提供指导。我们还提供了我们的基准数据集供下载,以便于将来对噬菌体鉴定工具进行比较。

视频摘要

介绍

噬菌体(噬菌体)和古细菌病毒在全球范围内普遍存在,种类繁多,并且在大多数生物群系中通常超过其原核宿主。噬菌体通过促进共同进化关系[2,3,4],必需营养素的生物地球化学循环[5,6,7],以及通过水平基因转移促进微生物进化[8,9,10],在微生物群落中发挥关键作用。尽管噬菌体对所有微生物生态系统都有丰富的影响,但它们仍然是复杂微生物群中研究和了解最少的成员之一。噬菌体是专性寄生虫,需要寄主的机器通过细胞裂解进行复制和传播。它们可以是溶性或温带噬菌体,前者只能遵循溶性生命周期,而温带噬菌体可以遵循溶性或溶原性生命周期。在裂解周期中,噬菌体劫持宿主细胞机制以产生新的病毒颗粒。在溶原循环中,噬菌体可以将它们的基因组作为线性DNA或自我复制的自主质粒整合到它们的细菌或古细菌宿主基因组染色体中。此外,第三个生命周期被称为假溶原性,在这个生命周期中,要么溶解噬菌体感染被停止,要么没有噬菌体形成[13]。

传统上,噬菌体鉴定和表征依赖于分离和培养技术,这是耗时的,往往需要大量的专业知识。这也是不切实际的,因为许多宿主及其噬菌体不能在实验室条件下培养。高通量新一代测序技术的到来使得来自不同环境的宏基因组数据得以常规生成。宏基因组测序允许直接鉴定和分析样本中的所有遗传物质,而不考虑可培养性。

宏基因组研究可以选择对整个群落宏基因组进行测序,然后计算分离病毒序列,或者在文库制备之前物理分离病毒片段以产生元病毒组。由于噬菌体与细胞部分的关联,后一种方法有消除大部分噬菌体的风险。这是由于噬菌体作为前噬菌体[16]被整合到宿主基因组中,附着在宿主表面[17]上,或者当它们处于假溶原状态时[18,19,20]。纯化方法也可以去除某些类型的噬菌体,例如氯仿可以灭活脂质包膜和/或丝状噬菌体[21,22],增加采样偏差。这一过程也导致DNA产率低,导致一些转移病毒组研究不得不使用多次位移扩增(MDA)来获得足够数量的DNA以生成文库[11]。MDA通过优先扩增小的环状ssDNA噬菌体,如来自微病毒科[25]的噬菌体,已被证明对病毒组成产生显著的偏倚[23,24]。尽管存在这些缺陷,但纯化步骤产生宿主污染很少的转病毒体,尽管很难产生不含任何细胞物质[26]的病毒部分。由于排除了大约100倍大的细菌和古细菌基因组,元病毒组还具有能够在相同测序深度下识别低丰度噬菌体的优势。或者,整个群落宏基因组测序可以同时提供对宿主和病毒组分的见解,允许分析宿主-噬菌体动力学。整合噬菌体,或已发现在某些环境中普遍存在的噬菌体,可以被识别,因为宿主基因组也在这一过程中测序。在这项研究中,我们专注于从整个群落宏基因组中计算提取噬菌体序列,因为与宿主相比,这些序列通常只占测序数据的一小部分。

在过去的5年中,已经开发了许多用于从混合元基因组和病毒组组装中识别病毒序列的工具(表1)。VirSorter[28]是其中的一个,之前的工具主要用于噬菌体预测(PhiSpy [29], Phage_Finder [30], PHAST/PHASTER [31], proinder[32])或病毒分析(metaVir2) [33], virome[34])。VirSorter通过检测与参考数据库具有同源性的病毒标志基因和基于不同指标(病毒样基因、Pfam基因、未表征基因、短基因和链转换)建立概率模型来识别噬菌体序列,这些指标测量每个预测的置信度。自VirSorter发布以来,其他基于基因同源的工具,如VIBRANT[35]和VirSorter2[36],已经被开发出来。基于隐马尔可夫模型(HMM)对多个数据库的蛋白质注释,使用多层感知器神经网络来恢复感染细菌和古细菌的多种噬菌体,包括整合的噬菌体。除此之外,它还具有鉴定后的辅助代谢基因和途径的特征。VirSorter2在其前身的基础上,将针对五种不同病毒群的五种不同的随机森林分类器合并到一个算法中,以提高其能够准确检测到的病毒的多样性。与上述工具不同,metaPhinder使用基于blast的同源性命中到自定义数据库来计算平均核苷酸身份和序列是病毒起源的可能性。

表1微生物生态系统中噬菌体序列鉴定和预测工具概述
全尺寸工作台

VirFinder是第一个利用k-mer签名的机器学习病毒识别工具。VirFinder的病毒序列回收率明显高于VirSorter,特别是在较短的序列(< 5 kbp)上,但在不同的环境中存在不同的性能问题,这可能是由于用于训练机器学习模型[50]的参考数据引入的偏差。DeepVirFinder[37]通过应用卷积神经网络对VirFinder进行改进,该神经网络是在包含来自环境元病毒组测序数据的病毒序列的扩大数据集上训练的。DeepVirFinder比其前身VirFinder在所有序列长度上都具有更高的病毒识别能力,同时减轻了后者对易于在实验室培养的噬菌体的偏见。Kraken2是基于k-mer的宏基因组分类器[51],可用于病毒检测[52]。它将k-mers查询到一个数据库,该数据库将其与它们最低的共同祖先分类群相关联,然后用于分配分类标签。PPR-meta使用三个卷积神经网络来识别序列是噬菌体,质粒还是染色体起源[42]。序列特征由网络直接提取,而不是使用预先选择的特征,如k-mer特征或基因。这三个网络也在三组不同的序列长度上进行训练,以提高其在较短片段上的性能,这是一些基于基因的工具由于可用于分析的全长基因数量较少而遇到的困难。Seeker也使用神经网络,在这种情况下是长短期记忆(LSTM)模型,它不是基于预先选择的特征[43]。metaviralSPAdes[40]使用了一种完全不同的方法,利用病毒和细菌染色体在组装图中的深度变化。该工具分为三个独立的模块:一个基于metaSPAdes的专门汇编器(viralAssembly);病毒识别模块,使用朴素贝叶斯分类器(viralVerify)对病毒/细菌/不确定的组合进行分类;以及计算构建的病毒序列与已知病毒的相似性的模块(viralComplete)。

值得注意的是,机器学习工具具有识别新物种的潜力,这对于理论上仍然未知的噬菌体的巨大多样性尤其重要。随着众多工具和方法的发展,需要一个全面的比较和基准来评估哪些工具最适用于研究人员。每种方法的性能可以根据样本内容、组装方法、序列长度、分类阈值和其他自定义参数而变化。为了解决这些问题,我们对十种宏基因组病毒鉴定工具进行了基准测试,包括人工配置、模拟社区和真实样本。

结果

用RefSeq噬菌体和no非病毒人工配置

选取10种常用的元基因组病毒序列鉴定工具进行评价:DeepVirFinder、Kraken2、metaPhinder、PPR-meta、Seeker、VIBRANT、viralVerify、VirFinder、VirSorter和VirSorter2。所有这些工具都可以在本地运行,而不依赖于web服务器,接受宏基因组配置作为输入,并且在过去十年中已经发布。

我们首先在相同的统一数据集上评估所有程序。下载2020年1月1日至2021年8月12日期间存放在RefSeq中的所有完整噬菌体基因组,对其进行质量控制,并将其均匀碎片化至1至15 kbp的大小,以创建一组真正阳性的人工组份。从同一时期提交的所有RefSeq细菌和古细菌染色体和质粒中构建阴性集。采取了多个步骤确保这些数据集不包括假阳性或假阴性。首先,所有细菌、古细菌和噬菌体的RefSeq基因组在2020年之前沉积,并使用每个机器学习工具的训练集来重复本研究中使用的数据集,以去除可能导致某些工具过拟合的任何相似序列。此外,删除开放阅读框中与pVOG数据库有≥30%匹配的染色体和质粒序列,以排除任何剩余的病毒序列。由于阴性数据集比阳性数据集大得多,我们对阴性数据集进行了14.3倍的次采样,得到253条宿主染色体和309个宿主质粒(表2)。选择这个采样率来产生噬菌体:宿主比例(~ 1:19),这与人类肠道微生物组[53]中发现的相似。最后,我们使用两种最先进的噬菌体检测工具Phigaro[54]和PhageBoost[55]从染色体和质粒序列中去除整合的噬菌体,以防止它们被错误地识别为病毒组合。共从染色体组中分离出2088个噬菌体,从质粒组中分离出91个。所得到的序列随后被分割成人工组合,并使用不同的工具进行分析(图1)。

表2 RefSeq基准测试工作流程中每个阶段的序列数量
图1
figure 1

RefSeq基准测试工作流程概述下载了2020年1月1日至2021年8月12日期间存放在RefSeq数据库中的所有细菌和古细菌染色体、质粒和噬菌体基因组。噬菌体基因组用于创建阳性测试集,染色体和质粒用于阴性测试集。使用基准测试的每个机器/深度学习工具的训练集(红色突出显示)以及2020年之前沉积的任何RefSeq序列对序列进行去重复。对阴性组进行取样,产生约1:19的阳性:阴性比例,以复制典型的肠道微生物组。用Phigaro和PhageBoost对噬菌体进行鉴定和去除。所有大于30%的开放读帧与原核病毒同源群数据库相匹配的宿主序列都被移除。然后将所有序列均匀地分割成长度在1至15 kbp之间的人工片段。然后在人工配置集上运行所有识别工具

除Kraken2外,所有评估的程序都产生了病毒鉴定的阈值或置信范围。对于分配连续阈值(分数,身份或概率)的工具(DeepVirFinder, metaPhinder, PPR meta, Seeker, VirFinder和VirSorter2),绘制F1曲线,并确定最佳阈值(附加文件1)。对于VirSorter和viralVerify,使用返回最高F1分数的类别(附加文件2)。在大多数工具中,在精度和召回率之间存在权衡。这可能是由于放宽了阈值,允许检测到更多的病毒和非病毒序列,增加了召回率,同时降低了精度。对于VIBRANT和VirSorter,阳性数据集分别在病毒组模式和病毒组去污模式下运行,因为这可以通过调整工具的敏感性来提高主要由病毒序列组成的样本中的病毒回收率[28,35]。我们在该数据集上进行基准测试的工具在F1得分(0.44-0.93)、精度(0.47-1.00)和召回率(0.46-0.96)方面具有高度可变的性能(图2)。VirSorter2、VIBRANT和PPR-meta的F1得分最高,分别为0.93、0.93和0.92。VirSorter2以高精度(0.92)和高召回率(0.93)实现了这一点,VIBRANT具有较高的精度(0.97)和较低的召回率(0.89),PPR-meta的精度略低,为0.88,但召回率较高,为0.96。其余大多数工具表现良好,其中六个工具(DeepVirFinder、Kraken2、metaPhinder、VirFinder和VirSorter)的F1分数超过0.83。Kraken2的精确分数接近1,只有2个染色体片段和没有标记为病毒的质粒片段,同时仍然正确识别超过5000个噬菌体片段。这些染色体片段原本是白黄链霉菌J1074/R2和火茶糖菌BEU9完整染色体的一部分,Kraken2鉴定它们分别为Streptomyces phiC31噬菌体和Sulfolobus virus 1。Kraken2鉴定的片段与参考文献的整体比对显示,它们之间的同源性有限,并表明这些片段是真阳性和假阳性的预测(额外文件3)。ViralVerify的精度得分较低,为0.55,而其召回得分为0.88,与其他工具相当。与其他工具相比,Seeker在准确率(0.48)和召回率(0.41)方面的表现都较差。一般来说,基于k-mer的工具比参考的基于相似性/基因的工具表现得更好,尽管所调查工具的样本量太小,无法得出具有统计学意义的结论。在我们的基准中,阳性集合中只有0.06%(4/6665)的噬菌体contigs未被任何工具检测到,89.7%被超过一半的工具检测到,11.6%(775/6665)被所有10个工具发现(附加文件4)。在噬菌体鉴定中,工具之间没有明显的分类偏差,除了Seeker分类的阿克曼病毒科片段比其他工具少(附加文件5)。

图2
figure 2

人工RefSeq序列病毒鉴定工具的比较。在2018年1月1日至2020年7月2日期间,将NCBI参考序列数据库(RefSeq)中保存的完整细菌/古细菌/噬菌体基因组和质粒随机片段化,以均匀分布生成Contigs。然后分别在真阳性(噬菌体基因组片段)和阴性(细菌/古细菌染色体和质粒片段)数据集上运行每个工具。对于分数/概率阈值或类别可以手动调整的工具,根据最优F1分数选择值/类别

假阳性病毒预测的分类在不同的工具之间以及在染色体和质粒数据集之间有所不同。近一半(44.7%)的DeepVirFinder的染色体假阳性命中(fph)属于梭状芽孢杆菌,尽管该类仅占数据集的1.27%(额外文件6)。Seeker在染色体数据集上的假阳性预测偏向于杆菌(30.0%的fph对染色体数据集的13.8%)和放线菌(29.2%的fph对染色体数据集的16.1%)。对于除Kraken2外的所有其他工具,假阳性的分类概况大致与数据集中的总体分布相匹配。Seeker在质粒数据集上的fph主要是属于Halobacteria类的片段(75.6%的fph与3.63%的质粒数据集)(附加文件6)。VirFinder也偏向于Halobacteria质粒片段,占其fph的20.0%。除了Kraken2和ViralVerify(病毒)外,Mollicutes类的质粒片段在每个工具的假阳性命中中都被过多代表,它们分别有0和2个假阳性命中。这种过度代表性在PPR-meta中尤为明显(10.5%的FPHs);充满活力(5.66%的fph);VirSorter所有类别(其fph的11.4%);VirSorter 1、2、4和5类(14.2%的fph);和VirSorter2(占其fph的8.2%),而Mollicutes仅占总质粒片段的0.36%。这些百分比必须在每个工具的整体性能的背景下考虑,因为工具之间的fph总数差异很大-在染色体数据集中从2到4331,在质粒数据集中从0到313。

随机打乱序列的基准测试工具

作为进一步的阴性对照,阳性RefSeq基准序列在核苷酸水平上随机洗牌,以产生不应被任何工具识别为病毒的序列。在测试的工具中,四个识别零洗牌的contigs (metaPhinder, viralVerify, VirSorter和VirSorter2), Kraken2分类了三个contigs, DeepVirFinder和VirFinder分别检测到742和1070,其余工具识别超过2500个洗牌的contigs为病毒,包括PPR-meta,它错误地将99.2%(6608/6664)的洗牌contigs分类为病毒(表3)。

表3工具在随机洗牌人工噬菌体上的表现

在之前的基准测试中生成的人工RefSeq噬菌体片段被随机洗牌,同时使用来自HMMER3套件的esl-shuffle保留二核苷酸分布。然后使用噬菌体检测工具对洗牌后的片段进行检测,并记录任何阳性命中。黑体字的工具表示使用机器/深度学习模型的方法。

使用模拟社区霰弹枪的基准测试工具tagenomes

接下来,我们试图在真实的社区霰弹枪宏基因组配置上比较这些工具。因此,我们获得了一个由Kleiner等人创建的不均匀模拟群落的测序数据,该群落包含32个来自生命树的物种,包括5个噬菌体,细胞丰度范围很大(0.25-21.25%;集成的噬菌体被Phigaro和PhageBoost删除,以防止它们被标记为病毒。在此过程中,没有属于这五个噬菌体的contigs被移除。

这使我们能够评估我们的工具在实际数据上的性能,同时保留对地面真相(样本组成)的了解,并确定每种工具对低丰度物种的检测极限。除了viralVerify和VirSorter之外,RefSeq基准测试中发现的优化参数被用于每个工具。这些工具的分类阈值会极大地改变已识别的病毒contigs的特征,并且阈值之间的F1分数非常接近,因此在该基准测试中进一步分析了参数。总的来说,这些工具在这个数据集上的F1分数比RefSeq人工配置低得多,与RefSeq基准相比,F1分数平均下降了40.6%(图3)。Kraken2的F1分数为0.86,比排名第二的DeepVirFinder高0.3。这是由于它的精度高达0.99,同时识别76%的噬菌体组合。DeepVirFinder的召回率很高(0.80),但与RefSeq基准不同的是,它的准确率较低,为0.42。其他几个工具也有类似的结果,PPR-meta、metaPhinder和VirFinder都分别达到了0.91、0.98和0.83的高精度,但精度分数相对较低(分别为0.24、0.28和0.35)。VirSorter(类别1、2、4和5)和Seeker的性能与上述工具相当,它们的F1分数都是0.44。在这个基准测试中,Seeker和Kraken2是唯一获得与RefSeq相似的F1分数的工具。在RefSeq基准测试中表现最好的VirSorter2,其F1分数(0.36)低于其前身,主要是因为它的假阳性命中次数较多。在之前的基准测试中也表现良好的VIBRANT,同样具有较差的精度和中等召回率,导致F1得分为0.32。viralVerify (NBC数据库,仅包含病毒)具有较低的准确率(0.24)和召回率(0.22),导致该测试中F1得分最低(0.22)。K-mer工具的平均F1得分高于参考相似性/基于基因的工具,但由于本研究中工具数量的样本量较小,这种差异没有统计学意义(p = 0.094,单尾Welch t检验)。

图3
figure 3

不同模拟群落样本病毒鉴定工具的比较。模拟社区读取从先前的研究[56]中检索,并与metaSPAdes组合。在使用基于先前基准测试的最佳阈值运行每个识别工具之前,使用Phigaro和PhageBoost检测和删除噬菌体(viralVerify和VirSorter除外)。F1分数,精度和召回指标显示为单独的面板。每个样本被绘制为每个工具的单个点,用箱线图表示所有三个样本的四分位数范围、极值和平均值

在组装体中发现的四种DNA噬菌体中(由于噬菌体F2是一种ssRNA病毒,因此没有对其进行测序),只有metaPhinder能够检测到M13,而F0是最丰富的噬菌体,在所有样品中通过所有基准工具检测到。metaPhinder在两个样本中鉴定出了属于所有四种噬菌体物种的contigs,在另一个样本中鉴定出了三个。PPR-meta、VIBRANT、VirSorter和VirSorter2能够在所有三个样本中识别出属于三个物种的contigs。viralVerify和VirFinder能够在三个样本中的两个样本中识别出这三种噬菌体,但遗漏了属于噬菌体ES18的片段。DeepVirFinder和Kraken2在三个样本中的一个样本中分类了属于三个噬菌体物种的病毒contigs,并在其他样本中检测到两个物种。Seeker只能识别出属于最丰富的噬菌体F0的contigs。F1评分与检出噬菌体菌种数无相关性(Rs =−0.371,p = 0.29)。然而,在鉴定出更多病毒来源contigs(真阳性+假阳性)的工具与鉴定出的噬菌体菌株数量之间观察到正相关,但不具有统计学意义(Rs = 0.604, p = 0.06)。

工具预测对多样性度量估计的影响

为了测试这些工具对多样性估计的影响,从[57]中检索了四个模拟的模拟社区元病毒组,平均包含719个病毒基因组。将Reads映射到每个工具鉴定为病毒的contigs (bbb1kbp)上,然后将这些映射的Reads映射到一组种群contigs上,以估计它们在每个样本中的丰度。作为对照,原始读数也直接映射到种群序列。然后通过映射的组长度和样本库大小对读取计数进行归一化,Roux等人发现这是一种可靠的归一化方法。然后使用归一化的种群计数计算多样性估计指标。与初始种群相比,所有工具返回的每个样本的基因组都更少,尽管工具之间存在显着差异。PPR-meta、metaPhinder和Kraken2检索到的基因组比例最高,分别为86.8%、89.1%和83.7%(图4A)。除了Seeker和viralVerify,其他所有工具都能检索到50%以上的基因组,它们分别只能检索到32.4%和41.3%的群体基因组。从每个工具的计数矩阵中计算出的Shannon’s alpha多样性都在默认总体的10%以内,除了Seeker,其H分数平均低27.0%(图4B)。Simpson α -多样性指数显示出类似的表现,所有工具的多样性得分都在初始种群的1%以内,但DeepVirFinder和Seeker除外,它们分别为1.4%和5.1%(图4C)。PPR-meta是唯一的工具来估计相对较高的α多样性比默认种群。对于beta多样性,除Seeker外,样本内的两两布雷-柯蒂斯差异较小,其相似性分析(ANOSIM)显示与其他工具相比存在显着差异(r = 0.3926, p = 0.0019,在FDR = 0.05时进行多次比较的Benjamini-Hochberg校正)(图4D;附加文件8)。

图4
figure 4

工具预测病毒群多样性指标的估计。为了评估每种工具对种群多样性的影响,下载了Roux等人的4个模拟病毒体组装。然后运行每个程序来确定预测的病毒组合的子集。读取被映射到这些组簇子集,然后映射的读取随后被映射到种群组簇池。所有多样性指标均由R包“素食主义者”计算。每个图中的“默认”表示每个样本的原始组装。从每个工具的读取映射到预测的病毒contig群,观察到许多基因组。B每个工具的病毒组子集估计香农多样性指数的比较。估计是基于读取计数被归一化的组态大小和测序深度的病毒。C每种工具的病毒组亚群的Simpson多样性指数比较。D各病毒鉴定工具预测的病毒亚群Bray-Curtis差异的非度量多维尺度(NMDS)排序图。椭圆表示每个样本簇质心的95%置信区间。样本用相同的符号和椭圆线型表示;工具用颜色表示

运行时和计算每个工具的负载

我们还在高性能集群(16个VCPU, 108 gb RAM)上记录了每个工具在refseq阳性数据集上的运行时间(图5)。Kraken2和VirFinder是最快的工具,在5分钟内完成。DeepVirFinder, metaPhinder, PPR-meta和Seeker在半小时内完成,而viralVerify刚好完成。VirSorter和VirSorter2在该数据集上运行的时间最长(分别为2.9小时和3.8小时)。

图5
figure 5

正RefSeq人工配置集上工具运行时间的比较。在没有GPU加速的16 VCPU、108 gb RAM和Linux高性能集群上记录每个工具在模拟社区样本上的Wall runtime。工具使用16个线程运行,可以将其设置为参数(除了metaPhinder、PPR-meta和Seeker之外的所有工具)。refseq阳性集合包含约5340万个bp


目录


摘要
介绍
结果
讨论
结论
材料与方法
数据和材料的可用性
缩写
参考文献
致谢

作者信息
道德声明

补充信息



#####

讨论

噬菌体是地球上几乎所有生态系统中微生物群落的重要成员,负责控制宿主种群规模,并对群落功能产生更广泛的影响。设计用于从混合群落宏基因组和病毒组样本中恢复病毒序列的工具是研究噬菌体在其更广泛环境中的作用的基础。这一领域的进步已经产生了一套广泛的病毒识别工具,每个工具都声称提高了类似工具的性能。因此,在这些工具中选择最适合数据集的工具并不简单,特别是每个新工具通常只针对两到三个其他现有工具进行基准测试。为此目的开发的大多数工具,特别是近年来发布的工具,都利用机器/深度学习对序列进行分类,而其他工具则依赖于基于序列与数据库中序列的相似性对序列进行分类。随着新发现的病毒基因组被添加到训练数据集和数据库中,这两种方法都有可能随着时间的推移而得到改进。

在这里,我们比较了从三个数据集的宏基因组中鉴定病毒序列的十种方法。我们首先在正数据集和负数据集上对工具进行基准测试,以评估它们在理想配置集(大小在1到15 kbp之间,没有错集)上的性能,并确定近似的最佳阈值。大多数工具在这里表现良好,检测大多数噬菌体序列,同时保持低假阳性。PPR-meta、VIBRANT和VirSorter2都使用不同的机器学习方法,在所有工具中表现最佳。一般来说,k-mer工具优于参考相似性和基于基因的工具。虽然我们确定的最佳阈值可能并不一定适用于所有其他数据集,但我们相信它们可以作为进一步使用这些工具的基础,因为在每种情况下,它们产生的结果都比默认参数好得多。因此,我们鼓励研究人员在其预期数据集的背景下应用这些阈值和参数。噬菌体分类学不影响病毒预测工具,而染色体和质粒分类学在一些工具中显示出偏差。因此,在分析宏基因组时,重要的是要考虑在采样的微生物生态系统中可能存在哪些噬菌体宿主,以及所选择的病毒预测工具的偏差是否会影响下游分析。

在分析机器学习和更传统的工具对随机洗牌噬菌体序列的识别时,可以看到它们之间的鲜明对比。在利用机器/深度学习方法的6个工具中,有5个工具识别出了很大一部分序列是病毒,只有VirSorter2是唯一的例外,这可能是由于它的分类器在一系列序列和基因特征上进行了训练。其他四个工具中有三个返回零误报,而Kraken2只返回三个误报。这突出表明,虽然这些机器/深度学习算法具有检测新型噬菌体的能力,但当暴露于具有与我们的训练集不同特征的新数据时,它们的性能可能是不可预测的。

当在真实的宏基因组数据上测试时,大多数工具的表现都比RefSeq基准测试差得多,除了Kraken2和Seeker。一般来说,与参考相似性/基于基因的工具相比,k-mer工具在RefSeq基准上的F1分数下降幅度较小。这可能是由于从宏基因组reads中组装的噬菌体序列相对较短,提供的遗传背景较少,从而对参考相似性和基于基因的工具的算法产生负面影响。metaPhinder是唯一能够检测每个样品中组装的噬菌体M13的一个contig的工具。不幸的是,这可能是由于metaPhinder在此基准测试中的低精度导致许多误报调用。大多数其他工具都能够识别属于其他三个噬菌体的contigs,除了Seeker只能识别每个样品中最丰富的噬菌体F0。这表明,在分析噬菌体物种可能处于低丰度的宏基因组数据集时,基于k-mer的工具(如Kraken2和DeepVirFinder)是最佳选择,当精确度特别重要时,前者是受欢迎的选择,而后者的使用是适当的,如果发现新的噬菌体是由于其深度学习算法而感兴趣。

我们还测量了任何潜在的偏差和这些工具可能对其预测病毒种群的多样性产生的影响。大多数工具的α -多样性指数在默认种群的10%以内表现良好,但由于Seeker最初预测的病毒种群基因组数量非常少,因此返回的值相当低。PPR-meta等工具预测的alpha多样性高于默认种群。这是由于这些工具在预测中遗漏了一些高丰度的基因组,导致了更均匀的多样性分布。在评估beta多样性时,Seeker是唯一一个产生与其他工具有显著差异的结果的工具,并且没有与其他程序聚类,这也是由于它在该数据集中恢复的基因组比例低。因此,即使只有一半的基因组被恢复,这里检查的工具的β -多样性趋势(除了Seeker)对原始种群是准确的。

运行时间和计算负载也是需要检查的重要因素,因为如果大型样本需要花费数小时或数天来分析,这些因素可能会成为实际限制。大多数工具都相当快,尽管VirSorter及其后继版本VirSorter2需要几个小时才能完成运行。值得注意的是,VIBRANT、VirSorter和VirSorter2注释已识别的病毒基因组并预测以运行时间为代价的前噬菌体,尽管这些对某些应用程序可能有用。Kraken2是迄今为止最快的工具,在我们的数据集上运行不到一分钟。然而,与其他基准测试工具相比,Kraken2需要非常高的内存使用,因此对于计算能力有限的研究人员来说可能不可行。

虽然这些基准全面比较了最先进的工具的性能,但我们的研究有许多局限性。首先,虽然我们使用RefSeq基因组和模拟宏基因组社区对这些工具进行基准测试,但我们没有解决工具识别属于不同噬菌体家族的病毒序列的能力。其次,我们使用默认数据库或每个工具提供的原始训练模型。虽然为每个工具提供相同的数据库或数据集进行训练,可能会对底层算法进行更公平的比较,但这超出了我们的研究范围。我们注意到,大多数常规用户也不太可能在使用这些工具之前对其进行再培训。第三,我们没有评估组合多个工具的性能,当只使用一个工具时,可能会遗漏一些有意义的见解,就像Marquet et al.[58]中作者将多个工具组合到单个工作流中一样。第四,许多工具都有额外的功能,我们没有在这里进行基准测试,但可能会影响研究人员对工具的选择,如噬菌体预测(VIBRANT, VirSorter, VirSorter2);质粒预测;分类鉴定(Kraken2);和功能注释(dynamic)。最后,我们在研究过程中发现的一些最近开发的工具没有包括在我们的基准测试中,因为(1)需要使用自己的web服务器,因此无法扩展(VIROME, VirMiner),(2)缺乏明确的安装/运行说明(ViraMiner),或(3)尝试使用工具时出现错误,我们无法解决(PhaMers, VirNet, VirMine)。

结论

我们对十种目前可用的宏基因组病毒/噬菌体鉴定工具进行了比较分析,为其他研究人员提供了有价值的指标和见解。使用模拟社区和人工数据集,可以计算这些工具的精度、召回率和偏差。通过调整每个工具的病毒识别过滤阈值并比较F1分数,我们能够在每种情况下优化性能。在人工RefSeq配置基准测试中,大多数工具表现良好,其中PPR meta、VIBRANT和VirSorter2的F1得分最高。在模拟不均匀社区数据集中,工具通常表现较差,但Kraken2除外,其表现包括几乎完美的精度和高于平均水平的召回分数。除了Seeker之外,所有工具都能够产生与原始病毒种群具有相似指数的多样性概况,因此适用于噬菌体生态学研究。我们建议,在目前可用的宏基因组噬菌体鉴定工具中,当研究人员试图鉴定先前表征的噬菌体时,应该考虑Kraken2。当需要新的噬菌体检测时,Kraken2应与VirSorter2和DeepVirFinder等工具结合使用。

材料与方法

使用RefSeq数据集进行基准测试

完整的细菌和古细菌染色体和质粒序列,以及自2020年1月1日(含)以来储存在RefSeq[59]中的噬菌体基因组,于2021年8月12日下载,以构建基准集。使用dedupe.sh[60]删除该集合中与2020年前RefSeq序列和DeepVirFinder、PPR meta、Seeker、VIBRANT、VirFinder和VirSorter训练数据集一致性≥95%的序列,以减少任何潜在的过拟合。然后使用reformat.sh(来自BBTools套件)[60]以14.3的倍数随机取样染色体和质粒序列,得到宿主与噬菌体的比例约为19:1。在染色体和质粒序列上依次运行Phigaro (v2.3.0,默认设置)[54]和PhageBoost (v0.1.7,默认设置)[55],去除前噬菌体序列。开放阅读帧(orf)与pVOG数据库的HMM命中率≥30%的宿主序列被作为污染去除。然后使用自定义python脚本(可在https://github.com/sxh1136/Phage_tools上获得)将所有序列均匀地分割为1到15 kbp,以创建人工配置。然后使用默认设置在三组contigs(染色体、质粒和噬菌体)上运行每个病毒预测工具,除了dynamic和VirSorter,其中噬菌体衍生的contigs集使用它们的病毒去污模式额外运行,因为它们在主要由病毒片段组成的数据集中可能提高性能。可以在https://github.com/sxh1136/Phage_tools/blob/master/manuscript_tools_script.md上找到用于运行每个工具的命令及其版本号。RefSeq基准测试数据集可在https://figshare.com/articles/dataset/RefSeq_Datasets_for_benchmarking_-_Ho_et_al_/19739884下载和使用。

对于分数/概率阈值可以手动调整的工具(DeepVirFinder、metaPhinder、PPR meta、Seeker和VirFinder),绘制F1曲线(100个数据点),并通过最大F1分数确定最佳阈值。对于具有噬菌体鉴定分类阈值的viralVerify和VirSorter,绘制了f1得分最高的类别集;对于VirSorter,比较了两个类别集,类别1、2、4和5以及所有类别,因为这些都是常用的。ViralVerify还使用Pfam-A 33.0和提供的病毒/染色体特异性hmm数据库运行,这些在该工具的GitHub使用指南中列出。Kraken2使用预先构建的Kraken2微生物数据库运行,该数据库可在https://lomanlab.github.io/mockcommunity/mc_databases.html上获得。从NCBI RefSeq下载Kraken2鉴定为病毒的两条染色体片段的参考基因组(NC_CP059254.1和NZ_CP07771.1)。片段和参考基因组用pharokka进行注释[61],并用熟料进行全局对齐和可视化[62]。使用R包taxonomizr v 0.10.3检索工具真阳性和假阳性病毒预测的分类[63]。

本数据集上各工具的运行时间记录使用由Cloud Infrastructure for Big Data Microbial Bioinformatics (CLIMB-BIG-DATA)提供的Linux虚拟机,配置CPU: Intel®Xeon®Processor E3-12xx v2 (8vcpu), GPU: Cirrus Logic GD 5446,内存:64gb多位ECC。工具使用16个线程运行,可以将其设置为参数(除了metaPhinder、PPR-meta和Seeker之外的所有工具)。

随机洗牌序列的基准测试

使用HMMER3套件(v3.3.2, -d)中的esl-shuffle,在保留二核苷酸分布的同时,对先前基准中创建的所有人工噬菌体组合进行随机洗牌[64]。然后使用RefSeq基准测试中确定的优化阈值,在随机打乱的序列上运行每个识别工具,并记录假阳性。

以模拟社区为基准tagenomes

从欧洲核苷酸档案(BioProject PRJEB19901)中检索了一个不均匀模拟群落[56]的三个散弹枪元基因组测序重复。这些群落包含5种噬菌体菌株:ES18 (H1)、F0、F2、M13和P22 (HT105)。使用FASTQC (v0.11.8,默认设置)[65]检查数据质量,并使用cutadapt (v2.10, -max-n 0)删除过度表示的序列[66]。然后用metaSPAdes (v3.14.1,默认设置)组装清洗过的对端读取[67],并删除< 1 kbp的contigs。按照RefSeq染色体和质粒数据集,依次运行Phigaro和PhageBoost从contigs中去除噬菌体。然后使用之前确定的最佳参数在三组配置上运行每个工具,但viralVerify和VirSorter除外,其中所有分类阈值都被重新评估。使用metaQUAST (v5.0.2,默认设置)[68]将contigs映射到参考噬菌体基因组并计算覆盖率。

模拟模拟病毒群落的基准测试

四个模拟群落(样本13、2、7和9)包含500到1000个病毒基因组,由Roux等人创建。这些样本属于四个不同的beta多样性组,并且不共享其中最丰富的50种病毒中的任何一种。使用Trimmomatic[69]对1000万对末端reads的每次模拟进行质量控制,并使用Roux等人的metaSPAdes进行组装。然后下载这些配置以进行基准测试。与之前一样,长度< 1 kbp的contigs被移除,然后输入到每个病毒鉴定程序中。然后提取每个工具的阳性病毒配置集,并用BBMap[60]将读取的片段映射到这些片段上,并随机将模糊映射的读取片段分配给配置(ambiguous = random),如Roux et al.[57]所示。然后使用SAMtools (v1.11)提取映射到相同组(options - f0 × 20 × 904)的主映射读取,并将其映射到一个非冗余填充组池。该池是通过将所有四个样本与nucmer (v3.1)[70]聚类而成,在≥80%的长度范围内,ANI(平均核苷酸同一性)≥95%。每个工具的丰度矩阵通过按序列长度和总库大小规范化读取计数来计算,以产生每百万计数(CPM)。然后使用这些丰度矩阵计算使用纯素包(v2.5.7)的香农、辛普森和布雷-柯蒂斯不相似指数[71]。非度量多维尺度(NMDS)和相似度分析(ANOSIM)也与素食计算。用Benjamini-Hochberg方法校正ANOSIM p值[72]。如果可能,种子和排列分别设置为123和9999。所有图均使用ggplot2 (v3.3.2)[73]生成,并使用ggpubr (v0.4.0)中的ggarrange进行排列。[74]。

补充信息

Additio文件1:补充图1。

工具的优化曲线导致得分/概率阈值。

Additio文件2:补充图2。

提供分类阈值的工具的f1得分图。

Additio文件3:补充图3。

人工RefSeq序列与参考基因组的比较。

Additio文件4:补充图4。

各工具预测RefSeq噬菌体人工组合体为病毒的颠覆图。

Additio文件5:补充图5。

各工具对RefSeq噬菌体片段病毒预测的分类分析。

Additio文件6:补充图6。

RefSeq染色体和质粒片段假阳性病毒预测的分类。

Additio文件7:补充表1。

Kleiner等[j].不均匀模拟群落的噬菌体菌株组成。

Additio文件8:补充表2。

非度量多维尺度工具间相似性分析。



Download: