基于混合策略的de novo序列拼接算法構造

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

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

0 引言

以人類基因組計劃(human genome project,HGP)的實施為界,可以將生物信息學(Bioinformatics)的發展大致分為2個階段,分別為前基因組時代和后基因組時代.HGP的主要貢獻是提供了一項新的技術——基因測序[1].測序技術是進行生物研究最基本的技術之一,近幾年來,隨著測序技術的高速發展,這既降低了測序成本,又推動了生物信息學的發展,對精準醫療和基因檢測等研究產生了促進作用.

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

因為當前測序讀取的最大長度的序列也無法將整個基因組的遺傳信息表達準確,所以需要將經測序得到的短基因序列通過相關的技術拼接起來,得到一條完整的基因序列,這就是序列拼接(sequences assembly),其流程如圖1所示.

圖1 序列拼接流程圖

按照是否有參考序列參與的標準,序列拼接可分為2類:第1類是de novo(從頭)拼接,不需要參考序列,僅依靠短序列之間的重疊關系拼接成長序列;第2類是參考序列拼接,需要本物種或相似物種的基因組序列作為參考序列,以降低拼接過程的復雜度和成本.在實際應用中對一個新物種進行的測序往往沒有合適的參考序列,在這種情況下只能采取de novo方式進行序列拼接.

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

本文在分析同一拼接策略算法實現原理的基礎上,進一步研究在相同策略下的不同算法在具體實現過程中的差異性,得到算法在領域內的共同特征和差異特征.利用形式化方法及形式化平臺方面的優勢,結合領域分析建模和產生式編程的方法,形式化開發算法在領域內的序列拼接算法構件.對于共同特征可以生成單種策略下的必須執行構件,差異特征可以生成在單種策略下的可選執行構件[13].根據用戶對可選構件的選擇和輸入文件的特征自動組裝序列拼接算法.根據(OLC+DBG)→OLC混合模式將不同策略算法在領域內生成的拼接算法進行組裝,構造得到基于混合策略的ODO拼接算法.ODO算法在綜合不同策略優勢、提高拼接效率和準確度的同時也保證了結果的可靠性.

1 方法和技術

1.1 單一策略拼接

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

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

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

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

1.2 多策略混合拼接

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

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

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

圖3 ODO算法流程圖

2 算法構造

實現相同策略功能的算法可以是不同的,這是因為可按照使用需求生成同種單一策略下的不同算法.軟件形式化是實現軟件自動化的關鍵,PAR是一種統一的軟件形式化方法,包含了自定義泛型算法設計語言Radl、泛型抽象程序設計語言Apla、統一的算法程序設計方法以及PAR平臺(Radl到Apla程序生成系統,Apla到Java、C++等系列程序生成系統)和其他關鍵技術[28].根據序列拼接算法的需求設計出相應的程序規約,并在此基礎上進行程序設計,得到最終的算法和程序.

首先在領域層次上建立領域模型,然后利用產生式編程對算法構件進行交互設計,再形式化實現算法的構件庫,最后在不同策略的構件庫中自動生成不同單一策略下的拼接算法.利用生成的序列拼接算法按照(OLC+DBG)→OLC模式組裝后可得ODO算法,同時利用形式化方法保證了算法構造的高效性和可靠性.下面敘述了基于DBG策略拼接算法的構造過程,闡述了基于OLC策略拼接算法的構造過程.

2.1 DBG策略拼接算法生成

基于DBG策略開發的算法有SOAPdenovo、Velvet、ABySS、EULER等,將這些序列拼接算法聯合起來,構成了基于DBG策略的算法拼接領域.通過對在領域內算法的分析,找出算法的共同特征和可變特征并進行分類,選擇和定義需要分析和解決的領域,收集相關的領域信息,并形成領域模型[29].根據領域分析的結果,抽取算法共性,用參數表示算法的差異性,設計出Apla語言描述的拼接算法構件,為實現完整的算法構件庫,進一步分析不同構件之間的交互模型和構件之間的約束關系.

2.2 OLC策略拼接算法生成

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

圖4 OLC領域算法流程圖

圖5 OLC領域算法特征模型

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

圖6 算法構件交互模型

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

|[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表示標準核查序列B_seq和序列{ A,G,C,T}互為置換);

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

在上述的形式化規約中,輸入(in)和輸出(out)是在PAR平臺中定義的2個關鍵字標識符,序列類型(list)、整型(integer)和布爾型(boolean)是在PAR平臺中預定義的數據類型,AQ和AR分別是構件的前置條件和后置條件.

2.2.3 Apla形式化實現 抽象程序設計語言Apla可以直接使用抽象數據類型和抽象過程編寫程序,可抽象描述算法,并進行正確性驗證,這就保證了生成算法的可靠性,并且易于生成各種可執行程序設計語言程序Java、C++等.在Apla中的標識符、關鍵字、符號表達方式和在Radl語言中涉及的一致.

在使用形式化方法開發圖6中算法構件的Apla實現后,利用PAR平臺Apla→Java程序生成系統得到對應的Java代碼.

下面僅列出Seq_check構件的Apla實現.

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.

再得到相應的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策略拼接算法組裝 在構件設計完成后,根據使用需求和輸入文件的特征可選擇相對應的算法構件組裝ODO算法.本文選擇Seq_check、Error _correction(rocks and stones)、Sequence _alignment、Contigs_assembly等構件構造出了OLC_assembly_1算法.按照相同的方法,從OLC和DBG領域算法構件庫中選擇符合使用需求的算法構件,生成OLC_assembly_2和DBG_ assembly算法,進一步組裝得到ODO算法,由此將盡可能多的創造性勞動轉化為非創造性勞動.

3 材料和實驗

本次實驗選擇了3種單一策略(貪婪策略、OLC策略和DBG策略)的算法(SSAKE、Edena和Velvet)與ODO構造算法進行對比實驗.

3.1 實驗材料

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

表1 實驗樣本信息

表2 實驗運行平臺

3.2 coverage depth和k值對序列拼接的影響

coverage depth增加表示基因組復制的次數增加,經過鳥槍法(shot gun)剪切,可得到更多不同的reads,建立更完善的Hash表和更完整的有向圖.但若無限的增加coverage depth的值,則過多的reads會導致文庫過大,消耗大量的內存和計算時間,同時還會增加在有向圖中的分支和錯誤節點,降低拼接準確度.本次實驗將coverage depth的值設為50.

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

3.3 實驗結果與分析

構建了OLC和DBG領域構件庫,生成了基于OLC的拼接算法OLC_assembly_1、OLC_assembly_2和基于DBG的拼接算法DBG_ assembly,進一步組裝出(OLC+DBG)→OLC混合模式算法,簡稱ODO算法.為評估實驗效果,在表2的實驗平臺上用SSAKE(貪婪策略)、Edena(OLC策略)、Velvet(DBG策略)和構造生成的ODO算法分別對在表1中的3種基因樣本(Clostridiumperfringens、Acetobacteraceti、Escherichiacoli)進行序列拼接實驗.

由實驗結果(見表3)可知:構造的ODO算法無論是在N50大小上還是在基因組覆蓋度(coverage depth)上都優于其他的單一策略.如對于Clostridiumperfringens這一樣本基因,ODO算法拼接出來的結果N50為1 625 bp,基因覆蓋度為93.9%,優于單一策略進行序列拼接的結果.

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

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

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

4 結論

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