基于de Bruijn圖的序列拼接算法研究與實現

李飛菲

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

基于de Bruijn圖的序列拼接算法研究與實現

李飛菲

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

0　引言

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

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

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

1　序列拼接算法的研究現狀

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

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

2　de Bruijn圖簡介

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

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

圖1　de Bruijn圖的構建過程

這樣的結構可以保證每個k-mer在圖中只出現一次，read之間所有的連接關系都被轉化到了k-mer上。根據k-mer之間的重疊信息，尋找一條經過所有邊一次且僅一次的路徑，即可得到一個近似原DNA序列的拼接結果。

序列在處理過程中可能產生丟失、變異、重復等問題，會使de Bruijn圖產生三種結構：Tips結構、Bubbles結構和Repeats結構。

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

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

圖2　Tips結構

圖3　Bubbles結構

Repeats其實不算是一種錯誤結構，它是由序列片段內部的重復產生的，會使de Bruijn圖產生循環結構，無法確定遍歷方向。如圖4所示，當走到ATC節點時，由于重復區域TCG、CGA的存在，無法直接確定應該去往GAA還是GAT，走到TTC節點時也有同樣的困擾。

圖4　Repeats結構

這些結構會影響de Bruijn圖的遍歷，導致無法找到正確的路徑，在實現拼接算法時應該進行相應的處理。

3　算法實現

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

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

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

圖5　de Bruijn圖的結構關系

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

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

（1）消除Tips結構

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

（2）消除Bubbles結構

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

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

（3）消除Repeats結構

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

圖6　Repeats結構的解決方法

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

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

4　結果分析

實驗環境為處理器：Intel Core i3-3220 CPU@ 3.30GHz；內存：4.0GB；操作系統：Windows 7。

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

表1　實驗數據集信息

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

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

表2　實驗結果

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

5　結語

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

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

[2]王旭.基于deBruijn圖的DNAcontig生成算法[D].哈爾濱工業大學:王旭,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-），女，北京人，本科，研究方向為計算機圖形學、虛擬現實

2015-11-24

2015-12-28

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

新一代測序技術；高通量測序；基因拼接；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.