肖存威,石海鶴,王 嵐,程柏良

(江西師范大學(xué)計(jì)算機(jī)信息工程學(xué)院,江西 南昌 330022)

0 引言

以人類基因組計(jì)劃(human genome project,HGP)的實(shí)施為界,可以將生物信息學(xué)(Bioinformatics)的發(fā)展大致分為2個(gè)階段,分別為前基因組時(shí)代和后基因組時(shí)代.HGP的主要貢獻(xiàn)是提供了一項(xiàng)新的技術(shù)——基因測(cè)序[1].測(cè)序技術(shù)是進(jìn)行生物研究最基本的技術(shù)之一,近幾年來(lái),隨著測(cè)序技術(shù)的高速發(fā)展,這既降低了測(cè)序成本,又推動(dòng)了生物信息學(xué)的發(fā)展,對(duì)精準(zhǔn)醫(yī)療和基因檢測(cè)等研究產(chǎn)生了促進(jìn)作用.

在測(cè)序技術(shù)的發(fā)展過(guò)程中,一共產(chǎn)生了3代基因測(cè)序技術(shù).第1代基因測(cè)序技術(shù)在分子生物學(xué)研究中發(fā)揮了重要作用,比較知名的基因測(cè)序方法是雙脫氧法[2].第2代基因測(cè)序技術(shù)是在HGP的研究過(guò)程中被逐步開(kāi)發(fā)出來(lái),主要有Roche公司的454[3]、Illumina公司的Solexa[4]和ABI公司的Solid[5]等技術(shù)方法.第2代基因測(cè)序技術(shù)是通過(guò)對(duì)DNA片段復(fù)制擴(kuò)增,把堿基的信號(hào)強(qiáng)度放大,加入dNTP,在每個(gè)接頭上聚合新的堿基,就可以得到在整個(gè)接頭上的DNA序列信息[6].目前出現(xiàn)了以單分子測(cè)序?yàn)樘攸c(diǎn)的第3代基因測(cè)序技術(shù),使得在測(cè)序成本和測(cè)序通量上有了較大的改善.第3代基因測(cè)序技術(shù)主要以生物科學(xué)公司的HeliScope[7]、牛津納米孔技術(shù)公司的Nanopore[8]以及太平洋生物科學(xué)公司的PacBio[9]等技術(shù)方法為主流[10].與第2代基因測(cè)序技術(shù)不同的是,第3代基因測(cè)序技術(shù)既無(wú)需打斷整條DNA鏈,又無(wú)需對(duì)PCR擴(kuò)增,就可實(shí)現(xiàn)對(duì)DNA分子的單獨(dú)測(cè)序[11].測(cè)序讀取的最大長(zhǎng)度可達(dá)到12～15 kbp,從而減少了拼接次數(shù),降低了成本.不進(jìn)行PCR擴(kuò)增,避免了在擴(kuò)增過(guò)程中出現(xiàn)錯(cuò)誤,然而相比于測(cè)序的準(zhǔn)確率來(lái)說(shuō),第3代基因測(cè)序技術(shù)的平均錯(cuò)誤率約為15%,這遠(yuǎn)遠(yuǎn)高于第2代基因測(cè)序技術(shù)的平均錯(cuò)誤率(1%).

因?yàn)楫?dāng)前測(cè)序讀取的最大長(zhǎng)度的序列也無(wú)法將整個(gè)基因組的遺傳信息表達(dá)準(zhǔn)確,所以需要將經(jīng)測(cè)序得到的短基因序列通過(guò)相關(guān)的技術(shù)拼接起來(lái),得到一條完整的基因序列,這就是序列拼接(sequences assembly),其流程如圖1所示.

圖1 序列拼接流程圖

按照是否有參考序列參與的標(biāo)準(zhǔn),序列拼接可分為2類:第1類是de novo(從頭)拼接,不需要參考序列,僅依靠短序列之間的重疊關(guān)系拼接成長(zhǎng)序列;第2類是參考序列拼接,需要本物種或相似物種的基因組序列作為參考序列,以降低拼接過(guò)程的復(fù)雜度和成本.在實(shí)際應(yīng)用中對(duì)一個(gè)新物種進(jìn)行的測(cè)序往往沒(méi)有合適的參考序列,在這種情況下只能采取de novo方式進(jìn)行序列拼接.

隨著物種的范圍越來(lái)越大,序列拼接的準(zhǔn)確度和長(zhǎng)度要求越來(lái)越高,這給序列拼接工作帶來(lái)挑戰(zhàn).在拼接序列檢測(cè)方面,目前拼接得到的大多數(shù)序列是基于de novo拼接的,沒(méi)有參考基因組,無(wú)法從基因組層面對(duì)結(jié)果進(jìn)行比對(duì)檢測(cè),即使有參考基因組參與比對(duì)驗(yàn)證,在拼接過(guò)程中堿基也可能會(huì)發(fā)生突變,從而降低拼接的準(zhǔn)確度.在測(cè)序方面,目前的測(cè)序技術(shù)無(wú)法同時(shí)滿足低錯(cuò)誤率和讀長(zhǎng)較長(zhǎng)的要求,如第2代基因測(cè)序準(zhǔn)確度較高,但是讀長(zhǎng)較短、測(cè)序時(shí)間長(zhǎng)、費(fèi)用高,而第3代基因測(cè)序讀長(zhǎng)較長(zhǎng)、時(shí)間短、測(cè)序成本低,但準(zhǔn)確度較低.在序列拼接算法方面,大多算法是基于單一策略開(kāi)發(fā)的,具有一定的缺陷,如貪婪策略算法是利用啟發(fā)式搜索方式[12],選擇與種子reads重疊最多堿基的reads進(jìn)行拼接,最后得到局部最優(yōu)解或次優(yōu)解,卻無(wú)法保證得到全局最優(yōu)解.基于OLC(overlap-layout-consensus)策略和DBG(De Bruijngraph,De Bruijn圖)策略開(kāi)發(fā)的算法,在單獨(dú)使用時(shí)存在缺陷.對(duì)于復(fù)雜度較高的基因組,若單獨(dú)采用OLC策略,則雜合reads之間的同源性降低,容易導(dǎo)致重疊;若單獨(dú)采用DBG策略,則所形成的有向圖會(huì)存在較多的分叉,導(dǎo)致contigs長(zhǎng)度減小,無(wú)法進(jìn)一步拼接.因此,可以嘗試研究在混合策略下的序列拼接問(wèn)題,期望得到更好的序列拼接結(jié)果.

本文在分析同一拼接策略算法實(shí)現(xiàn)原理的基礎(chǔ)上,進(jìn)一步研究在相同策略下的不同算法在具體實(shí)現(xiàn)過(guò)程中的差異性,得到算法在領(lǐng)域內(nèi)的共同特征和差異特征.利用形式化方法及形式化平臺(tái)方面的優(yōu)勢(shì),結(jié)合領(lǐng)域分析建模和產(chǎn)生式編程的方法,形式化開(kāi)發(fā)算法在領(lǐng)域內(nèi)的序列拼接算法構(gòu)件.對(duì)于共同特征可以生成單種策略下的必須執(zhí)行構(gòu)件,差異特征可以生成在單種策略下的可選執(zhí)行構(gòu)件[13].根據(jù)用戶對(duì)可選構(gòu)件的選擇和輸入文件的特征自動(dòng)組裝序列拼接算法.根據(jù)(OLC+DBG)→OLC混合模式將不同策略算法在領(lǐng)域內(nèi)生成的拼接算法進(jìn)行組裝,構(gòu)造得到基于混合策略的ODO拼接算法.ODO算法在綜合不同策略優(yōu)勢(shì)、提高拼接效率和準(zhǔn)確度的同時(shí)也保證了結(jié)果的可靠性.

1 方法和技術(shù)

1.1 單一策略拼接

基于貪婪策略開(kāi)發(fā)的序列拼接算法是在整個(gè)序列拼接研究過(guò)程中最早提出的一類算法,也被稱為貪心算法,是從問(wèn)題的某一個(gè)初始值出發(fā),逐步從兩邊逼近目標(biāo)值,以盡快地求得一個(gè)較好的解,當(dāng)算法運(yùn)行到某一步不能繼續(xù)前進(jìn)時(shí),算法停止.運(yùn)行結(jié)束后會(huì)得到一個(gè)局部最優(yōu)解或次優(yōu)解,在運(yùn)行過(guò)程中這一類算法所做的每個(gè)選擇都是在當(dāng)前狀態(tài)下最好的選擇(貪心選擇),貪心策略拼接原理示意圖如圖2(a)所示.常見(jiàn)的基于貪婪策略開(kāi)發(fā)的序列拼接算法有SSAKE[14]、VCAKE[15]、SHARCGS[16]、NPSCARF[17]等.

OLC策略主要基于reads之間的重疊關(guān)系進(jìn)行拼接,是最能直觀體現(xiàn)序列拼接原理思想的策略.OLC策略在序列拼接時(shí)將所有待拼接的reads構(gòu)造成一個(gè)有向圖,圖中的每個(gè)節(jié)點(diǎn)都代表特定的reads;若2個(gè)節(jié)點(diǎn)之間存在邊,則說(shuō)明這2個(gè)reads的重疊部分的堿基數(shù)量大于設(shè)置的閾值.對(duì)圖進(jìn)行簡(jiǎn)化和分割后,可以在一系列子圖上尋找經(jīng)過(guò)每個(gè)節(jié)點(diǎn)一次且僅一次的路徑,最后得到所需的序列.由此把序列拼接問(wèn)題轉(zhuǎn)化為Hamilton路徑問(wèn)題,OLC策略拼接原理示意圖如圖2(b)所示.基于OLC策略開(kāi)發(fā)的序列拼接算法有CABOG[18]、Phrap[19]、Celera Assembly[20]、Edena[21]等.

DBG策略的特點(diǎn)是不需要進(jìn)行序列比對(duì),可大幅度減少系統(tǒng)內(nèi)存的消耗,De Bruijn圖的大小只與基因組大小和算法復(fù)雜度有關(guān).在DBG算法中,把reads按照長(zhǎng)度k分割成k-mer,并將這些k-mer存入Hash表中,不同k-mer在Hash表中只存儲(chǔ)1次,然后將這些k-mer作為點(diǎn),構(gòu)建De Bruijn圖.對(duì)De Bruijn圖進(jìn)行簡(jiǎn)化和分割,并在圖中找到1條歐拉路徑構(gòu)成contig,再通過(guò)Pair-end文庫(kù)可以確定contig之間的位置關(guān)系,指導(dǎo)contig形成長(zhǎng)度更長(zhǎng)的scaffold,DBG策略拼接原理示意圖如圖2(c)所示.基于DBG策略開(kāi)發(fā)的序列拼接算法有SOAPdenovo[22]、Velvet[23]、ABySS[24]、EULER[25]、ALLPATHS[26]等.

圖2 單一策略拼接原理示意圖

1.2 多策略混合拼接

貪婪策略在拼接序列時(shí)總是做出在當(dāng)前情況下最好的選擇,不從整體加以考慮,且并非對(duì)所有問(wèn)題都適用.OLC策略在使用時(shí)因PCR擴(kuò)增而會(huì)存在大量的重復(fù)數(shù)據(jù),并且對(duì)長(zhǎng)reads拼接有偏向性.DBG策略在使用時(shí)將短reads重新分割成長(zhǎng)度更短的k-mer進(jìn)行序列拼接,增加了計(jì)算時(shí)間和拼接的復(fù)雜度.由此可見(jiàn)這些策略在單獨(dú)使用時(shí)均存在不足.

為了解決單一策略在序列拼接時(shí)易出現(xiàn)的問(wèn)題,綜合各個(gè)策略的優(yōu)點(diǎn),這可以使用混合模式來(lái)實(shí)現(xiàn)策略的混合.如(貪婪策略+DBG)→OLC、貪婪策略+OLC、貪婪策略+DBG、DBG+OLC、(OLC+DBG)→DBG、(OLC+DBG)→OLC、(OLC+OLC)→DBG以及(DBG+DBG)→OLC等模式.A.K. Kusuma等[27]對(duì)(OLC+DBG)→OLC進(jìn)行了研究,在實(shí)驗(yàn)過(guò)程中,首先使用Velvet算法拼接reads,然后使用Edena算法對(duì)相同的reads再次拼接,最后使用Minimus assembler對(duì)拼接產(chǎn)生的contigs再次拼接可得到最終結(jié)果.使用上述方法,A.K. Kusuma等[27]在3個(gè)真實(shí)數(shù)據(jù)集和4個(gè)模擬數(shù)據(jù)集上進(jìn)行了實(shí)驗(yàn).該方法既避免了DBG策略拼接導(dǎo)致的錯(cuò)誤率過(guò)高的現(xiàn)象發(fā)生,又解決了OLC策略拼接帶來(lái)的重復(fù)序列和占用內(nèi)存過(guò)大的問(wèn)題.

本文充分利用在軟件開(kāi)發(fā)方法上及其在平臺(tái)方面具有的優(yōu)勢(shì),對(duì)(OLC+DBG)→OLC混合模式的算法進(jìn)行了研究.使用形式化、產(chǎn)生式編程等方法和技術(shù),分別構(gòu)造了基于OLC策略和基于DBG策略的算法域構(gòu)件庫(kù),借助在PAR平臺(tái)上的Apla→Java程序生成系統(tǒng)生成了基于OLC的拼接算法OLC_assembly_1、OLC_assembly_2和基于DBG的拼接算法DBG_assembly,進(jìn)一步拼裝出混合模式的ODO算法,從而將盡可能多的創(chuàng)造性勞動(dòng)轉(zhuǎn)化為非創(chuàng)造性勞動(dòng),提高了算法開(kāi)發(fā)過(guò)程的可理解性、算法的可靠性以及算法構(gòu)件庫(kù)在相關(guān)領(lǐng)域中的可重用性,ODO算法流程圖如圖3所示.

圖3 ODO算法流程圖

2 算法構(gòu)造

實(shí)現(xiàn)相同策略功能的算法可以是不同的,這是因?yàn)榭砂凑帐褂眯枨笊赏N單一策略下的不同算法.軟件形式化是實(shí)現(xiàn)軟件自動(dòng)化的關(guān)鍵,PAR是一種統(tǒng)一的軟件形式化方法,包含了自定義泛型算法設(shè)計(jì)語(yǔ)言Radl、泛型抽象程序設(shè)計(jì)語(yǔ)言Apla、統(tǒng)一的算法程序設(shè)計(jì)方法以及PAR平臺(tái)(Radl到Apla程序生成系統(tǒng),Apla到Java、C++等系列程序生成系統(tǒng))和其他關(guān)鍵技術(shù)[28].根據(jù)序列拼接算法的需求設(shè)計(jì)出相應(yīng)的程序規(guī)約,并在此基礎(chǔ)上進(jìn)行程序設(shè)計(jì),得到最終的算法和程序.

首先在領(lǐng)域?qū)哟紊辖㈩I(lǐng)域模型,然后利用產(chǎn)生式編程對(duì)算法構(gòu)件進(jìn)行交互設(shè)計(jì),再形式化實(shí)現(xiàn)算法的構(gòu)件庫(kù),最后在不同策略的構(gòu)件庫(kù)中自動(dòng)生成不同單一策略下的拼接算法.利用生成的序列拼接算法按照(OLC+DBG)→OLC模式組裝后可得ODO算法,同時(shí)利用形式化方法保證了算法構(gòu)造的高效性和可靠性.下面敘述了基于DBG策略拼接算法的構(gòu)造過(guò)程,闡述了基于OLC策略拼接算法的構(gòu)造過(guò)程.

2.1 DBG策略拼接算法生成

基于DBG策略開(kāi)發(fā)的算法有SOAPdenovo、Velvet、ABySS、EULER等,將這些序列拼接算法聯(lián)合起來(lái),構(gòu)成了基于DBG策略的算法拼接領(lǐng)域.通過(guò)對(duì)在領(lǐng)域內(nèi)算法的分析,找出算法的共同特征和可變特征并進(jìn)行分類,選擇和定義需要分析和解決的領(lǐng)域,收集相關(guān)的領(lǐng)域信息,并形成領(lǐng)域模型[29].根據(jù)領(lǐng)域分析的結(jié)果,抽取算法共性,用參數(shù)表示算法的差異性,設(shè)計(jì)出Apla語(yǔ)言描述的拼接算法構(gòu)件,為實(shí)現(xiàn)完整的算法構(gòu)件庫(kù),進(jìn)一步分析不同構(gòu)件之間的交互模型和構(gòu)件之間的約束關(guān)系.

2.2 OLC策略拼接算法生成

2.2.1 領(lǐng)域特征模型分析 將基于OLC策略的序列拼接算法聯(lián)合起來(lái),構(gòu)成OLC算法領(lǐng)域,在領(lǐng)域內(nèi)的序列拼接算法在進(jìn)行序列拼接時(shí)的過(guò)程大致可以分為3步:首先,通過(guò)序列比對(duì)找到reads之間重疊長(zhǎng)度超過(guò)設(shè)置閾值的重疊關(guān)系;其次,通過(guò)這些重疊信息將reads建立一種新的組合關(guān)系,可得到重疊群contigs,再將contigs進(jìn)一步排列,生成較長(zhǎng)的scaffold重疊群;最后,在經(jīng)過(guò)簡(jiǎn)化和分割后的有向圖中找到一條最優(yōu)的路徑(Hamilton路徑)所對(duì)應(yīng)的序列,這即是需要的拼接序列.接下來(lái)對(duì)在OLC算法領(lǐng)域中的算法CABOG、Celera Assembly和Edena進(jìn)行分析,找出這3個(gè)算法的共同特性和可變特性.如在實(shí)現(xiàn)序列核查、序列比對(duì)、構(gòu)建有向圖等功能中,這3種算法使用的方法和步驟是類似的.但在糾正測(cè)序錯(cuò)誤這個(gè)步驟中,Edena使用快速檢測(cè)虛假reads來(lái)實(shí)現(xiàn)對(duì)測(cè)序錯(cuò)誤的糾正,CABOG使用rocks和stones技術(shù)來(lái)實(shí)現(xiàn)對(duì)測(cè)序錯(cuò)誤的糾正,Celera Assembly使用PacBioToCA自糾算法來(lái)實(shí)現(xiàn)對(duì)測(cè)序錯(cuò)誤的糾正.在對(duì)這3個(gè)算法進(jìn)行領(lǐng)域分析后,找到了3個(gè)算法的共同特性和不同特性,歸納得到了OLC領(lǐng)域算法流程圖和領(lǐng)域特征模型(見(jiàn)圖4和圖5).

圖4 OLC領(lǐng)域算法流程圖

圖5 OLC領(lǐng)域算法特征模型

2.2.2 領(lǐng)域構(gòu)件設(shè)計(jì) 根據(jù)對(duì)OLC領(lǐng)域的分析和建模,可抽取出其中的特征,用參數(shù)表示其中的差異性,設(shè)計(jì)出可使用Apla語(yǔ)言描述的算法構(gòu)件[29].在構(gòu)建得到ODO算法后,測(cè)序所得的序列作為輸入,算法將會(huì)對(duì)輸入的序列做核查,避免由硬件問(wèn)題而出現(xiàn)字符錯(cuò)誤的現(xiàn)象發(fā)生,在糾正測(cè)序錯(cuò)誤后,再進(jìn)行序列比對(duì),構(gòu)建有向圖,去除錯(cuò)誤分支,生成contigs和scaffolds.因此,這種OLC算法的主要構(gòu)件是序列核查構(gòu)件(Seq_check)、糾錯(cuò)構(gòu)件(Error_correction)、序列比對(duì)構(gòu)件(Sequence_alignment)、有向圖生成構(gòu)件(Graph_structure)、圖化簡(jiǎn)去分支構(gòu)件(Graph_simplification)、Contigs拼接構(gòu)件(Contigs_assembly)、Scaffolds拼接構(gòu)件(Scaffolds_assembly)和最長(zhǎng)路徑構(gòu)件(Longest_path).算法構(gòu)件關(guān)系圖如圖6所示,這里使用Radl語(yǔ)言為Seq_check構(gòu)件的功能做形式化描述.

圖6 算法構(gòu)件交互模型

Seq_check構(gòu)件:此構(gòu)件對(duì)輸入的基因序列進(jìn)行核查,如在DNA基因序列中僅含有4種堿基(A,G,C,T),若序列中含有其他堿基,則會(huì)停止運(yùn)行并提示錯(cuò)誤.

|[in A_seq:list of char;B_seq:list of char;n:integer;out_flag:boolean]|;

AQ:n=2∧|A_seq|>0∧|B_seq|=4∧perm(

B_seq,{A,G,C,T});

(perm表示標(biāo)準(zhǔn)核查序列B_seq和序列{ A,G,C,T}互為置換);

AR:flag=(?i:0≤i≤|A_seq|:(?j:0≤j≤3:A_seq [i]=B_seq [j])).

在上述的形式化規(guī)約中,輸入(in)和輸出(out)是在PAR平臺(tái)中定義的2個(gè)關(guān)鍵字標(biāo)識(shí)符,序列類型(list)、整型(integer)和布爾型(boolean)是在PAR平臺(tái)中預(yù)定義的數(shù)據(jù)類型,AQ和AR分別是構(gòu)件的前置條件和后置條件.

2.2.3 Apla形式化實(shí)現(xiàn) 抽象程序設(shè)計(jì)語(yǔ)言Apla可以直接使用抽象數(shù)據(jù)類型和抽象過(guò)程編寫(xiě)程序,可抽象描述算法,并進(jìn)行正確性驗(yàn)證,這就保證了生成算法的可靠性,并且易于生成各種可執(zhí)行程序設(shè)計(jì)語(yǔ)言程序Java、C++等.在Apla中的標(biāo)識(shí)符、關(guān)鍵字、符號(hào)表達(dá)方式和在Radl語(yǔ)言中涉及的一致.

在使用形式化方法開(kāi)發(fā)圖6中算法構(gòu)件的Apla實(shí)現(xiàn)后,利用PAR平臺(tái)Apla→Java程序生成系統(tǒng)得到對(duì)應(yīng)的Java代碼.

下面僅列出Seq_check構(gòu)件的Apla實(shí)現(xiàn).

function Seq_check(A_seq:list(char);B_seq:list(char,4)):boolean;

vari,j,sum:integer;

begin

i,j,sum:=0,0,0;

do (0

do (j<3)→if(A_seq[i]==B_seq[j])→

sum:=sum+1;

fi;

j:=j+1;

od;

i:=i+1;

od;

if(sum==|A_seq|) →Seq_check:=true

[]→Seq_check:=false;

fi;

end.

再得到相應(yīng)的Java程序如下:

public boolean Seq_check(String A_seq,String B_seq){

inti=0;

intj=0;

int sum=0;

for (i=0;i

for(j=0;j<3;j++){

if(A_seq[i]==B_seq[j]){

sum++;

}

}

}

if(sum == A_seq.length ()){

return true;

}else{

return false;

}

}.

2.2.4 OLC策略拼接算法組裝 在構(gòu)件設(shè)計(jì)完成后,根據(jù)使用需求和輸入文件的特征可選擇相對(duì)應(yīng)的算法構(gòu)件組裝ODO算法.本文選擇Seq_check、Error _correction(rocks and stones)、Sequence _alignment、Contigs_assembly等構(gòu)件構(gòu)造出了OLC_assembly_1算法.按照相同的方法,從OLC和DBG領(lǐng)域算法構(gòu)件庫(kù)中選擇符合使用需求的算法構(gòu)件,生成OLC_assembly_2和DBG_ assembly算法,進(jìn)一步組裝得到ODO算法,由此將盡可能多的創(chuàng)造性勞動(dòng)轉(zhuǎn)化為非創(chuàng)造性勞動(dòng).

3 材料和實(shí)驗(yàn)

本次實(shí)驗(yàn)選擇了3種單一策略(貪婪策略、OLC策略和DBG策略)的算法(SSAKE、Edena和Velvet)與ODO構(gòu)造算法進(jìn)行對(duì)比實(shí)驗(yàn).

3.1 實(shí)驗(yàn)材料

本文從GenBank中下載了3種基因組較小的樣本來(lái)進(jìn)行實(shí)驗(yàn),實(shí)驗(yàn)樣本信息如表1所示,基因名稱分別是Clostridiumperfringens(產(chǎn)氣莢膜梭菌)、Acetobacteraceti(醋酸桿菌)、Escherichiacoli(大腸桿菌),該實(shí)驗(yàn)在如表2所示的實(shí)驗(yàn)平臺(tái)上運(yùn)行.

表1 實(shí)驗(yàn)樣本信息

表2 實(shí)驗(yàn)運(yùn)行平臺(tái)

3.2 coverage depth和k值對(duì)序列拼接的影響

coverage depth增加表示基因組復(fù)制的次數(shù)增加,經(jīng)過(guò)鳥(niǎo)槍法(shot gun)剪切,可得到更多不同的reads,建立更完善的Hash表和更完整的有向圖.但若無(wú)限的增加coverage depth的值,則過(guò)多的reads會(huì)導(dǎo)致文庫(kù)過(guò)大,消耗大量的內(nèi)存和計(jì)算時(shí)間,同時(shí)還會(huì)增加在有向圖中的分支和錯(cuò)誤節(jié)點(diǎn),降低拼接準(zhǔn)確度.本次實(shí)驗(yàn)將coverage depth的值設(shè)為50.

Velvet算法和ODO算法在構(gòu)建De Bruijn圖時(shí)會(huì)將文庫(kù)中的reads切割成長(zhǎng)度為k的k-mer.k值越大越容易得到唯一的序列,但在拼接過(guò)程中越容易產(chǎn)生gaps;隨著k值減小,可以增加De Bruijn圖的連通性,但會(huì)提高算法的復(fù)雜度,增加處理時(shí)間.同時(shí),為了避免正反鏈混淆產(chǎn)生回文序列的現(xiàn)象發(fā)生,k值應(yīng)選奇數(shù).本次實(shí)驗(yàn)所使用的k值取為31.

3.3 實(shí)驗(yàn)結(jié)果與分析

構(gòu)建了OLC和DBG領(lǐng)域構(gòu)件庫(kù),生成了基于OLC的拼接算法OLC_assembly_1、OLC_assembly_2和基于DBG的拼接算法DBG_ assembly,進(jìn)一步組裝出(OLC+DBG)→OLC混合模式算法,簡(jiǎn)稱ODO算法.為評(píng)估實(shí)驗(yàn)效果,在表2的實(shí)驗(yàn)平臺(tái)上用SSAKE(貪婪策略)、Edena(OLC策略)、Velvet(DBG策略)和構(gòu)造生成的ODO算法分別對(duì)在表1中的3種基因樣本(Clostridiumperfringens、Acetobacteraceti、Escherichiacoli)進(jìn)行序列拼接實(shí)驗(yàn).

由實(shí)驗(yàn)結(jié)果(見(jiàn)表3)可知:構(gòu)造的ODO算法無(wú)論是在N50大小上還是在基因組覆蓋度(coverage depth)上都優(yōu)于其他的單一策略.如對(duì)于Clostridiumperfringens這一樣本基因,ODO算法拼接出來(lái)的結(jié)果N50為1 625 bp,基因覆蓋度為93.9%,優(yōu)于單一策略進(jìn)行序列拼接的結(jié)果.

表3 4種策略算法在不同樣本上的拼接結(jié)果

在評(píng)估指標(biāo)中N50表示在拼接得到的Contigs長(zhǎng)度從大到小排列后最中間Contigs長(zhǎng)度的值,N50的值越大說(shuō)明拼接得到的Contigs長(zhǎng)度越長(zhǎng),即拼接效果越好.在對(duì)Escherichiacoli實(shí)驗(yàn)樣本進(jìn)行拼接時(shí),SSAKE、Edena、Velvet和ODO得到的N50值分別為328 bp、545 bp、673 bp和894 bp.由實(shí)驗(yàn)結(jié)果可知:相比于其他的單一策略下的算法,ODO算法取得了較優(yōu)的結(jié)果.產(chǎn)生這一結(jié)果的主要原因有2個(gè)方面:(i)先對(duì)輸入的reads文庫(kù)進(jìn)行2次拼接,然后將2次拼接的結(jié)果混合再次作為第3階段拼接算法的輸入,得到最終的Contigs輸出;(ii)在ODO算法的構(gòu)造過(guò)程中使用形式化方法開(kāi)發(fā)算法領(lǐng)域構(gòu)件,可根據(jù)用戶的需求和輸入文件的特征選擇不用的構(gòu)件組裝得到策略拼接算法,在提高算法準(zhǔn)確度和效率的同時(shí),保證了結(jié)果的可靠性.

在評(píng)估指標(biāo)中Contigs number表示拼接結(jié)果的Contigs數(shù)量,在基因組堿基數(shù)量大致相同的情況下,Contigs數(shù)量越少表示每條Contigs的長(zhǎng)度越長(zhǎng),即拼接結(jié)果越好.在對(duì)Acetobacteraceti實(shí)驗(yàn)樣本進(jìn)行拼接時(shí),SSAKE、Edena、Velvet和ODO得到的Contigs number值分別為106、87、54和67.由評(píng)估指標(biāo)可知:ODO算法的結(jié)果比SSAKE和Edena算法的結(jié)果更優(yōu),但是離Velvet算法的結(jié)果還有一定的差距.產(chǎn)生這一結(jié)果的可能原因是:在序列中還有較多的gaps,當(dāng)堿基比對(duì)遇到gap時(shí)算法會(huì)默認(rèn)當(dāng)前匹配不成功,并停止匹配,輸出結(jié)果,這將導(dǎo)致Contigs number反向增大.

4 結(jié)論

本文利用形式化方法及形式化平臺(tái)方面的優(yōu)勢(shì),結(jié)合領(lǐng)域分析建模和產(chǎn)生式編程的方法,形式化開(kāi)發(fā)了序列拼接算法構(gòu)件,并構(gòu)造了(OLC+DBG)→OLC混合模式算法,簡(jiǎn)稱ODO算法.實(shí)驗(yàn)結(jié)果顯示:在對(duì)3種基因樣本進(jìn)行序列拼接對(duì)比實(shí)驗(yàn)時(shí),構(gòu)造生成的ODO算法比單一策略算法所取得的結(jié)果在N50和Coverage等參數(shù)上有一定的優(yōu)勢(shì).ODO算法在Escherichiacoli樣本的序列拼接實(shí)驗(yàn)中基因組覆蓋度達(dá)到97.1%,N50的值達(dá)到894 bp.研究結(jié)果表明:利用產(chǎn)生式編程和形式化方法開(kāi)發(fā)的ODO算法在綜合不同策略優(yōu)勢(shì)、提高拼接效率和準(zhǔn)確度的同時(shí)也保證了結(jié)果的可靠性,為深入研究算法構(gòu)造在序列拼接上的影響提供了參考.