買阿麗， 楊雯雯

(運城學(xué)院 應(yīng)用數(shù)學(xué)系 山西 運城 044000)

關(guān)于基因重組中OLC算法的改進研究

買阿麗， 楊雯雯

(運城學(xué)院 應(yīng)用數(shù)學(xué)系 山西 運城 044000)

針對基因組組裝問題，從數(shù)據(jù)預(yù)處理，利用KMP算法在O(m+n)的時間上快速確定某兩個堿基片段的最大重復(fù)度，將讀長序列依據(jù)Overlap圖連成Contigs鏈以及ContigsN50的確定4個環(huán)節(jié)，改進現(xiàn)有的OLC拼接技術(shù)，并給出優(yōu)化后的模型和算法，較好地解決了基因組組裝問題.

基因重組； KMP算法； Overlap圖； Contigs鏈； ContigsN50

0 引言

確定基因組堿基對序列的過程，稱為測序.目前，測序技術(shù)正向著高通量、低成本的方向發(fā)展[1-3].但能直接讀取的堿基對序列長度遠小于基因組序列長度，因此需要利用一定的方法將測序得到的短片段序列組裝成更長的序列.

常用的組裝算法主要基于OLC(Overlap/Layout/Consensus)方法、貪婪圖方法及De Bruijn圖方法等[1-2,4].一個好的算法應(yīng)具備組裝效果好、時間短和內(nèi)存小等特點[5-7].新一代測序技術(shù)的性能還有較大的改善空間，曾培龍分析基因組序列拼接所面臨的主要挑戰(zhàn)，比如大量重復(fù)片段的存在，reads數(shù)據(jù)海量、長度短及含有測序錯誤等，探討了當前基因組序列拼接所采用的貪心策略、交疊-排序-生成共有序列(OLC)策略和De Bruijn圖策略等，總結(jié)了不同算法的優(yōu)勢及不足，提出了序列拼接算法的改進方向，率先在拼接算法中提出了基于信息累計和數(shù)據(jù)特征相結(jié)合的評分方法[4].鄭緯名等提出了歐拉超路拼接算法，由于該算法要求構(gòu)造一個復(fù)雜的De Bruijn圖,因此用歐拉超路算法拼接大規(guī)模全基因組存在存儲瓶頸問題[8].王旭提出了一種新的重疊群生成算法SRGA,該算法基于De Bruijn圖,將從頭測序問題轉(zhuǎn)化成一個四叉樹的搜索問題,并采用啟發(fā)式搜索策略,能夠快速地處理海量測序數(shù)據(jù),能得到質(zhì)量較高的重疊群[9].

基因組組裝就是將得到的有ATCG4個字符組成的字符串，拼接成更長的字符串，本質(zhì)上就是字符串匹配問題，利用KMP算法可以較快速地實現(xiàn)基因重組.KMP算法如圖1所示[10].本文以2014年深圳杯全國大學(xué)生數(shù)學(xué)建模B題(以下簡稱2014B題)附件為例在現(xiàn)有的OLC拼接技術(shù)上提出改進.2014B題中數(shù)據(jù)為一個全長約為120 000個堿基對的細菌人工染色體(BAC)，采用Hiseq2000測序儀進行測序，基因組中每個位置平均被測到約70次[11].

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

2014B題得到的長數(shù)據(jù)存儲在fastq文件中，首先將數(shù)據(jù)讀取到Matlab中，并依據(jù)以下原則對數(shù)據(jù)進行篩選，去除不準確的數(shù)據(jù).

1) 整條read堿基中無N堿基.N代表模糊堿基，可能是由于測序熒光強度不夠造成，將其刪除[12].

2) 整條read中質(zhì)量值低于20的堿基不超過20%.根據(jù)文獻[12]中給出的公式：字符的ASCII碼-64=質(zhì)量值.也即將read_zhiliang矩陣中ASCII碼值中小于84的數(shù)據(jù)占該行比例小于20%的行刪除，對應(yīng)將reads1、reads2中的該堿基片段刪除.

3) 整條read中質(zhì)量值低于13的堿基不超過10 %.也即將read_zhiliang矩陣中ASCII碼值中小于77的數(shù)據(jù)占該行比例小于10 %的行刪除，對應(yīng)將reads1、reads2中的該堿基片段刪除[12].

4) 整條read的堿基平均值不小于20.將不符合的片段刪除[12].

依據(jù)刪除數(shù)據(jù)的4個標準，reads1共刪掉2 805行，reads2共刪掉4 487行.得到43 982×89的reads1_used矩陣和42 360×89的reads2_used矩陣，其中第89列記錄該行在reads1、reads2中的位置.

圖1 KMP算法Fig.1 KMP algorithm

2 建立Overlap有向圖

若有兩個短片段序列分別為AGTACCTTGCTAGCGT 和GCTAGCGTAGGTCTGA.則有可能基因組序列中包含有AGTACCTTGCTAGCGTAGGTCTGA這一段.事實上，一個片段可能與多個片段都有重復(fù)部分，但是重復(fù)的堿基數(shù)量不同，記此數(shù)量為重復(fù)度.認為重復(fù)越大，這兩個片段越有可能連在一起.因此，找到reads1中每個片段與之重復(fù)度最大的片段，記下此時的重復(fù)度，命為最大重復(fù)度.整理這些最大重復(fù)度，建立一個Overlap有向圖.

2.1 計算兩個片段最大重復(fù)度

若逐一比較整段基因片段，則需將兩個基因片段從不同的位置開始依次比較，這個過程比較復(fù)雜，效率很低.因此，在比較兩段基因之前先進行查找，快速鎖定可能的重復(fù)位置，根據(jù)可能的重復(fù)位置再從最大的可能重復(fù)的位置開始比較.以片段a：TCGATACTAG 和片段b：ACTAGGCTAG為例說明.

1) 在b片段中利用KMP算法查找a片段最后3個字符串TAG，稱為子串，并記下子串在b中出現(xiàn)的每個起始地址.

3 7

2) 先從位置7開始，倒序向前逐一比較每個字符.

倒序第5位對應(yīng)不相等，于是再從位置3開始比較.

此時重疊部分相同，重疊部分有5個堿基，認為a、b的最大重復(fù)度為5. 稱此時子串的起始地址為重復(fù)地址，即a與b重復(fù)地址為3.

2.2 在多個片段尋找某個片段的后繼片段

以a：TCGATACTAG; b：ACTAGGCTAG; c：TAGCCAATTA; d：GATACTAGAC4個片段為例，尋找a的后繼片段：

3 7

6

1) 用上述1中的方法得出b與a的重復(fù)地址為3，最大重復(fù)度為5.

2) 在c片段大于等于3的位置上利用KMP算法尋找“TAG”，不是從c的第一個位置開始找，這樣就減少了查找次數(shù)，提高了執(zhí)行效率.此例中，發(fā)現(xiàn)沒有找到.此時仍然認為b最有可能是a的后繼片段.

3) 在d片段大于等于3的位置上利用KMP算法尋找“TAG”，發(fā)現(xiàn)位置6出現(xiàn)“TAG”,再用2.1中的方法比較d與a，得出重復(fù)度為8，大于b與a的重復(fù)度，于是認為在b、c、d中d最可能為a的后繼片段.

2.3 建立Overlap圖算法

Overlap圖模型算法1：等長reads重復(fù)度計算.

算法描述：對每個ri片段在剩余的片段中尋找后繼片段，并計算最大重復(fù)度.

算法輸入：(1) reads片段集合,即ri集合； (2) 最小重復(fù)度d.

算法輸出：每個ri片段最大重復(fù)度，及其后繼片段在ri集合中的序號.

P1初始化,i=1.

P2 若i>n，則結(jié)束，否則轉(zhuǎn)P3.

P3 1) 取出ri的最后一個長度為d的子串，記為rd；

2) 令max_p=1及j=1；

3) 若j>n，則轉(zhuǎn)P4；否則轉(zhuǎn)4)；

4) 若j≠i，轉(zhuǎn)5，否則轉(zhuǎn)8)；

5) 令p=max_p；

6) 在ri中從p開始調(diào)用KMP算法找到所有子串rd出現(xiàn)的位置，記為B={b-1,b-2,…b-k}其中bi+1>bi,k表示在rd中rj出現(xiàn)的次數(shù)；

7) 令l=k,若l<1,則轉(zhuǎn)8)，否則執(zhí)行以下操作：

① 依次比較rj中第bi之前位置是否與ri中對應(yīng)字符相等，若直到rj中第一個字符均相等，則令max_p=bi+1，j=j+1，轉(zhuǎn)3).否則轉(zhuǎn)②

②l=l-1，轉(zhuǎn)①；

8)j=j+1，轉(zhuǎn)3).

P4 記錄ri的最大重復(fù)度為max_p+d-1及對應(yīng)后繼片段rj；i=i+1，轉(zhuǎn)P2.

當rd越大，兩個reads片段重復(fù)的概率越小，算法執(zhí)行越快，但會使許多reads無法連接到一起.另一方面，若rd越小，兩個reads片段連接的概率越大，算法執(zhí)行越慢，同時錯連概率增大.因此，先取rd=30建立Overlap圖，之后對無后繼節(jié)點取rd=10，完善Overlap圖. Matlab編程得到最后的44 042×2的矩陣readsl_used_youxiangtu，第ri行第1列表示第ri個片段的后繼片段是rj，第2列表示ri與rj的最大重復(fù)度.

3 Contigs的建立

基于 reads 之間的 Overlap 有向圖，拼接獲得的序列稱為 Contigs(重疊群)，Contigs的建立是連接成更長堿基鏈的過渡環(huán)節(jié).

3.1 模型建立

先列舉10個片段進行說明，假設(shè)這10個片段建立了如圖2的Overlap圖，第二列表示兩個片段的最大重復(fù)度.

Overlap圖說明：1的后繼片段是7，7的后繼片段是10，10沒有后繼片段.所以得出Contigs單鏈；2的后繼片段是3，4的后繼片段也是3，3的后繼片段是5，5沒有后繼片段，所以得出Contigs分叉鏈；6的后繼片段是8，8的后繼片段是9，9的后繼片段是8，所以得出Contigs環(huán)鏈,如圖3所示.

產(chǎn)生分叉和環(huán)的原因是在整條基因中存在很多堿基序列相同的區(qū)域，如圖4所示，陰影部分為堿基序列相同的兩個區(qū)域，虛線是切割位置，與4片段重復(fù)度最大的是3片段，于是錯把3片段當作4片段的后繼片段,因此產(chǎn)生了分叉.環(huán)的產(chǎn)生也類似.

實際上，基因中存在許多堿基序列一樣的區(qū)域，在這種情況下出現(xiàn)分叉和環(huán)的概率很大，而且還會出現(xiàn)分叉和環(huán)相組合的情況.

圖2 10個片段的Overlap圖Fig.2 Overlap figure of ten Congtigs

圖3 Contigs鏈Fig.3 Chain of Contigs

圖4 有重復(fù)片段基因圖Fig.4 Figure of duplication Contigs

3.2 數(shù)據(jù)結(jié)構(gòu)

基于以上分析，為了從Overlap有向圖找到Contigs，構(gòu)造如下的數(shù)據(jù)結(jié)構(gòu)Contigs.

1) 第i行表示第i條Contigs的信息;

2) 每一行有5個元素，第1列為該條Contigs的起點，第2列為終點，第3列為該條Contigs的類型,0表示單鏈，1表示分叉，2表示環(huán)；

3) 如果為分叉，第4列表示分叉接入的Contigs編號，第5列表示接入的位置;

4) 如果為環(huán)，第2列表示接入位置的前驅(qū)位置，第4列為0，第5列表示環(huán)接入的位置.

3.3 模型算法

Contigs連接算法描述：對建立的Overlap有向圖，尋找每一個有向路徑，即連接reads片段形成單支、環(huán)及分叉.

算法輸入：Overlap有向圖.

算法輸出：Contigs矩陣.

1) 將44 042×2的矩陣reads1_used_youxiangtu對應(yīng)建立一個44 042×2的矩陣biao_hao，將有后繼的節(jié)點標記為0，無后繼的節(jié)點標記為-1，即若biao_hao(i,1)=0表示第i個片段有后繼,biao_hao(i,1)=-1表示第i個片段無后繼.

2) 找到第1個標號為0的行，記為zerosl，并令Contigs_num=1.

3) 從zerosl之后開始，在矩陣biao_hao中找到第1個標號為0的行，若找到，則該行為第i行，做為第Contigs_num條Contigs的開頭片段；令zerosl=i，轉(zhuǎn)4)；否則結(jié)束.

4) 在對應(yīng)到Overlap圖中沿后繼往后找，同時將經(jīng)過的每個節(jié)點對應(yīng)的biao_hao改為Contigs_num.

如果節(jié)點biao_hao(i,1)=-1，Contigs_num=Contigs_num+1，則轉(zhuǎn)3)；如果節(jié)點標號biao_hao(i,1)=contigs_num，則轉(zhuǎn)5)；否則轉(zhuǎn)6).

5) (形成環(huán))記錄Contigs中Contigs_num行為環(huán)及相應(yīng)元素值；

6) (形成分叉)記錄Contigs中Contigs_num行為分叉及相應(yīng)元素值；

將形成的矩陣reads1_used_contigs按其類型分為矩陣danzhi_xulie和矩陣fencha_xulie，對形成的環(huán)Contigs從接入位置展開，合并到danzhi_xulie矩陣中. Matlab以此算法編程得到矩陣reads1_used_danzhi_xulie.

從矩陣reads1_used_danzhi_xulie中發(fā)現(xiàn)有的堿基片段無后繼片段，形成單個孤立的片段.單個孤立片段之所以屬于矩陣danzhi_xulie，是因為基因被切割后很多堿基序列是完全一樣的，這些片斷有后繼而且彼此之間的最大重復(fù)度是88，從接入位置展開后依然是片段本身.對于這些長度仍為88 bp的單鏈,Contigs將其從矩陣中刪除.Matlab編程不斷調(diào)試后發(fā)現(xiàn)rd=10時，即在ri片段中選取倒數(shù)10個堿基作為子串，建立Overlap圖和Contigs鏈效果最佳，reads1形成了3 574條Contigs單鏈，reads2形成3 415條Contigs單鏈.

Matlab依據(jù)Contigs算法編程得到矩陣reads1_used_fencha_xulie，如圖5所示.

圖5 Contigs矩陣Fig.5 Contigs matrix

對于形成的矩陣fencha_xulie，會形成比較復(fù)雜的情況，但都嫁接在單鏈上， 考慮選取一簇中最長的一條，Matlab編程實現(xiàn)找到每簇分叉中最長的單鏈，最后將reads1形成的3 574條和reads2形成的3 415條單鏈分別存入矩陣reads1_contigs_max_dazhi和reads2_contigs_max_dazhi中，再將其轉(zhuǎn)為字符串堿基存入Matlab軟件cell中，其中一條長度為282 bp的reads1_contigs堿基,如圖6所示：

圖6 單鏈Contigs Fig.6 Single stranded Contigs

圖7 堿基片段Fig.7 Basic group

如圖7, reads1、reads2是同一基因堿基片段互補鏈上的堿基，但reads2與reads1的測序方向相反，故將reads2形成的單鏈Contigs中的堿基先倒序然后再依據(jù)堿基互補原則(與堿基A結(jié)合的堿基必為T，與堿基C結(jié)合的堿基必為G)翻譯，并與reads1形成的單鏈Contigs合并，形成一個總共含有6 989個Contigs鏈的Contigs庫，對所有Contigs鏈進行一次Overlap圖和Contigs鏈建立，最終在10重復(fù)度下形成1 500條Contigs鏈.

4 ContigsN50的確定

4.1 結(jié)果分析

reads 拼接后會獲得一些不同長度的 Contigs.將所有的 Contigs長度相加，能獲得一個 Contigs 總長度.將所有的 Contigs 按照從長到短進行排序，再按照這個順序依次相加，當相加的長度達到 Contigs 總長度的一半時，最后一個加上的 Contigs 長度即為 ContigsN50.舉例：Contigs1+Contigs2+ Contigs3+Contigs4=Contigs 總長度的1/2 時，Contigs4 的長度即為 ContigsN50[3].ContigsN50 可以作為基因組拼接的結(jié)果好壞的一個判斷標準.

4.2 ContigsN50的確定

統(tǒng)計得到的1 500條Contigs鏈的長度，所有長度相加得到Contigs的總長為371 311 bp，將所有的 Contigs 按照從長到短進行排序，將這Contigs的長度從大到小依次累加，最終加到第462條Contigs時，達到總長的一半，第462條Contigs長度為284 bp，所以ContigsN50=284 bp，最長Contigs為1 226 bp .ContigsN50和最長Contigs的堿基序列如圖8、圖9所示.

圖8 ContigsN50的堿基序列Fig.8 Base sequence of ContigsN50

圖9 最長Contigs的堿基序列Fig.9 the longest base sequence of Contigs

5 模型計算時間

Overlap圖的運行時間：在子串最大重復(fù)度為30的情況下，對于reads1和reads2建立Overlap圖的運行時間分別在9 638.8 s、9 416 s；最大重復(fù)度為10的情況下， reads1和reads2分別在1 283.4 s、1 336 s建立了Overlap圖.

r1_r2_Contig的運行時間：重復(fù)度30運行時間為498.71 s;重復(fù)度10運行時間為485.701 s.

6 結(jié)語

對于海量的堿基片段，對數(shù)據(jù)依據(jù)相關(guān)文獻進行預(yù)處理，可以有效避免錯誤數(shù)據(jù)對后續(xù)拼接的干擾.處理中人工干預(yù)較少.在傳統(tǒng)的OLC算法上，運用KMP算法將原本的O(m×n)時間縮短為O(m+n)，在保證準確性和連續(xù)性不會降低的情況下提高了堿基段的拼接速度.進一步拼接形成3種Contigs鏈，對于單鏈和環(huán)給出了解決方法，但是出于拼接完整性的考慮，對分叉型Contigs處理還需進一步研究，另外重復(fù)度的選取也需進一步探討.

[1] 駱志剛,方小永，丁凡.DNA 序列拼接的研究進展及挑戰(zhàn)[J].計算機工程與科學(xué)，2007，29(8):127-127.

[2] 孫海汐，王秀杰.DNA測序技術(shù)發(fā)展及其展望[J].科研信息化技術(shù)與應(yīng)用，2009(3):18-18.

[3] 徐魁,陳科,徐君,等. CGDNA:基于簇圖的基因組序列集成拼接算法[J].計算機科學(xué), 2015, 42(9):235-239.

[4] 曾培龍.基于 reads 引導(dǎo)的基因組序列拼接算法[J]. 智能計算機與應(yīng)用，2015, 5(3)：23-25.

[5] 汪勇,張新,徐瓊,等.基因重組算法設(shè)計及多目標旅行商問題求解[J].系統(tǒng)工程, 2015(2):68-73.

[6] 耿麗,張仁杰.對基因組組裝算法的分析和研究[J].世界最新醫(yī)學(xué)信息文摘,2015(88):169-170.

[7] 毛華, 趙小娜, 史田敏,等. 多部圖的最大匹配算法[J]. 鄭州大學(xué)學(xué)報(理學(xué)版), 2013, 45(1):27-29.[8] 鄭緯民，林皎，羅水華.DNA 序列拼接中歐拉超路算法的新并行策略[J].計算機學(xué)報,2006，29(1): 139-139.

[9] 王旭.基于De Bruijn圖的DNA contig生成算法[D].哈爾濱:哈爾濱工業(yè)大學(xué),2011.

[10]嚴蔚敏.數(shù)據(jù)結(jié)構(gòu)(C語言版)[M].北京：清華大學(xué)出版社,2007.

[11]數(shù)學(xué)建模全國組委會.全國大學(xué)生數(shù)學(xué)建模競賽[EB/OL].[2014-04-19]. http：//www.Mcm.edu.cn.

[12]美吉生物網(wǎng).Illumina測序reads過濾[EB/OL]. [2014-05-19]. http：//www.majorbio.com.

(責(zé)任編輯：方惠敏)

The Improvement and Research about the OLC Algorithm in the Genetic Recombination

MAI Ali， YANG Wenwen

(DepartmentofAppliedMathematics,YunchengUniversity,Yuncheng044000,China)

It was great significance to obtain the genetic information of the organism quickly and accurately, then to obtain the sequence information of target creature genome for life science research. For genome assembly problem, four steps were proposed .They were the data preprocessing, the use of KMP algorithm in the time ofO(m+n) to quickly determine the maximum of certain segments of two bases the duplication, Contigs chain according to the Overlap figure, and the determination of ContigsN50. They could considerably improve existing OLC splicing technology. The optimized model and algorithm was presented, and the problem of the genome assembly could be solved in an easy way.

genetic recombination; KMP algorithm; Overlap figure; the chain of Contigs; ContigsN50

2015-11-17

國家自然科學(xué)基金資助項目(11526183);山西省基礎(chǔ)研究項目(2015021015);運城學(xué)院數(shù)學(xué)學(xué)科研究項目(XK-2014035, XK-2014030);運城學(xué)院博士啟動項目(YQ-2014011).

買阿麗(1981—),女,山西運城人,副教授,博士,主要從事微分方程理論及應(yīng)用研究,E-mail: maialiy@126.com.

買阿麗，楊雯雯．關(guān)于基因重組中OLC算法的改進研究[J]．鄭州大學(xué)學(xué)報(理學(xué)版)，2016，48(2)：34-39．

Q784

A

1671-6841(2016)02-0034-06

10.13705/j.issn.1671-6841.2015266