宏病毒组组装软件如何选？

上一期小编给大家介绍了宏病毒组的一般分析思路，这期我们就一起来考察一下分析流程中重要的“组装”环节常用的众多软件。

通过我们这段时间来的宏病毒专题介绍，想必大家也都发现了，自然界微生物群落中的病毒组分在驱动细菌多样性、促进养分转化和形成群落组成方面发挥着重要作用。尽管它们很重要，但绝大多数病毒序列注释较差，特别是未知组成复杂的环境样本，与参考数据库几乎没有太高的同源性。因此，对病毒宏基因组（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。

下面我们来看看数据的基本处理步骤
对上述数据集，用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本研究中使用的组装软件列表

大部分组装软件显示群落内的标准化基因组丰度与所有组装软件中基因组覆盖百分比呈显著正相关，并通过线性模型的验证，单SOAPdenovo2的组装结果呈现负相关。除了Velvet的标准化丰度与碎片化程度（contigs 数量）呈正相关，Geneious相关性不显著外，其他组装软件的标准化丰度与碎片化程度（contigs 数量）都呈负相关关系。在正常丰度低于30%的基因组中，没有一个基因组的平均恢复率超过75%，进一步说明了低丰度测序对覆盖率的影响。然而，高丰度也并不能始终如一地提高基因组覆盖率，在标准化丰度top30的172个基因组中，20个基因组的平均基因组覆盖低于50%。对数转换（由于极端值）归一化丰度与平均丰度的距离与所有组装软件中基因组覆盖度呈负相关（相关系数-0.42，p值

计算每个基因组的反向重复序列，回文重复序列，串联重复序列和基因组总重复序列的比例。除了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。

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)之外，所有的组装软件都有一个最小的局部组装错误。

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.

了解详细：原文解读