宏病毒组组装软件如何选？
上一期给大家介绍了宏病毒组的一般分析思路，这期一起来考察一下分析流程中重要的“组装”环节常用的众多软件。
通过我们这段时间来的宏病毒专题介绍，想必大家也都发现了，自然界微生物群落中的病毒组分在驱动细菌多样性、促进养分转化和形成群落组成方面发挥着重要作用。尽管它们很重要，但绝大多数病毒序列注释较差，特别是未知组成复杂的环境样本，与参考数据库几乎没有太高的同源性。因此，对病毒宏基因组（virome）的研究在很大程度上依赖于短测序reads的从头组装来恢复组成和功能信息。宏基因组组装对于病毒组数据尤其具有挑战性，通常导致组装碎片化和病毒群落成员恢复不完全。今天我们就通过2019年发表在Microbiome上的迄今为止对宏病毒组装相关软件测评最全面的一项研究[1]，来一起了解宏病毒组数据组装过程中的挑战与策略。
老规矩，软件评测我们先看测试使用到的计算平台。
除了Geneious和CLC，其他组装软件都是使用联想x3650 M5服务器，处理器为intel Xeon E5-2690 v3，内存512Gb，操作系统为ubuntu 14.04.5。CLC和Geneious运行在64位Win10电脑上，CPU为i5-4690, 内存为16Gb。
数据来源
mock群落A、B的测序reads来自于[2]，simulated病毒组数据集来自于[3]，用于比较测序深度对时间和RAM使用影响的reads来自于[4], 人类病毒组（包含107 PFU的乳球菌噬菌体Q33标准品）数据来源于[5]。
下面我们来看看数据的基本处理步骤
对上述数据集，用FastQC v0.11.5评估Raw read 质量，用cutadapt v1.9.1去除测序接头。使用Trimmomatic v0.36在特定的参数下对数据集进行修剪和过滤。两个mock群落的reads的滑动窗口大小为4，最小Phred得分为30，最小长度为60 bp。从健康人肠道噬菌体reads中去除前端15bp和后端60bp，滑动窗口为4bp，最小Phred评分为20。Q33 -spiked的病毒组reads中删除了前端的10 bp和后端的100 bp，滑动窗口大小为4 bp，最小Phred得分为30。所有数据除去短于60bp的reads。
对来自Q33-spiked数据集的高质量reads进行组装。利用Blastn将Contigs与已公开的Q33进行比对(e值阈值为1e-20)。Q33基因组的比对结果（要求比对长度大于800个碱基，并且至少具有95%的相似度）使用QUAST（v4.4， --unique mapping flag）进行进一步分析。使用Mauve (v. 20150226，build 10)进行Q33组装结果的进一步比较和可视化。使用MetaQUAST (v. 4.4)进行mock和simulated数据集与参考基因组集合的比对（--unique mapping，最小contig长度为500 bp，最小比对长度为65 bp，最小相似度阈值为95%）。使用Spearman方法进行相关性分析，并使用R（v.3.4.3）中的ggplot2（v 3.0.0）软件包进行结果的可视化。利用R中的线性模型验证了上述相关性结果。对于非正态分布的数据，进行Log转换。
表1本研究中使用的组装软件列表
研究结果
Simulated病毒组数据集
对选取的Simulated病毒组数据集（共涉及572个病毒成员）使用多种组装软件进行组装，从组装结果中统计每个病毒成员的基因组上contigs覆盖情况，评估回比的基因组覆盖状态和碎片化程度（图1）。其中MetaVelvet由于7天时间尚未分析完毕，因此未包含在此部分分析结果中。群落中大约一半的基因组的平均基因组覆盖度小于75%，并且在所有组装体中表现出高度片段化（平均每个基因组>10 contigs）。在572个基因组中，有87个基因组在所有组装软件中平均基因组覆盖度不到20%，其中有84个以低丰度存在（Simulated病毒组数据集中的基因组丰度使用基因组长度做过标准化处理）。
图1 每个基因组组装后产生的contigs数目以及在基因组上的覆盖百分比(基因组比例)
大部分组装软件显示群落内的标准化基因组丰度与所有组装软件中基因组覆盖百分比呈显著正相关，并通过线性模型的验证，单SOAPdenovo2的组装结果呈现负相关。除了Velvet的标准化丰度与碎片化程度（contigs 数量）呈正相关，Geneious相关性不显著外，其他组装软件的标准化丰度与碎片化程度（contigs 数量）都呈负相关关系。在正常丰度低于30%的基因组中，没有一个基因组的平均恢复率超过75%，进一步说明了低丰度测序对覆盖率的影响。然而，高丰度也并不能始终如一地提高基因组覆盖率，在标准化丰度top30的172个基因组中，20个基因组的平均基因组覆盖低于50%。对数转换（由于极端值）归一化丰度与平均丰度的距离与所有组装软件中基因组覆盖度呈负相关（相关系数-0.42，p值<2.2e -16）。MIRA和Geneious用比其他组装软件用更少的contigs覆盖了更大比例的低丰度基因组。MIRA组装的群落中13个高丰度的基因组(最高10%)显示出最高程度的碎片化，每个基因组产生401至2983个contigs。
计算每个基因组的反向重复序列，回文重复序列，串联重复序列和基因组总重复序列的比例。除了Ray Meta之外，其他组装软件中观察到的碎片化程度与每个基因组中预测的总重复区域的百分比正相关，并且除了ABySS（k-mer 63/127），Geneious和SOAPdenovo2外，与所有组装软件中基因组覆盖度负相关。
VICUNA, Ray Meta, SOAPdenovo2, Geneious, ABySS (both k-mer sizes) 和Velvet的总基因组覆盖率低于50%（群落中的所有基因组）。VICUNA一共产生了4个高水平错配的contigs(平均每100 kb有174个错配)，这可能与人工reads（artificial reads）的数据特征有关，因为在真实的测序数据中没有观察到这种情况。总体上基因组覆盖率最高的五个组装软件是SPAdes (default)，MEGAHIT，SPAdes (single cell)，SPAdes (single cell + careful)和CLC。单contig覆盖基因组90%以上的基因组总数如图2显示，排名靠前的包括SPAdes、CLC以及MEGAHIT。
图2 不同组装软件单个contig覆盖基因组（覆盖率>=90%）数目
Mock群落数据集
两个mock病毒群落被用来研究高丰度和低丰度ssDNA病毒对组装性能的影响。mock群落A(表2a)包含12个病毒基因组，其中10个具有相同丰度(占总菌群的9.82%)，2个ssDNA基因组(NC_001330和NC_001422)具有低丰度(0.92%)。该部分分析显示，尽管一些组装软件，如CLC、Geneious、SPAdes (single cell)和VICUNA检测到了全部12个基因组，但存在大量的假阳性(分别为1143、53、1513和4969)。Velvet和MetaVelvet没有产生假阳性，但未能组装三个基因组，ABySS (k-mers 63和k-mers 127) 产生了大量的假阳性结果，并分别未能组装4个和6个基因组。IDBA UD和Ray Meta的表现优于其他的组合，其基因组的contigs数量相等（12），其次是MEGAHIT，SPAdes(default) 和SPAdes (meta)，分别为13，14和14。mock群落B(表2b)也含有12个基因组，但ssDNA基因组NC_001330和NC_001422丰度较高(32.47%)。VICUNA对mock群落B的组装效果优于mock群落A，因为没有产生假阳性，而MIRA在mock群落A中的组装假阳性率在从群落B中的0增加到94。基于灵敏度和contigs数量，IDBA UD表现最佳，其次是SPAdes（默认），Ray Meta，MEGAHIT和SPAdes（Meta），而ABySS（k-mer大小）和SOAPdenovo2的灵敏度最低。尽管是一个仅仅由12名成员组成的相对简单的群落，但并非所有组装软件都能够组装覆盖所有病毒成员。与mock群落A中四个软件未能完整组装出实际成员相比，mock群落B中有六款软件未能组装出所有成员。ABySS(k-mer 63)、ABySS(k-mer 127)、Velvet和MetaVelvet在mock A中分别有6、4、3和3个基因组组装失败，在mock B中分别有6、4、1和1个基因组组装失败。此外，MIRA和SOAPdenovo2在mock群落B中分别有1个和2个基因组组装失败。
表2 每个组装软件对A mock群落a和B mock群落b产生的假阳性和假阴性contigs的数量以及它们各自的灵敏度
MEGAHIT除了ssDNA NC_001422的基因组在两个contigs中只恢复了56.5%之外，在mock群落A中至少组装了98%的基因组，在mock群落B中也组装了超过98%的基因组。
Q33
有5个组装软件：ABySS (k-mer 63)，ABySS(k-mer 127)， SOAPdenovo2， Velvet 和 MetaVelvet未能生成符合比对阈值的contigs，随后被排除在进一步的分析之外。除MIRA(8.5%)外，所有剩余的组装软件都恢复了超过90%的Q33基因组。一共有六个组装软件在单contig中恢复了超过90%的Q33基因组，分别是SPAdes (meta) 99.74%，MEGAHIT (99.6%)，VICUNA (99.6%)， Ray Meta (99.6%)， CLC (99.5%)和Geneious (99.1%) (图3)。但是，只有MEGAHIT组装了Q33基因组，其contigs的长度与基因组本身相等。SPAdes（meta）和CLC产生的组装体比参考基因组短86和141个碱基。其他四款软件都产生了比参考基因组长的组装体：VICUNA (723)， Geneious(1765)和Ray Meta(9884)。此外，SPAdes (default)、SPAdes (single cell)、IDBA UD 和 SPAdes (single cell +careful)分别组装出2、3、4、5和5个contigs。而Ray-Meta和VICUNA每100 kb具有最低数量的错配和插入缺失；但是Ray-Meta表现出最高的错配率。除了SPAdes (meta)之外，所有的组装软件都有一个最小的局部组装错误。
图3 六个组装软件组装出的Q33参考基因组(单个contig达到>99%基因组覆盖)。
reads深度分析(所用时间和内存)
通过测试四个已发表的不同的测序深度的健康人肠道病毒的组装完成时间和内存(RAM)的最大消耗，比较组装软件的实用性。此处需要提及的是所有组装任务都分配了五个线程。
但是，指定线程数量并不会改变某些软件使用的线程数量。MetaVelvet运行七天后还未完成，因此不包括在这部分分析中。CLC和Geneious是在台式WIN系统计算机上进行的，因此也不包括在这部分分析中。组装运行时，时间的消耗很大程度上受限于reads的数量，并且很大程度上是线性的影响，reads数越多，组装时间越长(图4a)。MIRA和VICUNA（图4a）是最慢的，MIRA比其他软件花费超过15倍的时间来组装350万条序列。SOAPdenovo2的完成时间最短，其次是IDBA UD和Velvet。除了MIRA和Ray Meta之外，其他大多数组装软件组装速度基本一致。MIRA，VICUNA和Velvet（图4b）具有最高的最大RAM使用量，而最低的是Ray Meta，IDBA UD和SPAdes（Meta）（图4b）。
reads深度分析N50和最长contig长度
对于N50(图4c)和最大contig长度(图4d)，组装软件之间大都存在较大的差异。随着序列深度的增加，若干组装软件（SPAdes (default)，SPAdes (meta)， MEGAHIT 和 ABySS (k-mer 127）)能产生更长的contigs。
图4 不同软件在组装不同深度的测序数据中消耗的时间、RAM以及contigs序列集的N50指数和最大contig长度的比较。
探讨
在本测试研究中使用的所有组装软件中，SPAdes（meta）与MEGAHIT显示出来不错的性能。然而，所有的组装软件都有局限性，并且受到病毒组数据各方面数据特征的影响。比如，低reads覆盖率和高基因组的重复区域导致组装结果全基因组的低覆盖及过高程度的片段化，组装本身就会变得不准确。因此未来病毒组研究的设计应仔细考虑测序深度的影响，因为read覆盖范围的极端情况将阻碍高丰度和低丰度病毒基因组的组装和检测。最后，作者也希望各位生信大佬们尽快开发出专门针对宏病毒组的专业组装软件，来提高组装效率和质量。
Reference
[1] Sutton TDS, Clooney AG, Ryan FJ, Ross RP, Hill C: Choice of assembly software has a critical impact on virome characterisation. MICROBIOME. 2019;7(1).
[2] Roux S, Solonenko NE, Dang VT, Poulos BT, Schwenck SM, Goldsmith DB, Coleman ML, Breitbart M, Sullivan MB. Towards quantitative viromics for both double-stranded and single-stranded DNA viruses. PeerJ. 2016;4:e2777.
[3] Hesse U, van Heusden P, Kirby BM, Olonade I, van Zyl LJ, Trindade M. Virome assembly and annotation: a surprise in the Namib Desert. Front Microbiol. 2017;8:13.
[4] Manrique P, Bolduc B, Walk ST, van der Oost J, de Vos WM, Young MJ. Healthy human gut phageome. Proc Natl Acad Sci. 2016;113(37):10400–5.
[5] Shkoporov AN, Ryan FJ, Draper LA, Forde A, Stockdale SR, Daly KM, McDonnell SA, Nolan JA, Sutton TD, Dalmasso M. Reproducible protocols for metagenomic analysis of human faecal phageomes. Microbiome. 2018; 6(1):68.