王 淋,鐘 誠

(廣西大學(xué) 計(jì)算機(jī)與電子信息學(xué)院,南寧 530004) (廣西高校并行分布與智能計(jì)算重點(diǎn)實(shí)驗(yàn)室,南寧 530004)

1 引 言

序列比對在序列同源性分析、生物序列保守位點(diǎn)識別和系統(tǒng)發(fā)育分析等領(lǐng)域具有重要作用[1,2].多序列比對(Multiple Sequence Alignment,MSA)問題是指把兩條以上的生物序列進(jìn)行比對,使相同殘基盡可能排在同一列上,以找出這些序列的同源區(qū)域[3,4].

從算法設(shè)計(jì)策略角度,多序列比對方法可以分為樹比對、星比對、迭代比對、基于強(qiáng)化學(xué)習(xí)的比對等.樹比對方法先對序列進(jìn)行層次聚類形成引導(dǎo)樹,再從引導(dǎo)樹的葉節(jié)點(diǎn)自底向上逐步對準(zhǔn)序列形成最終的比對結(jié)果,代表算法有ClustalW[5]、Kalign[6]、ClustalO[7]和FAMSA[8]等.星比對方法先選取一條中心序列,然后將其他序列分別和中心序列執(zhí)行對準(zhǔn)以形成比對結(jié)果,相關(guān)算法有MASC[9]和HAlign[10]等.樹比對方法和星比對方法都無法糾正序列比對前期的錯誤,其比對結(jié)果質(zhì)量依賴于比對順序和中心序列的選取[11].迭代式比對方法通過多次迭代糾正序列比對前期引入的錯誤,以提升比對的準(zhǔn)確率,相關(guān)代表算法有MUSCLE[12]、MAFFT[13]、MO-SAStrE[14]、SAPSO[15]和Sequoya[16]等.為進(jìn)一步提升比對的準(zhǔn)確率,Jafari等人研究了將深度強(qiáng)化學(xué)習(xí)應(yīng)用于求解多序列比對問題[17].文獻(xiàn)[18]提出了一種基于強(qiáng)化學(xué)習(xí)的序列比對算法RLALIGN.迭代式比對算法運(yùn)行速度較慢,基于強(qiáng)化學(xué)習(xí)的比對算法訓(xùn)練時間長且結(jié)果難收斂,兩者都難以適用于大規(guī)模長序列數(shù)據(jù)集的比對.

第三代測序技術(shù)的發(fā)展產(chǎn)生了越來越長的序列,這給大規(guī)模多序列比對帶來了挑戰(zhàn)[19].為減少多序列比對算法的運(yùn)行時間,一種方法是將多序列比對串行算法并行化.針對多序列比對算法常用的并行化技術(shù)包括多CPU并行[20,21]、CPU和GPU協(xié)同并行[22,23]等.另一種方法是基于劃分的思想,將輸入序列集拆分比對.算法VDGA將輸入序列集中每條長序列垂直分割成多個子序列,分別比對每個子序列,再組合比對結(jié)果[24].算法FAME利用求解公共子串的方法直接確定長序列的拆分位置[25],省去了算法VDGA拆分前的比對操作.

為進(jìn)一步有效比對大規(guī)模長序列,本文基于劃分思想,提出一種融合水平聚類分簇和垂直分組的多序列比對算法,以在維持較高比對精度的同時,顯著加快比對的完成.

2 方 法

設(shè)N條長序列的集合S={S0,S1,…,SN-1},序列Si的長度為Li,i=0,1,2,…,N-1.本文提出的融合聚類分簇和垂直分組的大規(guī)模長序列比對的思想是:第1階段,通過融合采用mBed方法和簡并字母表方法對S中每條長序列編碼,將長序列集S轉(zhuǎn)化為數(shù)值短向量集D,利用二分k-means算法聚類D,將S按序列相似度劃分為H個水平簇C0,C1,…,CH-1.第2階段,為每個簇Cj內(nèi)的長序列集構(gòu)建最長兼容鏈,并將簇Cj垂直分割為Vj個包含較短序列的垂直分組Cj,0,Cj,1,…,Cj,Vj-1,j=0,1,2,…,H-1.第3階段,對簇Cj中的每個垂直分組Cj,p,若其是非最長兼容鏈選中的垂直分組,則使用已有的適用于規(guī)模較小的多序列比對算法對準(zhǔn)Cj,p中的序列,以得到簇Cj的完整比對結(jié)果,j=0,1,2,…,H-1,p=0,1,2,…,Vj-1.第4階段,構(gòu)建水平簇集{C0,C1,…,CH-1}的最長兼容鏈,將序列集S劃分為NG個簇間分組CG0,CG1,CG2,…,CGNG-1,利用本文提出的帶有Gap類型推斷的動態(tài)規(guī)劃漸進(jìn)比對方法分別對準(zhǔn)每個非最長兼容鏈選中的簇間分組,以獲得最終的N條長序列的比對結(jié)果.

下面詳細(xì)闡述融合水平聚類分簇和垂直分組策略對N條長序列進(jìn)行比對的方法.

2.1 多序列聚類分簇

聚類分簇是指對長序列集S進(jìn)行聚類,按序列相似度將S劃分成H個水平簇C0,C1,…,CH-1,使得簇Cj中各條序列盡可能相似,j=0,1,2,…,H-1.聚類分簇使得簇內(nèi)序列擁有更多公共子片段且大部分公共子片段在垂直分組過程中直接對準(zhǔn),這將極大地降低第3階段需要比對的數(shù)據(jù)規(guī)模.聚類分簇工作首先對S中每條長序列編碼,將長序列集S轉(zhuǎn)化為數(shù)值短向量集D,然后通過對D進(jìn)行聚類,將長序列集S分簇,同時生成引導(dǎo)樹GT用于簇合并.

為降低數(shù)值向量的維數(shù),以加速大規(guī)模長序列的聚類,首先采用mBed方法[26]從序列集S中選取m條序列記為M0,M1,…,Mm-1;其次,采用簡并字母表方法[27]分別計(jì)算序列S0,S1,…,SN-1的k-mer編碼向量;然后,計(jì)算Si與Mj對應(yīng)的k-mer編碼向量的歐氏空間距離D(Si,Mj),j=0,1,2,…,m-1,i=0,1,2,…,N-1,m=(log2N)2.最后,將(D(Si,M0),D(Si,M1),…,D(Si,Mm-1))作為Si對應(yīng)的數(shù)值向量Di,i=0,1,2,…N-1.將S中每條序列編碼得到數(shù)值短向量集D.通過二分k-means算法[28]對D進(jìn)行聚類,以實(shí)現(xiàn)序列集S的水平劃分.

算法1.MS-kCluster

輸入：長序列數(shù)據(jù)集S={S0,S1,…,SN-1}

輸出：劃分S后的H個水平簇C0,C1,…,CH-1的集合Clusters、引導(dǎo)樹GT

Begin:

1.采用mBed方法從長序列集S中選取m=(log2N)2條長序列,記為M0,M1,…,Mm-1;

2.fori=0toN-1do

3. 通過簡并字母表方法計(jì)算長序列Si的k-mer編碼向量KV(Si);

4.endfor

5.fori=0toN-1do

6.forj=0tom-1do

7. 計(jì)算KV(Si)與KV(Mj)的歐氏空間距離D(Si,Mj);

8.endfor

9. 計(jì)算Si對應(yīng)的數(shù)值向量Di=(D(Si,M0),D(Si,M1),…,D(Si,Mm-1));

10.endfor

11.node←Node(D);//Node(D)表示初始化引導(dǎo)樹GT的新節(jié)點(diǎn),并將數(shù)值向量集合D作為新節(jié)點(diǎn)的值

12.初始化引導(dǎo)樹GT←node且樹當(dāng)前葉節(jié)點(diǎn)的集合CS←{node};

13.while(CS中節(jié)點(diǎn)數(shù)量

14.node←CS中具有最大SSE值的節(jié)點(diǎn);

15. 使用k-means聚類算法將節(jié)點(diǎn)node中包含的數(shù)值向量集合分割為兩個子集合Dleft、Dright;

16.node.left←Node(Dleft);//node.left表示引導(dǎo)樹節(jié)點(diǎn)node的左子節(jié)點(diǎn)

17.node.right←Node(Dright);//node.right表示引導(dǎo)樹節(jié)點(diǎn)node的右子節(jié)點(diǎn)

18. 刪除CS集中的node節(jié)點(diǎn),并將node.left和node.right加入到CS集;

19.endwhile

20.Clusters←{};

21.fornodeinCSdo

22.Clusters←Clusters∪{節(jié)點(diǎn)node中數(shù)值向量集對應(yīng)的序列集合};

23.endfor

End

2.2 簇內(nèi)序列垂直分組

經(jīng)聚類分簇后,長序列集S被劃分為H個簇C0,C1,…,CH-1.每個簇內(nèi)序列數(shù)量大幅減少,但簇內(nèi)序列的長度沒有改變,若直接采用現(xiàn)有的多序列比對算法對簇內(nèi)的長序列進(jìn)行比對,算法仍需要較高的計(jì)算開銷.為降低大規(guī)模長序列比對的時間,需在序列長度維度上對簇內(nèi)序列進(jìn)一步進(jìn)行垂直劃分.為提高垂直劃分的準(zhǔn)確性,以提升后續(xù)比對精度的同時降低后續(xù)比對的時間,本文在FAME算法[25]的基礎(chǔ)上,改進(jìn)了最長兼容鏈構(gòu)建算法以選擇具有最大加權(quán)長度和的公共子鏈集合,并在簇內(nèi)序列垂直分割時引入碎片分組擴(kuò)展,進(jìn)而設(shè)計(jì)簇內(nèi)序列垂直分組方法分別對每個簇Cj進(jìn)行劃分,j=0,1,2,…,H-1.本文給出的簇內(nèi)序列垂直分組方法包括尋找公共子片段、生成公共子鏈、構(gòu)建最長兼容鏈和簇內(nèi)序列垂直分割4個步驟.

2.2.1 尋找公共子片段

尋找公共子片段的目的是找到簇內(nèi)序列的潛在切割點(diǎn).設(shè)簇Cj中包含的長序列數(shù)量為Nj,Cj[u]表示簇Cj中的第u條序列,u=0,1,2,…,Nj-1,j=0,1,2,…,H-1.若有v個k-mer在簇Cj的每條序列中均出現(xiàn),則這些k-mer被稱為公共k-mer,記為km0,km1,…,kmv-1,公共k-mer的長度k=「log2L?,L為簇中序列的平均長度[25].令Npos(Cj[u],kmi)表示kmi在序列Cj[u]中出現(xiàn)的次數(shù),本文采用kmi在簇Cj中的冗余度Redundancy(Cj,kmi)[25]和冗余度閾值Threshold(Cj,kmi)[25]來確定簇Cj中長序列的公共子片段.冗余度Redundancy(Cj,kmi)和冗余度閾值Threshold的計(jì)算分別如式(1)和式(2)所示[25]:

(1)

(2)

簇Cj中序列公共子片段為冗余度小于閾值Thresholdj的公共k-mer.公共k-mer確定了簇內(nèi)序列垂直分割的潛在切割點(diǎn),但冗余度過大的公共k-mer在簇內(nèi)各序列中頻繁出現(xiàn),會導(dǎo)致垂直分割精度下降的同時也加大了計(jì)算開銷.為此,通過冗余度閾值過濾掉冗余度較大的那些公共k-mer.本文尋找簇Cj的公共子片段的過程如下:

第1步.尋找簇Cj中最短序列的下標(biāo)idxmin.

第2步.采用滑動窗口策略提取序列Cj[idxmin]中的k-mer,將k-mer及其在序列Cj[idxmin]中出現(xiàn)的次數(shù)作為鍵值對,存入Hash表Hres.

第3步.u←0.

第4步.若u=idxmin,則轉(zhuǎn)至第6步;否則通過滑動窗口提取序列Cj[u]中的k-mer,若Hres[k-mer]存在,則將此k-mer及其在序列Cj[u]中出現(xiàn)的次數(shù)存入Hash表Hu,并標(biāo)記Hres[k-mer].

第5步.遍歷Hash表Hres,對于Hres中的每個k-mer,若其已被標(biāo)記,則執(zhí)行Hres[k-mer]←Hres[k-mer] *Hu[k-mer];否則從Hres中刪除此k-mer.

第6步.u←u+1,重復(fù)第4～第6步,直至u=Nj為止.

第7步.遍歷Hres,對于Hres中每個k-mer,若Hres[k-mer] >Thresholdj,則從Hres中刪除此k-mer.

第8步.提取Hres中的k-mer作為簇Cj的公共子片段.

2.2.2 生成公共子鏈

確定公共子片段之后,通過將簇中不同序列的相同公共子片段進(jìn)行組合,以生成公共子鏈.簇Cj的公共子鏈的數(shù)目等于簇中各公共子片段的冗余度Thresholdj之和,j=0,1,2,…,H-1.

根據(jù)公共子鏈的位置關(guān)系,本文將公共子鏈分為兼容關(guān)系、相鄰關(guān)系、重疊關(guān)系和交叉關(guān)系.圖1給出了序列公共子鏈位置關(guān)系的示例.

圖1 序列公共子鏈位置關(guān)系示例Fig.1 Example for positional relationship of common subchains in sequences

從圖1可知,相鄰關(guān)系(b)和重疊關(guān)系(c)的公共子鏈之間無間隙Gap且合并時無沖突,可將其合并以降低最長兼容鏈構(gòu)建的時間開銷.合并具有相鄰關(guān)系和重疊關(guān)系的公共子鏈后,剩余的公共子鏈還可能存在兼容關(guān)系(a)和相交關(guān)系(d).若公共子鏈之間存在相交關(guān)系,則無法根據(jù)公共子鏈對簇內(nèi)序列進(jìn)行垂直分割.為此,本文通過構(gòu)建最長兼容鏈,選擇公共子鏈集合的一個兼容子集,以去除具有相交關(guān)系的公共子鏈.

2.2.3 構(gòu)建最長兼容鏈

兼容鏈?zhǔn)枪沧渔湴凑瘴恢门判虼傻囊粭l鏈,在這條鏈中,任意兩個公共子鏈均滿足兼容關(guān)系[25].最長兼容鏈?zhǔn)蔷哂凶疃喙沧渔湹募嫒萱淸25].這個最長兼容鏈的定義將所有公共子鏈平等看待,忽視了公共子鏈寬度對后續(xù)分割精度的影響.

為提高垂直分割的精度,以提升后續(xù)比對的精度和降低比對的時間,本文根據(jù)公共子鏈寬度對公共子鏈進(jìn)行加權(quán),設(shè)計(jì)一個新的算法以選擇具有最大公共子鏈加權(quán)長度和的兼容鏈作為最長兼容鏈.令K表示公共子鏈的數(shù)量,LST表示長度為K的列表,LST中每個元素都保存如下信息:公共子鏈subchain,subchain所在兼容鏈的權(quán)值weight,subchain所在兼容鏈的前驅(qū)節(jié)點(diǎn)在LST中的下標(biāo)previous.算法2描述本文提出的構(gòu)建最長兼容鏈的算法Build-LCC.

算法2.Build-LCC

輸入：含有兼容關(guān)系和相交關(guān)系的公共子鏈列表subchains

輸出：最長兼容鏈LCC

Begin:

1.按前后位置關(guān)系對列表subchains中的公共子鏈進(jìn)行排序;

2.LST←φ;

3.forchaininsubchainsdo

4.weight←公共子鏈chain的寬度;

5. 倒序遍歷LST,找到第i個位置,使得LST[i].chain與chain是兼容關(guān)系;若沒有找到,則i←-1;

6.previous←i;

7.if(i≠-1)then

8.weight←weight+LST[i].weight;

9.endif

10. 從第i個位置開始正向遍歷LST,找到第j個位置,使得LST[j].weight>weight;若沒有找到,則j←Length(LST);//Length(LST)為LST的長度;

11. 將{chain,weight,previous}插入到列表LST的第j個位置;

12.endfor

13.初始化鏈表LCC節(jié)點(diǎn)存儲公共子鏈chain;

14.pos←Length(LST)- 1;

15.while(pos≠-1)do

16. 將LST[pos].chain插入到鏈表LCC的頭部;

17.pos←LST[pos].previous;

18.endwhile

End

算法Build-LCC步驟1需時間O(K×log2K),步驟2～步驟12最壞情形下需時間O(K2),步驟13～步驟18最壞情形下需時間O(K).因此,算法Build-LCC的時間復(fù)雜度為O(K2).LST列表所需空間開銷為O(K),鏈表LCC所需空間開銷為O(K).因此,算法Build-LCC的空間復(fù)雜度為O(K).

2.2.4 簇內(nèi)序列垂直分割

最長兼容鏈構(gòu)建后,其公共子鏈間僅存在兼容關(guān)系,可根據(jù)鏈中的公共子鏈對簇Cj中的序列進(jìn)行垂直分割,將簇Cj劃分為多個垂直分組.圖2給出垂直分割簇Cj的示例.

圖2 垂直分割簇Cj示例Fig.2 Example of splitting cluster Cj vertically

圖2中,簇Cj劃分后的垂直分組分為最長兼容鏈選中的垂直分組和非最長兼容鏈選中的垂直分組.最長兼容鏈選中的垂直分組即公共子鏈,其內(nèi)的各序列片段已對準(zhǔn),無需進(jìn)一步比對.對于非最長兼容鏈選中的垂直分組,需使用已有的適用于規(guī)模較小的多序列比對算法進(jìn)行多序列比對.

對于非最長兼容鏈選中的垂直分組,若分組內(nèi)的序列片段太短,則無法為后續(xù)使用的多序列比對算法提供足夠的殘基信息,這可能導(dǎo)致比對精度下降.本文將此類含有較短序列的非最長兼容鏈選中的垂直分組稱為碎片分組.本文通過對碎片分組進(jìn)行擴(kuò)展,以提升后續(xù)比對的精度.

設(shè)最長兼容鏈中包含的公共子鏈數(shù)量為Length(LCC),Ri表示第i個最長兼容鏈選中的垂直分組,NRj表示第j個非最長兼容鏈選中的垂直分組,Length(Ri)表示Ri中的序列長度,Length(NRi)表示NRi中序列的平均長度,0≤i≤Length(LCC)-1,0≤j≤Length(LCC).設(shè)NRi是碎片分組,圖3給出擴(kuò)展碎片分組NRi的示例.

圖3 對碎片分組NRi擴(kuò)展的示例Fig.3 Example of extending fragmented partition NRi

對于碎片分組NRi兩端的公共子鏈都較窄(圖3(a))的情形,若Length(NRi-1)

經(jīng)簇內(nèi)序列垂直分割和碎片分組擴(kuò)展后,簇Cj被劃分為Vj個垂直分組Cj,0,Cj,1,…,Cj,Vj-1,j=0,1,2,…,H-1,兼容鏈選中的垂直分組內(nèi)的序列已對準(zhǔn),無需后續(xù)比對,但非最長兼容鏈選中的垂直分組內(nèi)的序列仍未對準(zhǔn),需進(jìn)一步使用已有的適用于規(guī)模較小的多序列比對算法將其對準(zhǔn).

2.3 對準(zhǔn)簇內(nèi)分組

若簇Cj中的垂直分組Cj,p為非最長兼容鏈選中的垂直分組,則使用已有的適用于規(guī)模較小的多序列比對算法對此垂直分組內(nèi)的序列進(jìn)行對準(zhǔn),并將比對結(jié)果中的Gap寫入簇Cj內(nèi)序列的相應(yīng)位置.將簇Cj劃分得到的每個非最長兼容鏈選中的垂直分組均對準(zhǔn)后,得到簇Cj內(nèi)多條長序列的完整比對結(jié)果,p=0,1,2,…,Vj-1,j=0,1,2,…,H-1.

2.4 簇間序列水平合并

簇Cj經(jīng)過簇內(nèi)序列垂直分組與簇內(nèi)分組對準(zhǔn)后,簇Cj內(nèi)的那些序列已對準(zhǔn),j=0,1,2,…,H-1,但簇與簇之間的序列尚未對準(zhǔn).為得到序列集S的完整比對結(jié)果,還需進(jìn)行簇間序列水平合并,即以簇為單位將簇C0,C1,…,CH-1進(jìn)行對準(zhǔn).

2.4.1 簇間序列垂直分組

由于對準(zhǔn)過程中添加的間隙Gap的影響,簇Cj的長度LCj大于等于簇Cj中原始序列的最大長度,j=0,1,2,…,H-1.對于大規(guī)模長序列集S,若直接通過比對算法比對各簇,則計(jì)算開銷仍然相當(dāng)大.為此,設(shè)計(jì)簇間序列垂直分組方法,通過尋找簇間公共子片段構(gòu)成的最長兼容鏈,將簇Cj垂直劃分為NG組,得到{Gj,0,Gj,1,…,Gj,NG-1},j=0,1,2,…,H-1,記各簇中的第k組序列組成的集合為簇間分組CGk={G0,k,G1,k,…,GH-1,k},k=0,1,2,…,NG-1,對劃分得到的各個簇間分組分別進(jìn)行對準(zhǔn),以減少簇合并的時間開銷.

為進(jìn)行簇間序列垂直分組,首先將簇Cj映射表示為相同長度LCj的序列Sj,j=0,1,2,…,H-1.對于簇Cj中的每一列,都用序列Sj中的一個字符表示,Sj中第i個字符表示簇Cj中的第i列,若簇Cj第i列中的每個字符都相同,則用簇Cj第i列的字符代替序列Sj中第i個字符,否則用符號‘?’代替序列Sj中的第i個字符,i=0,1,2,…,LCj-1.

獲得簇C0,C1,…,CH-1的代表序列S0,S1,…,SH-1之后,對序列S0,S1,…,SH-1依次進(jìn)行尋找公共子片段、生成公共子鏈、構(gòu)建最長兼容鏈和序列垂直分割.與簇內(nèi)序列垂直分組不同的是,在尋找公共子片段中提取每條序列的k-mer時,需剔除含有符號‘?’的k-mer;在簇間序列垂直切割中劃分Sj時,需在簇Cj的相同列上也進(jìn)行劃分,j=0,1,2,…,H-1.

經(jīng)過簇間序列垂直分組,簇Cj被垂直劃分為NG組Gj,0,Gj,1,…,Gj,NG-1,j=0,1,2,…,H-1,將各簇中第k組進(jìn)行組合得到簇間分組CGk={G0,k,,G1,k,…,GH-1,k},k=0,1,2,…,NG-1.

劃分后的簇間分組包括最長兼容鏈選中的簇間分組和非最長兼容鏈選中的簇間分組兩類,對于最長兼容鏈選中的簇間分組,其內(nèi)的序列已對準(zhǔn),無需進(jìn)一步比對;對于非最長兼容鏈選中的簇間分組,需做進(jìn)一步對準(zhǔn).

2.4.2 簇間分組對準(zhǔn)

簇間垂直分組之后,得到簇間分組CG0,CG1,…,CGNG-1,需依次對準(zhǔn)每個非最長兼容鏈選中的簇間分組以獲得序列集S的完整比對結(jié)果.

為提升簇間分組CGk={G0,k,,G1,k,…,GH-1,k}比對結(jié)果的準(zhǔn)確度,k=0,1,2,…,NG-1,本文設(shè)計(jì)一種新的帶有間隙Gap類型推斷的動態(tài)規(guī)劃漸進(jìn)比對方法.

本文使用基于仿射變換罰分的動態(tài)規(guī)劃漸進(jìn)比對方法對準(zhǔn)分組Gi,k和Gj,k,0 ≤i,j

計(jì)算比對路徑的得分即對準(zhǔn)結(jié)果SGi,j的比對得分AlignScore(i,j)[11],如式(3)所示:

AlignScore(i,j)=

(3)

式(3)中,LSG為對準(zhǔn)結(jié)果SGi,j的列數(shù),NGi、NGj分別表示分組Gi,k和Gj,k中包含的序列數(shù)量,β(Gi,k,b,a)表示對準(zhǔn)后組Gi,k中第b條序列的第a個字符(包括Gap),s(c1,c2)表示字符c1與字符c2的替換分?jǐn)?shù),s(c1,c2)由間隙Gap罰分方法及DNA/RNA或蛋白質(zhì)替換計(jì)分矩陣共同確定.

仿射變換罰分方法將間隙Gap分為4類:開放Gap、擴(kuò)展Gap、終端開放Gap和終端擴(kuò)展Gap,不同Gap類型具有不同的罰分,即不同的s(c1,c2)值[29].現(xiàn)有的多序列比對算法在利用動態(tài)規(guī)劃方法計(jì)算AlignScore(i,j)時,大多不考慮序列分組中已有的Gap與新插入Gap之間類型的相互影響,這導(dǎo)致計(jì)算的比對得分有所偏差.

本文在仿射變換罰分的動態(tài)規(guī)劃漸進(jìn)比對方法基礎(chǔ)上,通過在比對過程中引入Gap類型推斷計(jì)算出準(zhǔn)確的比對得分,以提高分組Gi,k和Gj,k對準(zhǔn)的精度.為進(jìn)行Gap類型推斷,需在比對時獲取各Gap在分組內(nèi)序列中的位置信息、上一比對狀態(tài)信息、當(dāng)前比對狀態(tài)信息和下一比對狀態(tài)信息.

為在比對時快速獲得Gap的位置信息,首先根據(jù)Gap在組內(nèi)序列中的位置,將Gap細(xì)分為8類.圖4給出根據(jù)Gap在組內(nèi)序列中的位置對Gap進(jìn)行分類的示例.

圖4 在組內(nèi)序列中的位置對Gap進(jìn)行分類的示例Fig.4 Example of classifying Gap based on its position in the sequence in group

接著將分組Gi,k轉(zhuǎn)換為采用規(guī)模為(ξalpha+ξgap)×LGi的二維矩陣Pi表示,其中ξalpha表示生物序列字符表的大小(對于DNA或RNA序列ξalpha=4,對于蛋白質(zhì)序列ξalpha=20),ξgap=8表示Gap種類數(shù)為8,矩陣Pi的第m列中,前ξalpha行元素表示組Gi,k第m列中不同字符出現(xiàn)的頻數(shù),后ξgap行元素表示不同Gap類型在組Gi,k第m列中出現(xiàn)的次數(shù).然后,按相同方法將組Gj,k轉(zhuǎn)換為采用規(guī)模為(ξalpha+ξgap)×LGj的二維矩陣Pj表示.以Pi、Pj分別替代分組Gi,k和Gj,k作為動態(tài)規(guī)劃比對方法的輸入.

根據(jù)動態(tài)規(guī)劃矩陣Pi及其當(dāng)前比對列編號a、Pj及其當(dāng)前比對列編號b、上一比對狀態(tài)信息、當(dāng)前比對狀態(tài)信息和下一比對狀態(tài)信息,即可推斷出當(dāng)前比對狀態(tài)對應(yīng)列的Gap的詳細(xì)信息.圖5給出進(jìn)行Gap類型推斷的示例.

在圖5中,若上一比對狀態(tài)為X,當(dāng)前比對狀態(tài)為M,下一比對狀態(tài)為Y,則僅根據(jù)矩陣Pi第a列可得矩陣Pi第a列Gap類型推斷的結(jié)果,僅根據(jù)矩陣Pj第b列可得矩陣Pj第b列Gap類型推斷的結(jié)果.圖5中箭頭表示推斷后各類Gap的數(shù)量與原始Gap數(shù)量之間的關(guān)系.

圖5 比對過程中進(jìn)行Gap類型推斷的示例Fig.5 Example of Gap inference during alignment

由Gap類型推斷的結(jié)果,根據(jù)相應(yīng)的罰分參數(shù)和替換計(jì)分矩陣可得當(dāng)前比對狀態(tài)的得分,將比對路徑中各狀態(tài)節(jié)點(diǎn)的得分相加即為比對路徑的得分.通過帶有間隙Gap類型推斷的動態(tài)規(guī)劃漸進(jìn)比對方法尋找得分最高的比對路徑,并執(zhí)行路徑上各狀態(tài)節(jié)點(diǎn)對應(yīng)的操作即得分組Gi,k和Gj,k的對準(zhǔn)結(jié)果SGi,j.

沿水平聚類分簇得到的引導(dǎo)樹GT,自底向上利用帶有間隙Gap類型推斷的動態(tài)規(guī)劃漸進(jìn)比對方法對準(zhǔn)CGk內(nèi)的兩兩分組后,即得到CGk的比對結(jié)果,k=0,1,2,…,NG-1.

對每個簇間分組CG0,CG1,CG2,…,CGNG-1依次進(jìn)行簇間分組對準(zhǔn)后,即可得到序列集S的完整比對結(jié)果.

2.5 比對算法

為方便理解,圖6給出了本文提出的面向大規(guī)模長序列的多比對算法的處理流程.

圖6 大規(guī)模長序列多比對流程Fig.6 Procedure of multiple alignment for large-scale long reads

算法3描述了本文提出的融合聚類分簇和垂直分組的大規(guī)模長序列多比對算法(簡稱GridMSA算法).

算法3.GridMSA

輸入：長序列集S={S0,S1,…,SN-1}

輸出：帶有Gap的多序列比對結(jié)果S

Begin:

1.fori=0toN-1do

2. 運(yùn)用mBed方法及簡并字母表方法將長序列Si編碼為長度為(log2N)2的數(shù)值向量Di;

3.endfor

4.通過二分k-means算法聚類數(shù)值向量D0,D1,…,DN-1,獲得引導(dǎo)樹GT及序列集S的聚類結(jié)果,將S按序列相似度劃分為H簇C0,C1,…CH-1;

5.forj=0toH-1do

6. 通過簇內(nèi)垂直分組算法尋找簇Cj中以序列公共子片段構(gòu)成的最長兼容鏈,將簇Cj垂直分割為Vj組Cj,0,Cj,1,…,Cj,vj-1;

7.fork=0toVj-1do

8. 若Cj,k為非最長兼容鏈選中的垂直分組,則執(zhí)行現(xiàn)有的適用于規(guī)模較小的多序列比對算法對準(zhǔn)Cj,k內(nèi)序列,并將比對結(jié)果中的Gap添加到S中原序列的對應(yīng)位置,以獲得簇Cj內(nèi)各序列的比對結(jié)果;

9.endfor

10.endfor

11.通過簇間序列垂直分組算法尋找簇C0,C1,…CH-1之間公共子片段構(gòu)成的最長兼容鏈,將序列集S劃分為NG個簇間分組CG0,CG1,…,CGNG-1;

12.fork=0toNG-1do

13. 若CGk為非最長兼容鏈選中的簇間分組,則執(zhí)行動態(tài)規(guī)劃漸進(jìn)比對方法對準(zhǔn)簇間分組CGk內(nèi)的各組G0,k,G1,k,…GH-1,k,并將比對結(jié)果中的Gap添加到S中原序列的對應(yīng)位置;

14.endfor

End

1http://github.com/naznoosh/msa

3 實(shí) 驗(yàn)

3.1 實(shí)驗(yàn)環(huán)境與數(shù)據(jù)

實(shí)驗(yàn)在廣西大學(xué)Sugon 7000A超級并行計(jì)算機(jī)系統(tǒng)進(jìn)行.使用節(jié)點(diǎn)的處理器為Intel Xeon Gold 6230,CPU主頻2.1GHz,動態(tài)加速頻率3.9GHz,節(jié)點(diǎn)的DDR4內(nèi)存頻率為2933MHz,內(nèi)存容量512GB,硬盤容量8×900GB.運(yùn)行64位CentOS 7.4操作系統(tǒng).采用C語言編程實(shí)現(xiàn)算法.

表1 實(shí)驗(yàn)使用的長序列數(shù)據(jù)集信息Table 1 Information of experimental long-read datasets used

本文同樣采用文獻(xiàn)[25]使用的數(shù)據(jù)集進(jìn)行實(shí)驗(yàn),該數(shù)據(jù)集可從Github網(wǎng)站1免費(fèi)獲得.數(shù)據(jù)集中包含9個不同的長序列集,各序列集的信息如表1所示.注:Github網(wǎng)站上未提供文獻(xiàn)[25]提到的超長序列集Sorangium cellulosum,因此表1中沒有包含此數(shù)據(jù)集.

3.2 結(jié)果實(shí)驗(yàn)

實(shí)驗(yàn)測試了本文算法GridMSA和同類算法MAFFT[13]、ClustalO[7]、ClustalW[5]、FAMSA[8]、FAME[25]在不同序列集上運(yùn)行的比對結(jié)果質(zhì)量(SPavg)[25]和時間(Time).

設(shè)序列集S對準(zhǔn)后每條序列的長度為L,Si,j表示比對完成后S中第i條序列的第j個字符,0≤i

(4)

本文算法GridMSA在水平聚類分簇時通過k-mer編碼將字符序列轉(zhuǎn)換為數(shù)值向量,以計(jì)算序列間的距離,參照文獻(xiàn)[28]選取k1=5對DNA序列和RNA序列進(jìn)行編碼、選取k1=2對蛋白質(zhì)序列進(jìn)行編碼.算法GridMSA在垂直分組時需依據(jù)k-mer尋找長序列的公共子片段,為此參照文獻(xiàn)[25]設(shè)置k2=「log2Lavg?,Lavg為S中序列的平均長度.本文使用與文獻(xiàn)[28]相同的替換計(jì)分矩陣.為設(shè)置空隙Gap罰分參數(shù),將本文8種不同Gap類型映射到文獻(xiàn)[28]中4種不同Gap類型,對應(yīng)關(guān)系如表2所示.

表2 Gap類型映射及罰分取值Table 2 Gap type mapping and penalty setting

本文將實(shí)現(xiàn)多序列比對算法FAME時分別執(zhí)行適用于規(guī)模較小的多序列比對算法MAFFT、ClustalO、ClustalW、FAMSA完成對垂直分組得到各分組中的序列進(jìn)行比對的版本分別記為FAME1、FAME2、FAME3和FAME4.同時,將實(shí)現(xiàn)多序列比對算法GridMSA時分別執(zhí)行適用于規(guī)模較小的多序列比對算法MAFFT、ClustalO、ClustalW、FAMSA完成對垂直分組得到各分組中的序列進(jìn)行比對的版本分別記為GridMSA1、GridMSA2、GridMSA3和GridMSA4.

實(shí)驗(yàn)首先測試了算法在大規(guī)模長序列數(shù)據(jù)集Titin、Variola、Mycoplasma、Streptococcus和Ecoli上的運(yùn)行時間(Time)和比對結(jié)果質(zhì)量(SPavg),結(jié)果如表3所示.

根據(jù)表3的實(shí)驗(yàn)結(jié)果,在長序列數(shù)據(jù)集Titin、Variola、Mycoplasma、Streptococcus和Ecoli上,本文算法實(shí)驗(yàn)結(jié)果得到的SPavg值大部分比其他同類算法實(shí)驗(yàn)得到的SPavg值更接近0,這表明本文算法GridMSA實(shí)現(xiàn)的4個版本GridMSA1、GridMSA2、GridMSA3和GridMSA4在進(jìn)行比對時均具有整體上較好的比對結(jié)果質(zhì)量和較低的運(yùn)行時間開銷,其中在超長序列集Streptococcus和Ecoli上運(yùn)行的效果優(yōu)勢更為明顯.這主要是因?yàn)樗惴℅ridMSA在垂直分組時選擇具有最大公共子鏈長度和的兼容鏈作為最長兼容鏈,一方面使得劃分后的待比對區(qū)域長度減小,從而降低了比對時間,另一方面最長兼容鏈選擇區(qū)域更加合理,從而提升了比對精度;此外算法GridMSA將垂直分割后的碎片分組擴(kuò)展,為后續(xù)多序列比對提供了足夠的殘基信息,這進(jìn)一步提升了比對精度.

為評估算法在不同規(guī)模的同一個數(shù)據(jù)集上運(yùn)行的比對效果,實(shí)驗(yàn)進(jìn)一步測試了算法在數(shù)據(jù)集MT-1X、MT-20X、MT-50X、MT-100X上的運(yùn)行時間(Time)和比對結(jié)果質(zhì)量(SPavg),結(jié)果如表4所示.

從表4的實(shí)驗(yàn)結(jié)果可以看到:在大規(guī)模長序列數(shù)據(jù)集MT-1X～MT-100X上,與其他算法相比,本文算法GridMSA實(shí)現(xiàn)的4個版本GridMSA1、GridMSA2、GridMSA3和GridMSA4進(jìn)行多序列比對在整體上保持比對精度的同時,所需的時間開銷更少.當(dāng)序列集數(shù)據(jù)規(guī)模越大時,與其他算法相比,本文算法GridMSA進(jìn)行多序列比對獲得的加速效果越為明顯.這主要是因?yàn)樗惴℅ridMSA通過融合水平聚類分簇和垂直分組,將原始大規(guī)模長序列比對問題分解為若干比對子問題進(jìn)行求解,從而降低了比對時間;同時算法GridMSA在簇間分組對準(zhǔn)時考慮到簇中已有Gap與新插入Gap之間類型的相互影響,在動態(tài)規(guī)劃漸進(jìn)比對過程中引入Gap類型推斷計(jì)算出了準(zhǔn)確的比對得分,從而確保了比對的精度.雖然算法FAME在大規(guī)模序列數(shù)據(jù)集上運(yùn)行也具有較好的比對結(jié)果質(zhì)量和較低的時間開銷,但由于算法FAME僅對序列集進(jìn)行垂直分割,劃分后各子序列集的序列數(shù)量并沒有減少,它對垂直分割得到的各子序列集中的序列進(jìn)行多比對較為耗時,所以其所需的比對時間高于本文算法GridMSA.

表3 算法在5種長序列集上運(yùn)行的時間和比對結(jié)果質(zhì)量Table 3 Required time and values of SPavg of running the algorithms on five long-read sets

綜合表3和表4的實(shí)驗(yàn)結(jié)果可知:采用本文算法GridMSA對大規(guī)模長序列數(shù)據(jù)集進(jìn)行多比對,整體上可以大大降低比對時間且獲得較好的比對結(jié)果質(zhì)量.序列集規(guī)模越大、序列長度越長,算法GridMSA加速效果優(yōu)勢越為明顯.算法GridMSA在垂直分組時選擇具有最大公共子鏈長度和的兼容鏈作為最長兼容鏈并對垂直分割后的碎片分組進(jìn)行擴(kuò)展,同時在簇間序列水平合并過程中引入間隙Gap類型推斷計(jì)算出準(zhǔn)確的對準(zhǔn)得分,因此算法GridMSA對大規(guī)模的長序列集進(jìn)行多比對,仍可以獲得較高的比對結(jié)果質(zhì)量得分.另一方面,在9組數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果表明,當(dāng)長序列集數(shù)據(jù)規(guī)模相對較小時,算法GridMSA執(zhí)行適用于規(guī)模較小的多比對算法MAFFT對垂直分組得到各分組中的序列進(jìn)行比對時,獲得整體上較好的多序列比對結(jié)果;當(dāng)長序列集數(shù)據(jù)規(guī)模較大時,算法GridMSA執(zhí)行適用于規(guī)模較小的多比對算法FAMSA對垂直分組得到各分組中的序列進(jìn)行比對時,獲得整體上較好的多序列比對結(jié)果.

4 總 結(jié)

本文的研究特色和新穎點(diǎn)是:基于劃分思想,通過融合水平聚類分簇和垂直分組將大規(guī)模長序列集劃分成多個規(guī)模相對較小、長度較短的子序列集;通過在構(gòu)建最長兼容鏈時對公共子鏈按其寬度加權(quán)、在簇內(nèi)序列垂直分割時進(jìn)行碎片分組擴(kuò)展,進(jìn)而設(shè)計(jì)了簇內(nèi)序列垂直分組方法,利用此方法垂直分割水平簇以提升后續(xù)比對的質(zhì)量和速度;提出針對水平簇集的簇間序列分組算法,在簇間序列水平合并時利用簇間序列分組算法將大規(guī)模長序列數(shù)據(jù)集垂直劃分為多個簇間分組分別進(jìn)行比對,以加快比對速度;利用水平簇內(nèi)序列多比對結(jié)果中包含的Gap信息,設(shè)計(jì)了一種新的帶有Gap類型推斷的動態(tài)規(guī)劃漸進(jìn)比對方法,利用此方法對劃分得到的各簇間分組分別進(jìn)行多比對,可以確保大規(guī)模長序列比對結(jié)果保持較高的精度.在大規(guī)模長序列數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果表明,本文算法GridMSA在整體上維持較高的多比對結(jié)果質(zhì)量的同時,大大地減少了大規(guī)模長序列多比對的計(jì)算時間;當(dāng)數(shù)據(jù)規(guī)模越大、序列長度越長的時候,本文算法加速效果越明顯.本文算法GridMSA的運(yùn)行時間和待比對的序列長度以及數(shù)據(jù)規(guī)模成正相關(guān),當(dāng)數(shù)據(jù)集規(guī)模更大、序列更長的時候,算法GridMSA仍需要較長的時間開銷.如何通過CPU+GPU協(xié)同并行計(jì)算加速算法GridMSA以適應(yīng)更長序列更大規(guī)模的多序列比對問題將是下一步研究方向.