張磊

摘要：研究基因組拼接算法，進(jìn)而更高效地實(shí)現(xiàn)全基因組拼接、獲取生物體遺傳信息，對(duì)于生命科學(xué)研究具有重要的意義。 拼接之前對(duì)數(shù)據(jù)進(jìn)行一定的預(yù)處理。拼接采用經(jīng)典的DBG算法，分析了DBG算法的實(shí)現(xiàn)步驟以及需要注意的問題。根據(jù)read切割K-mer，然后采用位運(yùn)算壓縮存入Hash表。以有效的K-mer為頂點(diǎn)，相鄰的K-mer之間連有向邊，建立DBG圖，然后在該圖中尋找Euler路徑。每一條Euler路徑對(duì)應(yīng)一條contig。最后進(jìn)行參數(shù)組合實(shí)驗(yàn)，并將結(jié)果與專業(yè)軟件CLC的結(jié)果進(jìn)行對(duì)比。兩者效率相近，說明了算法的高效性。

關(guān)鍵詞：拼接；DBG圖；K-mer；contig；scaffold

中圖分類號(hào)：Q812；TP311 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1009-3044（2015）33-0123-02

快速和準(zhǔn)確地獲取生物體的遺傳信息對(duì)于生命科學(xué)研究具有重要的意義。對(duì)每個(gè)生物體來說，基因組包含了整個(gè)生物體的遺傳信息，這些信息通常由組成基因組的DNA或RNA分子中堿基對(duì)的排列順序所決定。獲得目標(biāo)生物基因組的序列信息，進(jìn)而比較全面地揭示基因組的復(fù)雜性和多樣性，成為生命科學(xué)領(lǐng)域的重要研究內(nèi)容。

本文深入探討了基因組拼接算法的實(shí)現(xiàn)細(xì)節(jié)，為實(shí)現(xiàn)基因組組裝算法小平臺(tái)化提供了參考。

測試數(shù)據(jù)：一個(gè)全長約為120000個(gè)堿基對(duì)的細(xì)菌人工染色體（BAC）， 采用Hiseq2000測序儀進(jìn)行測序，測序深度（sequencing depth）約為70×。

1 背景介紹

基因組組裝是指以基因組測序所得的序列片段read數(shù)據(jù)為輸入，利用序列片段之間的交疊關(guān)系，再加上配對(duì)數(shù)據(jù)、堿基質(zhì)量等輔助信息，通過一定的方法手段重建目標(biāo)基因組序列的過程?；谛乱淮鷾y序數(shù)據(jù)的拼接算法，通常包括如下三個(gè)階段：

1） 數(shù)據(jù)的預(yù)處理：主要是對(duì)測序數(shù)據(jù)中的堿基進(jìn)行糾錯(cuò)；

2） 拼接contig：短read片段被拼接成長度較長的contig；

3） 組裝scaffold：利用配對(duì)數(shù)據(jù)，確定contig之間的相對(duì)位置關(guān)系和距離，對(duì)contig進(jìn)行組裝。

2 拼接算法

本文討論的是第二個(gè)階段——拼接contig。

2.1 DBG算法簡介

目前，基于DBG圖策略的拼接算法被最廣泛應(yīng)用到新一代測序數(shù)據(jù)的處理中，比如華大基因的SOAPdenovo [5]。基于DBG圖的拼接算法，非常巧妙地將具有交疊關(guān)系的read映射到一起，降低了計(jì)算交疊時(shí)的復(fù)雜度，減少了內(nèi)存消耗。本文也采用DBG算法進(jìn)行拼接。

在本文中，主要用到的經(jīng)典算法有Euler路徑套圈算法，以及Hash表的插入算法。讀者可以參考相關(guān)的算法書籍。

DBG圖以K-mer為頂點(diǎn)，相鄰的K-mer之間連有向邊，配對(duì)文件中的任一條read對(duì)應(yīng)一個(gè)K-mer路徑，將尋找contig的問題轉(zhuǎn)化為在該圖中尋找一條Euler路徑。

2.2 拼接算法大致流程

筆者實(shí)現(xiàn)的DBG算法，具體步驟如下：

S1：用戶輸入?yún)?shù)：K-mer的大小、K-mer有效的最低頻次、contig有效的最短長度。

S2：讀入并保存預(yù)處理后的文件中的read，并生成對(duì)應(yīng)的反向互補(bǔ)串。

S3：Hash表初始化。根據(jù)指定的K-mer，將所有的read及其反向互補(bǔ)串連續(xù)地切割成指定長度的串，相鄰的兩個(gè)K-mer交疊K-1個(gè)堿基。利用位運(yùn)算壓縮存儲(chǔ)，插入Hash表。數(shù)組p[c][i][j]記錄該K-mer對(duì)應(yīng)的去重后編號(hào)。

S4：根據(jù)數(shù)組p[c][i][j]構(gòu)建數(shù)組opp[]。

S5：將Hash表去除冗余，存入數(shù)組kmer[]。

S6：根據(jù)數(shù)組p[c][i][j]構(gòu)建DBG圖。如果去重后編號(hào)為p[c][i][j]和p[c][i][j+1]的K-mer都是有效的，則在兩個(gè)K-mer對(duì)應(yīng)的頂點(diǎn)之間連一條有向邊，由p[c][i][j]指向p[c][i][j+1]，記錄在數(shù)組map[][]里面。同時(shí)記錄每個(gè)頂點(diǎn)i的入度數(shù)組ind[]和出度數(shù)組outd[]。

S7：從每個(gè)有效而且ind[i]=0，outd[i]>0的頂點(diǎn)出發(fā)，尋找Euler路徑，每一條Euler路徑對(duì)應(yīng)一個(gè)contig。尋找的時(shí)候，總是出現(xiàn)頻次最大的后繼K-mer對(duì)應(yīng)的路徑優(yōu)先走。每經(jīng)過一條邊，便在DBG圖中刪除該邊以及相應(yīng)的反向互補(bǔ)串對(duì)應(yīng)的邊。每找完一個(gè)contig，記錄其長度，并根據(jù)contig對(duì)應(yīng)的K-mer序列，輸出該contig。

S8：根據(jù)各個(gè)contig的長度，并計(jì)算N50以及拼接序列總長度。

2.3 拼接算法的細(xì)節(jié)處理

關(guān)于DBG算法，需要說明以下幾點(diǎn)：

1）S2中生成反向互補(bǔ)串，這是由DNA實(shí)際結(jié)構(gòu)決定的：DNA的兩條鏈之間的關(guān)系為反向互補(bǔ)。

2）S3中位運(yùn)算采用二進(jìn)制編碼，堿基與堿基編碼對(duì)應(yīng)關(guān)系為：A-00，C-01，T-10，G-11。采用二進(jìn)制編碼，將存儲(chǔ)字符轉(zhuǎn)化為存儲(chǔ)二進(jìn)制數(shù)，大大減少了內(nèi)存的消耗。一般來說，對(duì)整數(shù)的位操作遠(yuǎn)遠(yuǎn)快于對(duì)序列堿基的字符操作。其中，K-mer的第k個(gè)堿基對(duì)應(yīng)的int的下標(biāo)p1計(jì)算如下：

[p1=ints-[（K-mer-k+1）×2-1]/bits] 其中：ints為用來存儲(chǔ)K-mer的int數(shù)據(jù)的個(gè)數(shù)，bits為一個(gè)int用來存儲(chǔ)堿基信息的位數(shù)（int數(shù)據(jù)的最高位為符號(hào)位，所以bits最多為30）。式中除號(hào)為整除。

K-mer的第k個(gè)堿基對(duì)應(yīng)int的數(shù)據(jù)位p2，p2+1，p2計(jì)算如下：

[p2=bits-[（K-mer-k+1）×2-1]%bits]其中：%表示求余運(yùn)算。

3）S3中采用的Hash函數(shù)為：

[p=（i=1intsbit[i]/ints）%hashsize]

其中：bit[i]為該K-mer對(duì)應(yīng)的第i個(gè)int數(shù)據(jù)，hashsize為Hash表的大小。

4）經(jīng)過實(shí)驗(yàn)，采用上述Hash函數(shù)時(shí)，當(dāng)Hash表的大小為最大可能值的3-4倍時(shí)，沖突便可以控制在理想的范圍以內(nèi)。

5）S3中p[c][i][j]，c=1表示原串，c=2表示反向互補(bǔ)串，i表示read序號(hào)，j表示K-mer在read中第一堿基的位置。該數(shù)組是為了計(jì)算數(shù)組opp[]和構(gòu)建DBG圖。

6）S4中opp[i]表示去重后編號(hào)為i的K-mer的反向互補(bǔ)串對(duì)應(yīng)的去重后編號(hào)。

7）S6中map[i][j]，i表示某個(gè)K-mer的去重后編號(hào)，j為某個(gè)堿基對(duì)應(yīng)的二進(jìn)制編碼。若map[i][j]=0，表示去重后編號(hào)為i的K-mer后面不可能緊鄰j對(duì)應(yīng)的堿基。若map[i][j]>0，則表示去重后編號(hào)為i的K-mer去掉第一個(gè)堿基，后面緊鄰j對(duì)應(yīng)的堿基形成的K-mer的去重后編號(hào)。

8）S7中條件ind[i]=0和outd[i]>0保證了i一定是某個(gè)read（包括反向互補(bǔ)串）的初始K-mer。

9）S7中選擇出現(xiàn)頻次最大的后繼K-mer對(duì)應(yīng)的路徑優(yōu)先走，是因?yàn)槌霈F(xiàn)頻次越大，該K-mer正確的可能性越大。這樣可以有效地避免遇到錯(cuò)誤堿基而導(dǎo)致的拼接中斷。

10）S7中如果不刪除相應(yīng)的反向互補(bǔ)串對(duì)應(yīng)的邊，可能會(huì)出現(xiàn)兩個(gè)contig或者其子串之間存在反向互補(bǔ)關(guān)系的現(xiàn)象，而這樣顯然是冗余的。

11）本文中未考慮重復(fù)片段的影響。位于重復(fù)區(qū)域的read會(huì)被集中拼接到一個(gè)區(qū)域，而在另一個(gè)區(qū)域卻無法拼接。

12）基于DBG圖的拼接算法的一個(gè)關(guān)鍵參數(shù)是K值，一般認(rèn)為對(duì)于原核生物的基因組拼接，K值在21-35之間是合適的[6-7]。提供的測試數(shù)據(jù)為某細(xì)菌的，而細(xì)菌屬于原核生物。

4 結(jié)論

目前的幾乎所有組裝軟件都需要較大內(nèi)存的計(jì)算平臺(tái)，但是本文的程序在普通的PC機(jī)上即可實(shí)現(xiàn)，具有很大的可推廣價(jià)值。但是，對(duì)測序錯(cuò)誤導(dǎo)致的tip、bubble結(jié)構(gòu)以及repeat的處理，沒有進(jìn)行深入考慮。另一方面，需要用戶進(jìn)行參數(shù)組合實(shí)驗(yàn)，不能自動(dòng)計(jì)算最優(yōu)參數(shù)組合。如果進(jìn)一步進(jìn)行研究，可以從采用多種配對(duì)文庫、長序列與短序列混合、并行化拼接等方面進(jìn)行考慮。

參考文獻(xiàn)：

[1] Lucian I， Farideh F， Silvana I. HiTEC： accurate error correction in high-throughput sequencing data[J/OL]. （2010-07-28）[2012-05-28].

[2] Jan S， Heiko S， Simon J P， et al. SHREC： a short-read error correction method[EB/OL]. Oxford：OxfordJournal，2009[2012-06-03].http：//bioinformatics.oxfordjournals.org/content.

[3] Yang X， Dorman K S， Aluru S. Reptile： representative tiling for short read error correction[EB /OL]. Bioinformatics， 2010， 285： 313-314[2012-06-10].

[4] Hernandez D， Francois P， Farinelli L. De Novo Bacterial Genome Sequencing： Millions of Very Short Reads Assembled on a Desktop Computer[J]. Genome Research 2008， 18（5）：802-809.

[5] Zerbino DR， Birney E. Velvet： Algorithms for De Novo Short Read Assembly Using De Bruijn Graphs[J]. Genome Research 2008， 18（5）：821-829.

[6] 趙瑋. 基于配對(duì)信息的 DNA 從頭測序組裝算法研究[D]. 哈爾濱： 哈爾濱工業(yè)大學(xué)， 2011.

[7] 韓東濤. 基于概率模型的基因組從頭測序算法研究[D]. 哈爾濱： 哈爾濱工業(yè)大學(xué)， 2012.

[8] 曾培龍. 基于 reads 引導(dǎo)的基因組序列拼接[D]. 哈爾濱： 哈爾濱工業(yè)大學(xué)， 2012.

[9] 黃勇. 基于高通量測序的微生物基因組學(xué)研究[D]. 北京： 中國人民解放軍軍事醫(yī)學(xué)科學(xué)院， 2013.

[10] 王東陽， 任世軍， 王亞東. DNA 序列拼接中 de Bruijn 圖結(jié)構(gòu)的研究[J]. 智能計(jì)算機(jī)與應(yīng)用， 2011， 1（2）： 20-25.