水平聚類分簇和垂直分組的大規(guī)模長(zhǎng)序列多比對(duì)

王 淋,鐘 誠(chéng)

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

1 引 言

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

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

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

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

2 方 法

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

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

2.1 多序列聚類分簇

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

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

算法1.MS-kCluster

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

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

Begin:

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

2.fori=0toN-1do

3. 通過(guò)簡(jiǎn)并字母表方法計(jì)算長(zhǎng)序列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對(duì)應(yīng)的數(shù)值向量Di=(D(Si,M0),D(Si,M1),…,D(Si,Mm-1));

10.endfor

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

12.初始化引導(dǎo)樹(shù)GT←node且樹(shù)當(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ù)值向量集合分割為兩個(gè)子集合Dleft、Dright;

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

17.node.right←Node(Dright);//node.right表示引導(dǎo)樹(shù)節(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ù)值向量集對(duì)應(yīng)的序列集合};

23.endfor

End

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

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

2.2.1 尋找公共子片段

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

(1)

(2)

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

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

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

第3步.u←0.

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

第5步.遍歷Hash表Hres,對(duì)于Hres中的每個(gè)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,對(duì)于Hres中每個(gè)k-mer,若Hres[k-mer] >Thresholdj,則從Hres中刪除此k-mer.

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

2.2.2 生成公共子鏈

確定公共子片段之后,通過(guò)將簇中不同序列的相同公共子片段進(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)的公共子鏈之間無(wú)間隙Gap且合并時(shí)無(wú)沖突,可將其合并以降低最長(zhǎng)兼容鏈構(gòu)建的時(shí)間開(kāi)銷.合并具有相鄰關(guān)系和重疊關(guān)系的公共子鏈后,剩余的公共子鏈還可能存在兼容關(guān)系(a)和相交關(guān)系(d).若公共子鏈之間存在相交關(guān)系,則無(wú)法根據(jù)公共子鏈對(duì)簇內(nèi)序列進(jìn)行垂直分割.為此,本文通過(guò)構(gòu)建最長(zhǎng)兼容鏈,選擇公共子鏈集合的一個(gè)兼容子集,以去除具有相交關(guān)系的公共子鏈.

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

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

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

算法2.Build-LCC

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

輸出：最長(zhǎng)兼容鏈LCC

Begin:

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

2.LST←φ;

3.forchaininsubchainsdo

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

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

6.previous←i;

7.if(i≠-1)then

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

9.endif

10. 從第i個(gè)位置開(kāi)始正向遍歷LST,找到第j個(gè)位置,使得LST[j].weight>weight;若沒(méi)有找到,則j←Length(LST);//Length(LST)為L(zhǎng)ST的長(zhǎng)度;

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

12.endfor

13.初始化鏈表LCC節(jié)點(diǎn)存儲(chǔ)公共子鏈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需時(shí)間O(K×log2K),步驟2～步驟12最壞情形下需時(shí)間O(K2),步驟13～步驟18最壞情形下需時(shí)間O(K).因此,算法Build-LCC的時(shí)間復(fù)雜度為O(K2).LST列表所需空間開(kāi)銷為O(K),鏈表LCC所需空間開(kāi)銷為O(K).因此,算法Build-LCC的空間復(fù)雜度為O(K).

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

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

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

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

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

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

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

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

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

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

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

2.4 簇間序列水平合并

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

2.4.1 簇間序列垂直分組

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

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

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

經(jīng)過(guò)簇間序列垂直分組,簇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.

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

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

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

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

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

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

AlignScore(i,j)=

(3)

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

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

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

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

圖4 在組內(nèi)序列中的位置對(duì)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表示生物序列字符表的大小(對(duì)于DNA或RNA序列ξalpha=4,對(duì)于蛋白質(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作為動(dòng)態(tài)規(guī)劃比對(duì)方法的輸入.

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

在圖5中,若上一比對(duì)狀態(tài)為X,當(dāng)前比對(duì)狀態(tài)為M,下一比對(duì)狀態(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 比對(duì)過(guò)程中進(jìn)行Gap類型推斷的示例Fig.5 Example of Gap inference during alignment

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

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

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

2.5 比對(duì)算法

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

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

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

算法3.GridMSA

輸入：長(zhǎng)序列集S={S0,S1,…,SN-1}

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

Begin:

1.fori=0toN-1do

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

3.endfor

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

5.forj=0toH-1do

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

7.fork=0toVj-1do

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

9.endfor

10.endfor

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

12.fork=0toNG-1do

13. 若CGk為非最長(zhǎng)兼容鏈選中的簇間分組,則執(zhí)行動(dòng)態(tài)規(guī)劃漸進(jìn)比對(duì)方法對(duì)準(zhǔn)簇間分組CGk內(nèi)的各組G0,k,G1,k,…GH-1,k,并將比對(duì)結(jié)果中的Gap添加到S中原序列的對(duì)應(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ì)算機(jī)系統(tǒng)進(jìn)行.使用節(jié)點(diǎn)的處理器為Intel Xeon Gold 6230,CPU主頻2.1GHz,動(dòng)態(tài)加速頻率3.9GHz,節(jié)點(diǎn)的DDR4內(nèi)存頻率為2933MHz,內(nèi)存容量512GB,硬盤(pán)容量8×900GB.運(yùn)行64位CentOS 7.4操作系統(tǒng).采用C語(yǔ)言編程實(shí)現(xiàn)算法.

表1 實(shí)驗(yàn)使用的長(zhǎng)序列數(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個(gè)不同的長(zhǎng)序列集,各序列集的信息如表1所示.注:Github網(wǎng)站上未提供文獻(xiàn)[25]提到的超長(zhǎng)序列集Sorangium cellulosum,因此表1中沒(méi)有包含此數(shù)據(jù)集.

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

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

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

(4)

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

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

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

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

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

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

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

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

4 總 結(jié)

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