基于reads引導(dǎo)的基因組序列拼接算法

曾培龍

摘 要：針對(duì)新一代測(cè)序技術(shù)數(shù)據(jù)讀取片段reads長(zhǎng)度短、準(zhǔn)確度低、數(shù)據(jù)海量等特點(diǎn)，本文提出了基于reads引導(dǎo)的基因序列拼接算法（SRGA），以整條reads為拼接單位，并首次提出了基于數(shù)據(jù)特征和拼接信息累計(jì)的評(píng)分機(jī)制。選取常用測(cè)試集，將本文中的算法與序列拼接領(lǐng)域中的經(jīng)典算法進(jìn)行對(duì)比和分析，取得了較好的效果。

關(guān)鍵詞：生物信息學(xué)；新一代測(cè)序技術(shù)；基因組序列拼接

中圖分類號(hào)：TP391 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：2095-2163（2015）03-

GENOME ASSEMBLY GUIDED BY READS

ZENG Peilong

（China Ship Development and Design Center， WuHan 430064，China）

Abstract：Due to next generation sequencing data of mass， short length and relatively low precision， this paper proposes a new genome assembly guided by reads， regarding one entire reads sequences as assembly unit. This algorithm firstly invents a scoring mechanism based on accumulated assembly information and data charactistics. Then the paper gives the metrics results of several algorithms on the test set， the proposed （SRGA） and several classical algorithm of genome assembly. Experimental results show SRGA can obtain satisfactory stereo matching results.

Key words：Bioinformatics， Next-generation Sequencing， Genome Assembly

0 引 言

新一代測(cè)序技術(shù)促進(jìn)了生命科學(xué)的快速發(fā)展，但其產(chǎn)生的基因讀取片段reads具有長(zhǎng)度短、準(zhǔn)確度低、數(shù)據(jù)海量等特點(diǎn)[1-2]，這就對(duì)序列拼接算法提出了相當(dāng)嚴(yán)峻的挑戰(zhàn)，傳統(tǒng)的序列拼接軟件已不再適用[3]。為此，即需針對(duì)新一代測(cè)序的數(shù)據(jù)特點(diǎn)，從實(shí)際應(yīng)用需求出發(fā)，研發(fā)新的優(yōu)質(zhì)高效的序列拼接軟件。

本文針對(duì)新一代測(cè)序數(shù)據(jù)的數(shù)據(jù)特點(diǎn)，提出了基于reads引導(dǎo)的基因組序列拼接算法（SRGA），并以整條reads為拼接單位，首次提出了基于數(shù)據(jù)特征和拼接信息累計(jì)的評(píng)分機(jī)制，從而減少不必要的重復(fù)計(jì)算，同時(shí)也提高了基因組序列拼接的質(zhì)量和速度。

1 reads數(shù)據(jù)預(yù)處理

在序列拼接過程中，reads數(shù)據(jù)的預(yù)處理具有重要意義。由于新一代測(cè)序數(shù)據(jù)精確度較低，以及數(shù)據(jù)海量，造成reads中含有大量的堿基錯(cuò)誤。理論上講，De Bruijn圖構(gòu)建時(shí)，圖的大小規(guī)模只與基因組的大小相關(guān)，與reads數(shù)據(jù)量無(wú)關(guān)。但在堿基錯(cuò)誤的影響下，De Bruijn圖[4]實(shí)際大小會(huì)隨著reads數(shù)據(jù)的增加呈現(xiàn)幾何型增長(zhǎng)。拼接時(shí)，reads中的堿基錯(cuò)誤也極易導(dǎo)致錯(cuò)誤拼接。因此，拼接前需要對(duì)數(shù)據(jù)進(jìn)行預(yù)處理，去除初始數(shù)據(jù)中的錯(cuò)誤堿基。

（1）新一代測(cè)序數(shù)據(jù)的準(zhǔn)確率較低，錯(cuò)誤堿基主要分布在reads 3端，并且越靠近3端錯(cuò)誤率越高，reads 3端錯(cuò)誤率更高，接近20%，而在5端則非常準(zhǔn)確[5]，如圖1所示。為降低錯(cuò)誤堿基對(duì)拼接的影響，拼接前需過濾掉出錯(cuò)率較高的堿基數(shù)據(jù)。處理方法為：以靠近3端二分之一reads長(zhǎng)度的堿基序列為基準(zhǔn)，計(jì)算該區(qū)域堿基序列的質(zhì)量平均值，若該值小于15，則過濾掉該條reads。該平均值對(duì)應(yīng)堿基的錯(cuò)誤率，計(jì)算公式為：

（1）

其中，Q為堿基質(zhì)量值， 為堿基出錯(cuò)率。

圖1 Solexa數(shù)據(jù)錯(cuò)誤分布

Fig.1 Error display of Solexa data

（2）測(cè)序過程中往往會(huì)產(chǎn)生許多人工數(shù)據(jù)[6]，這些reads數(shù)據(jù)會(huì)有許多標(biāo)識(shí)為A的堿基序列，需要去除。處理方法為：若某條reads中A含量>=0.9，則該reads被過濾掉。

（3）測(cè)序過程中有時(shí)會(huì)產(chǎn)生一些沒有被測(cè)出來的通常表示為“N”或“.”未知堿基[7]，需要去除。處理方法為：如果某條reads中含有未知堿基，則該條reads被過濾掉。

2 De Bruijn圖構(gòu)建

SRGA以reads為輸入，經(jīng)過數(shù)據(jù)預(yù)處理、De Bruijn圖構(gòu)建、長(zhǎng)序列片段contigs拼接等，生成基因組序列。以reads為拼接單位，考慮所有正在參與拼接的reads信息，采用基于數(shù)據(jù)特征和拼接累計(jì)信息的評(píng)分機(jī)制，選擇最佳的路徑進(jìn)行contigs擴(kuò)展，直到contigs擴(kuò)展結(jié)束。

2.1 De Bruijn圖的建立

建立De Bruijn圖時(shí)，為了節(jié)省存儲(chǔ)空間及便于信息的快速查找，分兩次解析數(shù)據(jù)文件，包括reads的堿基序列和堿基質(zhì)量數(shù)據(jù)。De Bruijn圖只存儲(chǔ)堿基序列，堿基質(zhì)量數(shù)據(jù)用來數(shù)據(jù)預(yù)處理。Reads上的連續(xù)子堿基串稱之為kmer，通常用K值表示kmer的長(zhǎng)度。

第一次解析reads初始數(shù)據(jù)文件：統(tǒng)計(jì)kmer數(shù)量，構(gòu)建無(wú)邊的De Bruijn圖。依次讀取reads數(shù)據(jù)文件中的每條reads上的kmer，通過哈希函數(shù)將kmer映射到De Bruijn圖中，并更新圖中kmer的數(shù)量。

第二次解析reads初始數(shù)據(jù)文件：存儲(chǔ)kmer堿基序列，并更新kmer位置、標(biāo)識(shí)等信息。依次解析數(shù)據(jù)文件中的每條reads， 并取得相應(yīng)的kmer堿基序列。對(duì)于每個(gè)kmer，首先計(jì)算其在De Bruijn圖中對(duì)應(yīng)的地址，然后將其堿基序列存儲(chǔ)在對(duì)應(yīng)的數(shù)組中。

2.2 De Bruijn圖的基本操作

序列拼接時(shí)，De Bruijn圖的主要操作包括kmer的查找、kmer的刪除、kmer的更新、kmer所在reads數(shù)據(jù)信息的更新等。kmer的堿基序列在De Bruijn圖中采用位存儲(chǔ)的方式存儲(chǔ)，其查找、刪除、更新等與位操作一一對(duì)應(yīng)。拼接時(shí)，需要將kmer的堿基序列經(jīng)過二進(jìn)制編碼轉(zhuǎn)換成整數(shù)，堿基序列所對(duì)應(yīng)的整數(shù)經(jīng)過位操作也可以轉(zhuǎn)換成kemr的堿基序列。位操作的方式一方面節(jié)省了kmer堿基序列存儲(chǔ)所需要的內(nèi)存空間，另一方面也提高了kmer的查找、刪除等操作的效率。

3 contigs拼接

拼接過程中，綜合考慮kmer的堿基質(zhì)量、拼接信息累計(jì)等各項(xiàng)信息，以整條reads為拼接單位，整體研究拼接問題，設(shè)計(jì)合理的評(píng)分機(jī)制，選擇得分最高的kmer進(jìn)行contigs擴(kuò)展，直至拼接結(jié)束，是基于reads引導(dǎo)的基因組序列拼接算法SRGA的主要思想。

3.1 kmer評(píng)分機(jī)制

Kmer的評(píng)分機(jī)制基于kmer的堿基質(zhì)量及所在reads的拼接信息累計(jì)。針對(duì)某特定的kmer，其得分計(jì)算方法為：首先計(jì)算拼接時(shí)當(dāng)前kmers所在的每條reads的得分，并將得分累加，即得該k-er的正向得分；然后取該kmer的反向互補(bǔ)序列，并計(jì)算其得分，兩個(gè)得分相加，即為該特定kmer的綜合得分。

將reads按kmer所出現(xiàn)的位置及reads的方向劃分成5個(gè)分值區(qū)域，其中的k-mers按照k-mer的位置和該read的方向劃分成5個(gè)區(qū)域，如圖2所示，“+”代表reads的正向得分區(qū)域，“-”代表reads的反向得分區(qū)域。

圖2 reads的得分區(qū)域示意

Fig.2 Score display of reads

kmer所在reads的得分計(jì)算公式為：

（2）

kmer綜合得分的計(jì)算公式為：

（3）

序列拼接過程中，一般會(huì)優(yōu)先選擇拼接即將完成的reads，但考慮到reads的數(shù)據(jù)區(qū)域特征，正向reads越靠近末端，其堿基的錯(cuò)誤率越高。為提高序列拼接的拼接質(zhì)量，將kmer拼接次數(shù)和堿基的錯(cuò)誤率結(jié)合起來，以評(píng)分的方法進(jìn)行量化，為kmer的導(dǎo)航選擇提供依據(jù)。

Kmer評(píng)分機(jī)制的數(shù)學(xué)模型如圖3所示，圖中陰影部分的面積即為kmer的得分。該得分綜合了reads的累計(jì)拼接信息和kmer所在reads區(qū)域的堿基質(zhì)量，reads拼接次數(shù)越多，堿基序列所在的區(qū)域分值越高，對(duì)kmer的得分貢獻(xiàn)越大。

圖3 kmer得分的數(shù)學(xué)模型

Fig.3 Mathematical model of kmer score

3.2 contigs的拼接

基于reads引導(dǎo)的contigs拼接過程如下：

step 1：選擇De Bruijn圖中出現(xiàn)次數(shù)為所有kmer出現(xiàn)次數(shù)平均值的kmer作為初始kmer，并將該kmer作為初始的contigs，開始第一輪拼接；

step 2：取得contigs擴(kuò)展待選的所有kmer，根據(jù)reads拼接信息及kmer所對(duì)應(yīng)的區(qū)域得分，對(duì)這些kmer進(jìn)行評(píng)分，選擇分值最高的kmer進(jìn)行拼接；

step 3：如果待選kmer為空，且拼接狀態(tài)處于第二輪拼接，則停止拼接，轉(zhuǎn)向step 6；如果仍處于第一輪拼接，則拼接狀態(tài)標(biāo)記為第二輪拼接，取contigs的反向互補(bǔ)，在其3端選擇kmer，進(jìn)行第二輪拼接；

Step 4：更新reads拼接信息，刪除拼接失敗或成功的reads信息；

Step 5：刪除De Bruijn圖中所有拼接成功reads中的kmer信息，并讀取reads的堿基序列，添加到contigs；

Step 6：拼接結(jié)束，保存contig及拼接成功的reads。

contigs分兩輪進(jìn)行，采用雙向擴(kuò)展。第一輪拼接從contig3端開始，直至擴(kuò)展終止；第二輪取contig的反向互補(bǔ)，其3端重新進(jìn)行擴(kuò)展，直至擴(kuò)展終止。contigs拼接過程如圖4所示。

圖4 contig雙向擴(kuò)展

Fig.4 Two directions extension of contigs

4 實(shí)驗(yàn)結(jié)果及分析

選取大腸桿菌和酵母菌的基因組作為測(cè)試數(shù)據(jù)，對(duì)SRGA進(jìn)行測(cè)試，將拼接結(jié)果與velvet[8]、Edena[9]和SOAPdenovo[10]等常用軟件進(jìn)行對(duì)比，并用Mauve Assembly Metrics軟件對(duì)各自的拼接結(jié)果進(jìn)行指標(biāo)分析、客觀評(píng)測(cè)。

選定的大腸桿菌E.coil（ERR001665）基因組數(shù)據(jù)，其數(shù)據(jù)量為4.2G ，reads長(zhǎng)度為36bp，目標(biāo)基因組序列為4.6M。

選定的酵母菌 S.pombe（ERR018885）基因組數(shù)據(jù)，其數(shù)據(jù)量為1.3G ，reads長(zhǎng)度為54bp，目標(biāo)基因組序列為12.6M。

兩組數(shù)據(jù)的測(cè)評(píng)結(jié)果分別如表1和表2所示。

表1 大腸桿菌的測(cè)評(píng)結(jié)果

Tab.1 Metrics of E.coil

如表1所示，針對(duì)大腸桿菌基因組，SRGA的Miscalled bases、N50、Missing bases、DCJ distance等指標(biāo)上優(yōu)于SOAPdenovo與velve，N90、N50、DCJ distance等優(yōu)于Edena。

表2 酵母菌的測(cè)評(píng)結(jié)果

Tab.2 Metrics of S.pombe

如表2所示，針對(duì)酵母菌基因組，SRGA的N90、N50、DCJ distance、Missing bases，等指標(biāo)Edena及velvet，Missing bases、DCJdistance、 Miscalled bases等優(yōu)于SOAPdenovo。

通過全面分析兩組基因組數(shù)據(jù)的測(cè)評(píng)結(jié)果可知，采用SRGA進(jìn)行的序列拼接，無(wú)論是拼接精度還是拼接的contigs長(zhǎng)度均表現(xiàn)不俗，且SRGA適用于大型基因組的序列拼接。

5 結(jié)束語(yǔ)

本文首次在基于reads引導(dǎo)的序列拼接算法中，提出了基于基因組序列數(shù)據(jù)特征和拼接信息累計(jì)的評(píng)分導(dǎo)航機(jī)制，避免了傳統(tǒng)的基于reads拼接的計(jì)算模型中數(shù)據(jù)質(zhì)量信息丟失的問題，為contigs的快速拼接提供了精確導(dǎo)航。然而，基因組序列拼接是一個(gè)NP難問題，需要不斷地創(chuàng)新和探索。本文的研究工作只是解決該難題的一小步，期待在此基礎(chǔ)有更大的進(jìn)展。

參考文獻(xiàn)：

[1] SANGER F， NICKLEN S， COULSON A R. DNA sequencing with chain-terminating inhibitors. 1977. [J]. Biotechnology， 1992， 24（12）：5463-5467.

[2] MAXAM A M， GILBERT W. A new method for sequencing DNA[J]. Proceedings of the National Academy of Sciences， 1977， 74（2）：560-564.

[3] POP M， SALZBERG S L. Bioinformatics challenges of new sequencing technology[J]. Trends Genet， 2008， 24（3）：142-149.

[4] PEVZNER P A， TANG H， WATERMAN M S. An Eulerian path approach to DNA fragment assembly[J]. Proceedings of the National Academy of Sciences of the United States of America， 2001， 98（17）：9748-9753.

[5] ERLICH Y， MITRA P P， Delabastide M， et al. Alta-Cyclic： a self-optimizing base caller for next-generation sequencing[J]. Nat Methods，2008， 5：679-682.

[6] SCHEIBYE-ALSING K， HOFFMANN S， FRANKEL A， et al. Sequence assembly[J]. Comput Biol Chem， 2009， 33（2）：121-136.

[7] DOHM J C， LOTTAZ C， BORODINA T， et al. Sharcgs， a fast and highly accurate short-read assembly algorithm for De Novo genomic sequencing[J]. Genome Research， 2007， 17（11）：1697-1706.

[8] ZERBINO D R， BIRNEY E. Velvet： Algorithms for De Novo short read assembly using De Bruijn graphs[J]. Genome Research， 2008， 18（5）：821-829.

[9] HERNANDEZ D， FRANCOIS P， FARINELLI L， et al. De Novo bacterial genome sequencing： millions of very short reads assembled on a desktop computer[J]. Genome Research 2008， 18（5）：802-809.

[10] LI R， ZHU H， RUAN J， et al. De Novo assembly of human genomes with massively parallel short read sequencing[J]. Genome Research， 2009， 20（2）：265-272.