四種de novo組裝軟件對柞蠶微孢子蟲全基因組組裝結果的比較

劉宗林 許金山＊ 周澤揚 ，3

（1.重慶師范大學 重慶市動物生物學重點實驗室，重慶 400047；2.重慶師范大學 活性物質(zhì)生物技術教育部工程研究中心，重慶 400047；3.西南大學 家蠶基因組生物學國家重點實驗室，重慶 400716）

微孢子蟲（Microsporidia）是一類專性寄生于細胞內(nèi)的單細胞真核生物，可以感染包括昆蟲、魚類和人類在內(nèi)的幾乎所有的動物［1］。感染柞蠶的柞蠶微孢子蟲（Nosemaantheraeae），是柞蠶微粒子病的病原，其可經(jīng)母蛾傳染子代，給柞蠶養(yǎng)殖業(yè)生產(chǎn)造成嚴重損失［2］。目前對柞蠶微粒子病的防治主要是通過胚種檢疫杜絕母體傳染和卵面與養(yǎng)蠶環(huán)境消毒。從分子水平探究柞蠶微孢子蟲對柞蠶的侵染機制，將有助于建立高效、準確的病害檢疫和防治方法。ABYSS、Velvet，SOPAdenovo和Ray是目前常用的四種de novo組裝拼接軟件。Abyss［3］是基于Bruijn算法的拼接軟件，Abyss拼接的優(yōu)點是可以進行平行運算，并且可以多線程以及多個庫同時運行，它往往適用于大基因組的短paired－end reads。有研究者曾利用Abyss對非洲男性的基因組測序得到的35億對reads進行組裝，獲得了覆蓋度達到了68%的全基因組框架圖譜［4］。SOAPdenovo［5］是華大基因開發(fā)的高通量從頭測序軟件，能構建人類基因組大小的從頭拼接草圖。SOAPdenovo使用的Bruijn算法，是專門為處理Illumina GA產(chǎn)生的短reads而開發(fā)的，其為構建參考基因組提供了新的途徑［6］。Velvet［7］是一款常用的基因組拼接軟件，采用的也是Bruijn算法，它能同時支持fasta、fastq格式的數(shù)據(jù)，同時支持多個文庫的數(shù)據(jù)同時使用。Velvet僅僅利用短reads和paired－ends信息，就可以產(chǎn)生可觀長度的重疊群，Vevlet是目前短序列組裝的應用較多的軟件，對細菌基因組十分適合。其工作的一般過程簡化為：輸入短read序列，排除錯誤，產(chǎn)生高質(zhì)量的contigs。然后用成對reads信息，檢索contigs之間的重復區(qū)域［8－9］。Ray［10］是需要使用MPI2.2來進行de novo組裝的軟件。Ray是構建在RayPlatform基礎上的組裝軟件，可以進行平行運算以及多個庫同時運行。它可以同時拼接不同測序軟件的reads，在對應的軟件下面有重要的參數(shù)說明。Ray可以同時拼接Roche/454和Illumina這三種測序工具測序得到的三種不同的基因組序列，結果表明這種混合拼接技術可以大大的減少拼接結果的錯誤［11］。目前本研究團隊完成了柞蠶微孢子蟲的染色體核型分析［12］，同時也初步完成了柞蠶微孢子蟲全基因組的測序，但對數(shù)據(jù)有效全基因組組裝，還處于空白。昆蟲微孢子蟲基因組一般具有大量的重復序列［13－14］，加大了組裝的難度，所以組裝軟件和組裝參數(shù)的選擇顯得十分重要。本文將通過比較不同的組裝軟件和組裝參數(shù)對柞蠶微孢子蟲的從頭組裝結果，通過比較N50的長度，contig的條數(shù)，contigs的總長度等來評價組裝結果，從中選擇最優(yōu)的組裝軟件和組裝參數(shù)。

1 材料和方法

1.1 基因組數(shù)據(jù)來源

實驗材料來源于重慶師范大學資源昆蟲及其病原微生物學研究室測序完成柞蠶微孢子蟲基因組數(shù)據(jù)。

1.2 組裝關鍵參數(shù)選擇

四種組裝軟件 Velvet（1.2.07）、Abyss（1.3.4）、SOAPdenovo（v1.04linux 32）Ray（v2.0.0）組裝時涉及到的參數(shù)很多，但是k－mer是其中最為關鍵的參數(shù)之一，因此，本實驗設計了不同的k－mer進行組裝結果的比較，而其他參數(shù)設為默認。通過perl程序來測試不同的k－mer值.對于Velvet，Abyss，SOAPdenovo，Ray四種軟件我們也都通過編寫perl腳本程序，完成其他統(tǒng)計分析。

2 結果分析

2.1 不同的k－mer下基因組組裝質(zhì)量比較

采用四種軟件組裝后，結果顯示隨著k－mer的增加，組裝所得到的基因組全長也是逐漸增加的，最大基因組總長度在5～6M之間（如圖1a所示）。我們進一步考察了四種軟件在各自不同k－mer值下最大的contig長度，如圖1b所示，對于Ray和SOAPdenovo兩種軟件而言，當k值設定在21～31的范圍內(nèi)，k－mer對最長的contig大小影響不大。而對于Abyss來說，當k＝23時最長的contig大小明顯變大，而k＞23時這種變化趨穩(wěn)。對于Velvet來說，不同的k值對最長contig大小影響最為顯著，具體來說，當k＝25時最長contig大小明顯高于k－mer＝21和k－mer＝23，而在k＝31時最為顯著。而比較四種軟件在k－mer數(shù)值最優(yōu)化后組裝的最長contig長度，Velvet所獲得的結果明顯優(yōu)于另外三種軟件。

N50是反映拼接效果最重要的參數(shù)，N50越大表明組裝結果也越好。如圖1c所示，在不同的k－mer下四種軟件得到的N50差別很大，N50Velvet＞N50Abyss＞N50Ray＞N50SOAPdenovo。從圖中我們可以看出，在設定的k－mer范圍內(nèi)N50Velvet總體上時隨著k－mer的增加而增加的，而N50Abyss是先增加后減小，在k－mer＝25時最大，而N50Ray、N50SOAPdenovo卻沒有明顯變化。Velvet的N50最大達到了47295bp，明顯高于其他三種軟件，而N50SOAPdenovo最小，小于1000bp。由此可見不同的k－mer值對不同軟件的拼接結果有顯著的影響，綜合以上分析，我們列出了不同軟件在最優(yōu)參數(shù)選擇下的組裝質(zhì)量結果，如表1所示。

圖1 不同組裝軟件下基因組組裝質(zhì)量

表1 四種軟件最優(yōu)k－mer化的拼接結果比較

2.2 四種軟件的硬件資源消耗

我們比較了四種軟件在最優(yōu)組裝條件下的系統(tǒng)資源消耗情況，如圖2a、圖2b所示，Velvet對CPU的消耗是相對小的，而內(nèi)存的消耗相對要高一點。而Abyss對硬件的消耗相對是最小的。拼接速度上來講，Abyss、Velvet速率相對較快，而Ray拼接最慢。

通過上面的分析可知，對于柞蠶微孢子蟲的基因組，通過比較N50的長度，最大的contigs大小，contigs的總長度等參數(shù)得知，Velvet的拼接結果相對最好，最適合柞蠶微孢子蟲的基因組組裝。在硬件消耗方面Velvet是最高的，但是拼接速度是最快的。

圖2 四種軟件的硬件資源消耗

3 討 論

本研究通過比較分析可知，不同的拼接軟件對柞蠶微孢子蟲的組裝結果差別非常大。其中Velvet軟件更適合進行柞蠶微孢子蟲的基因組組裝。究其原因分析可能是由于柞蠶微孢子蟲基因組自身的結構特殊性，即包含大量的重復序列，導致不同拼接的結果差異比較大，所以合理的選擇拼接軟件和參數(shù)對于拼接昆蟲微孢子蟲基因組十分重要的。通過本次研究，對今后柞蠶微孢子蟲基因組的分析具有十分重要的理論和實踐意義。

［1］ Keeling P.Five questions about Microsporidia［J］.PLoS Pathog，2009，5（9）：e1000489.

［2］ 中國農(nóng)業(yè)科學院蠶業(yè)研究所.中國養(yǎng)蠶學［M］.上海：上海科學技術出版社，1990：847－851.

［3］ http：//www.bcgsc.ca/platform/bioinfo/software/abyss.

［4］ Simpson J T ，Wong K，Jackman S D.Abyss：aparallel assembler for short read sequence data［J］.Genome Res 2009，19：1117－1123.

［5］ http：//soap.genomics.org.cn/soapdenovo.

［6］ Li R Q，Zhu H M，Ruan J，Qian W B，et al.Denovo assembly of human genomes with massively parallel short read sequencing［J］.Genome Res 2009，20：265－272.

［7］ http：//www.ebi.ac.uk/～zerbino/velvet/.

［8］ Zerbino D R，Birney E.Velvet：algorithms for de novo short read assembly using deBruijn graphs［J］.Genome Res 2008，18：821－829.

［9］ Zerbino D R.Using the Velvet de novo assembler for short－read sequencing technologies［M］.Curr Protoc Bioinformatics 2010.

［10］https：//github.com/sebhtml/ray.bastien.

［11］Se′bastien boisvert F J，L.Aviolettel，Corbeil J.Ray：Simultaneous Assembly of Reads from a Mix of High－Throughput Sequencing Technologies［J］，Computational Biology，2010，（17）：1519－1533.

［12］Xu J S，Wang L J，Zhou Z Y.The Nuclear Apparatus and Chromosomal DNA of the Microsporidian Nosema antheraeae［J］.Eukaryotic Microbiology 2011，58（2）：178－180.

［13］Xu J S，Wang M，Zhang X Y，et al.Identification of NbME MITE families：Potential molecular markers in the microsporidia Nosema bombycis［J］，Invertebrate Pathology，2010，（103）：48－52.

［14］Cornman R S，Chen Y P，C.Schatz M C，Street C，Zhao Y，et al.Genomic Analyses of the Microsporidian Nosema ceranae，an Emergent Pathogen of Honey Bees［J］，PLoS Pathogens，2009，5（6）：e1000466.