基于de Bruijn圖的序列拼接算法研究與實現(xiàn)

李飛菲

（四川大學計算機學院，四川　610065）

基于de Bruijn圖的序列拼接算法研究與實現(xiàn)

李飛菲

（四川大學計算機學院，四川610065）

0　引言

全基因組測序是生物信息學最為重要的研究方向，基因組記錄了某個物種的全部遺傳信息，測定基因序列（DNA序列）可以破解其中的信息，進一步了解各類生物，并為人類健康和物種發(fā)展帶來實際意義。DNA序列由四種不同的堿基組成，即腺嘌呤 （A）、胞嘧啶（C）、鳥嘌呤（G）、胸腺嘧啶（T），它們按照特定的編碼規(guī)則連接成鏈?；驕y序就是為了測定DNA序列中堿基的排列順序。

目前的測序技術(shù)一次實驗測得的長度不超過1000個堿基對（bp）[1]，而普通高等生物基因中的堿基數(shù)目往往十分龐大，人類基因組的總長度就達到了30億bp。因此，測序時首先將DNA序列復(fù)制成多份，然后將長序列打斷成小序列片段（read），通過測序儀精確地測出每個read的堿基排序，然后利用序列拼接軟件將read在計算機上拼接起來，恢復(fù)成原來的DNA序列中的一條或多條連續(xù)段（contig）[2]。在實際操作過程中，DNA序列可能會出現(xiàn)丟失、變異、重復(fù)等問題，對如此復(fù)雜的序列進行拼接，DNA序列拼接算法的重要性不言而喻。

本文對基于de Bruijn圖的序列拼接算法進行研究，分析拼接過程中的每個階段，并利用C++語言編程實現(xiàn)該算法，對運行實驗結(jié)果進行分析，與Velvet拼接軟件的結(jié)果進行對比。

1　序列拼接算法的研究現(xiàn)狀

基于第一代測序技術(shù)的序列拼接方法適合處理長度為500-1000bp的read，利用read之間的重疊關(guān)系建立重疊圖，這種基于重疊圖的算法主要有三個階段：（1）overlap階段：將所有read進行兩兩比對，找到有重疊的read；（2）layout階段：利用重疊信息建立重疊圖；（3）consensus階段：遍歷重疊圖，圖中的每條路徑都可能成為contig，將得到的所有contig進行組裝，得到最后的拼接結(jié)果。

隨著生物信息的飛速增長，生物組織中的多樣性和關(guān)聯(lián)性不斷增強，傳統(tǒng)的測序方法已經(jīng)不能滿足當下的需求。第二代測序技術(shù)面對的數(shù)據(jù)長度極短、數(shù)量龐大、覆蓋度較高，最短的只有25-50bp，不適合建立重疊圖?；赿e Bruijn圖的序列拼接算法是新一代測序技術(shù)最常用的方法，適用于拼接高通量的短read，已經(jīng)被應(yīng)用于各種序列拼接技術(shù)，如Velvet[3]、SOAPdenovo[4]、AbySS[5]等。它將read轉(zhuǎn)化成一組連續(xù)的k-mer，然后尋找kmer之間的重疊關(guān)系，建立de Bruijn圖，最后將問題轉(zhuǎn)化為尋找歐拉路徑問題。

2　de Bruijn圖簡介

假設(shè)給定的DNA序列集合是S={S1，S2,S3，…，Sn}，其中Si（1≤i≤n）表示一個長度為n的 read。對每個read，從頭開始依次取出長度為k（1≤k≤n）的子串，該子串稱為k-mer，一個read可以得到（n-k+1）個連續(xù)的k-mer，它的集合為K={k1，k2，…，kn-k+1}。

在de Bruijn圖結(jié)構(gòu)中，以k-mer作為節(jié)點，如果存在兩個k-mer在read中相鄰，即k1的后（k-1）個堿基與k2的前（k-1）個堿基重疊，那么圖中對應(yīng)的兩個節(jié)點之間存在一條由k1指向k2的邊。由于k-mer是根據(jù)read產(chǎn)生的，因此每一條read可以看作圖中的一條路徑。由此可以構(gòu)建出一個以k-mer為基本單位的有向圖，如圖1所示。

圖1　de Bruijn圖的構(gòu)建過程

這樣的結(jié)構(gòu)可以保證每個k-mer在圖中只出現(xiàn)一次，read之間所有的連接關(guān)系都被轉(zhuǎn)化到了k-mer上。根據(jù)k-mer之間的重疊信息，尋找一條經(jīng)過所有邊一次且僅一次的路徑，即可得到一個近似原DNA序列的拼接結(jié)果。

序列在處理過程中可能產(chǎn)生丟失、變異、重復(fù)等問題，會使de Bruijn圖產(chǎn)生三種結(jié)構(gòu)：Tips結(jié)構(gòu)、Bubbles結(jié)構(gòu)和Repeats結(jié)構(gòu)。

Tips是因為read在邊緣出現(xiàn)了錯誤，導(dǎo)致圖中的路徑在某一端產(chǎn)生分叉，這種結(jié)構(gòu)會影響到圖的遍歷。如圖2所示，當走到GAA節(jié)點時，無法確定接下來該去往AAT還是AAA。

Bubbles是因為read中間出現(xiàn)了錯誤，導(dǎo)致圖中路徑的開始節(jié)點和結(jié)束節(jié)點相同，但是中間的節(jié)點序列不完全相同，這種結(jié)構(gòu)同樣會影響遍歷結(jié)果。如圖3所示，當走到ATC節(jié)點時，無法選擇出正確的路徑。

圖2　Tips結(jié)構(gòu)

圖3　Bubbles結(jié)構(gòu)

Repeats其實不算是一種錯誤結(jié)構(gòu)，它是由序列片段內(nèi)部的重復(fù)產(chǎn)生的，會使de Bruijn圖產(chǎn)生循環(huán)結(jié)構(gòu)，無法確定遍歷方向。如圖4所示，當走到ATC節(jié)點時，由于重復(fù)區(qū)域TCG、CGA的存在，無法直接確定應(yīng)該去往GAA還是GAT，走到TTC節(jié)點時也有同樣的困擾。

圖4　Repeats結(jié)構(gòu)

這些結(jié)構(gòu)會影響de Bruijn圖的遍歷，導(dǎo)致無法找到正確的路徑，在實現(xiàn)拼接算法時應(yīng)該進行相應(yīng)的處理。

3　算法實現(xiàn)

基于de Bruijn圖的序列拼接算法可以分為三個階段：構(gòu)建圖、化簡圖、消除錯誤結(jié)構(gòu)。

首先逐條處理read，按照指定的k-mer長度將read劃分為一系列k-mer。為了保證序列拼接的完整性，對于每個k-mer求出其反向互補的k-mer’，比較兩個k-mer的大小，將大的k-mer稱為k-mer+。分別將兩個k-mer存儲到圖中，以k-mer+所在節(jié)點為代表，任何對某節(jié)點的處理都會影響到它的反向互補節(jié)點（TwinNode）。然后根據(jù)該k-mer在read中與其他kmer的重疊關(guān)系，在節(jié)點之間建立有向邊，如果節(jié)點A到節(jié)點B有一條有向邊，那么反向互補節(jié)點B’同樣應(yīng)該有一條有向邊指向A’。與此同時，還需要累計每個節(jié)點和每條邊出現(xiàn)的次數(shù)，作為該節(jié)點或邊的覆蓋率，覆蓋率越低表明這個節(jié)點越不可靠。

本文設(shè)計了三個數(shù)據(jù)類型，分別是節(jié)點Node，邊Arc和圖Graph。Node中存儲了k-mer序列、覆蓋率、邊等信息，并且有一個指向TwinNode的指針；Arc中存儲了節(jié)點的出邊，用指針來指向下一個節(jié)點；Graph中含有一個以k-mer為索引的哈希表，存儲了相應(yīng)的Node信息。三者之間的關(guān)系如圖5所示，從圖中已經(jīng)無法看出read之間的重疊關(guān)系，兩個節(jié)點只通過一條邊連接，這種結(jié)構(gòu)消耗的內(nèi)存少，而且搜索便捷。

圖5　de Bruijn圖的結(jié)構(gòu)關(guān)系

經(jīng)過對de Bruijn圖的構(gòu)建可以發(fā)現(xiàn)，圖中的節(jié)點非常多，每個k-mer都是單獨的一個節(jié)點，想要直接發(fā)現(xiàn)其中的關(guān)聯(lián)非常困難，因此可以對圖進行化簡，從而使圖中的結(jié)構(gòu)更加清晰。當節(jié)點A只有一個出邊指向節(jié)點B，并且B也只有一個來自A的入邊，那么就可以把A和B兩個節(jié)點合并成一個節(jié)點，同時合并它們的反向互補節(jié)點A’和B’。值得注意的是，節(jié)點與邊的覆蓋率也應(yīng)該進行相應(yīng)的合并。

化簡圖操作不會影響節(jié)點的信息，同時能夠有效地減少內(nèi)存占用。為了保證de Bruijn圖一直保持最簡，在消除錯誤結(jié)構(gòu)時需要迭代的對圖化簡。

（1）消除Tips結(jié)構(gòu)

首先設(shè)置一個大小為2kb的閾值，當Tips的長度小于閾值時才可以消除它。這個閾值遠大于read的長度，因為Tips中包含的read越多，它的長度就越長，出錯的可能性就越小。然后比較Tips節(jié)點與其他同源“兄弟節(jié)點”的覆蓋率大小，如果沒有比Tips的覆蓋率更小的節(jié)點，那么說明其他路徑都更加的普遍，是正確路徑的可能性更高。如果以上兩個條件都滿足，直接刪除該Tips節(jié)點即可。

（2）消除Bubbles結(jié)構(gòu)

首先根據(jù)Bubbles結(jié)構(gòu)的特點，使用Dijsktra廣度優(yōu)先算法，從任意節(jié)點出發(fā)，訪問它所能達到的其他節(jié)點，優(yōu)先訪問快的那條邊。從節(jié)點A到節(jié)點B所需的時間可以根據(jù)節(jié)點的長度和邊的覆蓋率計算得到，如公式1所示。

對訪問過的節(jié)點進行標記，當遇到了已標記的節(jié)點時開始回溯，找出兩條路徑最近的共同祖先節(jié)點，記錄下這兩條路徑。將這兩條路徑進行對齊處理，如果足夠相似就將它們合并成一條，否則選擇快的那條路徑，因為消耗時間少的邊覆蓋率高，所以是正確路徑的可能性更高。

（3）消除Repeats結(jié)構(gòu)

首先從Repeats結(jié)構(gòu)中的任意節(jié)點出發(fā)，確定該節(jié)點所在的read，將read中的節(jié)點進行標記，然后根據(jù)不同read在標記節(jié)點上的重疊情況擴展到下一個read，繼續(xù)對后續(xù)節(jié)點進行標記，直到這個Repeats結(jié)構(gòu)被處理完。最后將被標記的節(jié)點串連起來，可以得到一條確定的路徑，消除Repeat結(jié)構(gòu)產(chǎn)生的困擾。如圖6所示，重復(fù)區(qū)域是節(jié)點C、D，首先從節(jié)點A出發(fā)進行標記，由于read1能夠覆蓋A和C兩個節(jié)點，因此對C進行標記；然后從節(jié)點C出發(fā)，可以發(fā)現(xiàn)C還存在于read2中，而且read2覆蓋了C、D、F三個節(jié)點，按照read2的路徑繼續(xù)進行擴展和標記；最后將被標記的節(jié)點串連起來，可以確定該Repeats結(jié)構(gòu)的一種尋路方式為A→C→D→F。

圖6　Repeats結(jié)構(gòu)的解決方法

在消除以上幾種錯誤結(jié)構(gòu)后，可能還存在一些無法從圖中直接識別出來的連接錯誤。設(shè)置一個覆蓋率閾值，將小于這個閾值的節(jié)點全部刪除，消除錯誤連接的影響。這個閾值的大小可以根據(jù)實際情況自定義。

以上全部操作完成后，再對圖進行一次化簡合并，可以得到若干條較長的序列片段，稱為contig。一個contig包含多個k-mer甚至多個read。這些contig就是最終的拼接結(jié)果，根據(jù)實際需要還可以對它們做進一步的拼接處理，得到一個含有若干空隙（gap）的長序列，盡可能的還原到原始DNA序列。

4　結(jié)果分析

實驗環(huán)境為處理器：Intel Core i3-3220 CPU@ 3.30GHz；內(nèi)存：4.0GB；操作系統(tǒng)：Windows 7。

本文選取了三種不同規(guī)模的DNA序列模擬樣本進行實驗。Read長度為35，具體信息如表1所示。

表1　實驗數(shù)據(jù)集信息

本文以contig片段在原始序列上的恢復(fù)率為衡量標準，判斷算法的正確性與可行性。假設(shè)原始DNA序列長度為L，拼接算法得到的contig個數(shù)為n，每個contig的長度為ci（1≤i≤n），其中拼接錯誤的contig總長度為X，那么算法的恢復(fù)率如公式2所示。

由于gap的存在，contig的恢復(fù)率很難達到100%，為了有效地驗證實驗結(jié)果，本文使用Velvet算法運行相同的實驗數(shù)據(jù)，將得到的結(jié)果進行對比分析。K值為21，實驗結(jié)果如表2所示。

表2　實驗結(jié)果

從當前的實驗數(shù)據(jù)來看，本文算法得到的contig長度較長、正確率較高，與Velvet的恢復(fù)率相差不大，可以證明本文中的算法是正確的、可行的。

5　結(jié)語

本文基于de Bruijn圖，設(shè)計了一種數(shù)據(jù)結(jié)構(gòu)，通過化簡、消除錯誤結(jié)構(gòu)等操作對短基因序列進行拼接，并且得到了與原始序列很相似的拼接結(jié)果。但本文也存在一些不足之處，對錯誤結(jié)構(gòu)的處理不夠全面，在一定程度上影響了最終尋路的正確性，導(dǎo)致無法獲得更長更連續(xù)的contig，因此在消除錯誤結(jié)構(gòu)的算法上還需要繼續(xù)研究和改善。

[1]駱志剛等.DNA序列拼接的研究和挑戰(zhàn)[J].計算機工程與科學,2007.

[2]王旭.基于deBruijn圖的DNAcontig生成算法[D].哈爾濱工業(yè)大學:王旭,2011.

[3]Daniel R.Zerbino,Ewan Birney.Velvet Algorithms for de Novo Short Read Assembly Using de Bruijn Graphs[J].Genome Res,2008.

[4]Li R,Zhu H,Ruan J,et al.De Novo Assembly of Human Genomes with Massively Parallel Short Read Sequencing.GenomeRes，2010.

[5]Simpson JT,Wong K,Jackman SD,et al.ABYSS:A Parallel Assembler for Short Read Sequence Data[J].Genome Res,2009.

Next-generation Sequencing Technology;High-throughput Sequencing;DNA Sequence Assembly;de Bruijn Graph

Research&Implementation of Sequence Assembly Algorithm Based on de Bruijn Graphs

LI Fei-fei
（College of Computer Science,Sichuan University，Sichuan 610065）

1007-1423（2016）02-0003-05

10.3969/j.issn.1007-1423.2016.02.001

李飛菲（1991-），女，北京人，本科，研究方向為計算機圖形學、虛擬現(xiàn)實

2015-11-24

2015-12-28

序列拼接算法是DNA測序過程中的關(guān)鍵技術(shù)。隨著新一代測序技術(shù)的發(fā)展，如何實現(xiàn)高通量、高效率測序已經(jīng)成為生物信息學領(lǐng)域的重要挑戰(zhàn)，序列拼接算法也在逐漸改進以提高拼接效果。基于de Bruijn圖的序列拼接算法是目前使用最廣泛的方法之一，對其進行分析研究，利用C++編程實現(xiàn)該算法，并對實驗結(jié)果進行分析。

新一代測序技術(shù)；高通量測序；基因拼接；de Bruijn圖

Algorithms for DNA sequence assembly are the key process of gene sequencing.With the development of next-generation sequencing technologies,how to achieve high-throughput,high-efficiency sequencing has become an important challenge in bioinformatics.At the same time,sequence assembly algorithms also have been improved to get better results.The algorithms based on de Bruijn graph are the most popular DNA sequence assembly algorithm,analyses and researches on it,applies C++language to program this algorithm,and analyses the experiment results.