基于剖面隱馬氏模型的多序列比對(duì)

李成淵, 龍海俠, 孫俊, 須文波＊

(1.江南大學(xué)信息工程學(xué)院,江蘇無(wú)錫 214122;2.江南大學(xué)教育學(xué)院,江蘇無(wú)錫 214122)

基于剖面隱馬氏模型的多序列比對(duì)

李成淵1, 龍海俠2, 孫俊1, 須文波＊1

(1.江南大學(xué)信息工程學(xué)院,江蘇無(wú)錫 214122;2.江南大學(xué)教育學(xué)院,江蘇無(wú)錫 214122)

多序列比對(duì)被稱為NP完全問(wèn)題,是生物信息中最基本的問(wèn)題之一。目前,廣泛使用剖面隱馬爾可夫模型解決多序列比對(duì)問(wèn)題。作者在粒子群優(yōu)化算法的基礎(chǔ)上,提出了將量子粒子群優(yōu)化算法用于剖面隱馬爾可夫模型的訓(xùn)練過(guò)程,進(jìn)而構(gòu)建了一種基于剖面隱馬氏模型和量子粒子群優(yōu)化算法的多序列比對(duì)算法。從核酸序列和BaliBASE比對(duì)數(shù)據(jù)庫(kù)中選取了一些比對(duì)例子進(jìn)行了模擬實(shí)驗(yàn),并與其他算法進(jìn)行了比較,結(jié)果表明,所提出的算法能在有限的時(shí)間內(nèi)不僅能找到理想的隱隱馬爾可夫模型,而且能得到最優(yōu)的比對(duì)結(jié)果。

多序列比對(duì);剖面隱馬爾可夫模型;量子粒子群優(yōu)化算法

核苷酸或氨基酸的多序列比對(duì)或聯(lián)配是生物信息學(xué)中最重要、最具有挑戰(zhàn)性的任務(wù)之一。多序列比對(duì)問(wèn)題是一個(gè)將不等長(zhǎng)的多個(gè)序列通過(guò)插入空位變成等長(zhǎng)的過(guò)程,這些位置上的空位代表著相比對(duì)的序列從共同的祖先通過(guò)插入/刪除操作的進(jìn)化過(guò)程。利用多序列比對(duì)算法得到的最優(yōu)比對(duì),可用于找出蛋白質(zhì)家族的模體(motifs)或保守區(qū)域(conserved domains),可用于預(yù)測(cè)蛋白質(zhì)的結(jié)構(gòu)和功能,也可用于進(jìn)行系統(tǒng)發(fā)育的分析[1-2]。

目前主要有下列3種策略用于多序列比對(duì)。第一種策略是“漸進(jìn)比對(duì)”策略[3-4],其基本思想是:迭代地利用兩序列動(dòng)態(tài)規(guī)劃算法,先由兩條序列的比對(duì)開(kāi)始,逐漸添加新序列,直到所有序列都加入為止。第二種策略是使用隨機(jī)優(yōu)化算法,如模擬退火算法(SA)[5],遺傳算法(GA)[6];第三種策略基于概率模型的隱馬爾可夫模型[7-8]。

本研究中使用第三種策略。在多序列比對(duì)的過(guò)程中,隱馬爾可夫模型主要解決3個(gè)問(wèn)題:一是得分問(wèn)題,二是聯(lián)配問(wèn)題,三是訓(xùn)練問(wèn)題。將得分問(wèn)題用來(lái)評(píng)估模型的性能,聯(lián)配問(wèn)題用來(lái)實(shí)現(xiàn)多序列的比對(duì),訓(xùn)練問(wèn)題用來(lái)優(yōu)化模型的參數(shù)。最常用的訓(xùn)練隱馬爾可夫模型模型的方法是基于統(tǒng)計(jì)和重估的方法,比如Baum-Welch算法。但是Baum-Welch算法是一個(gè)局部最優(yōu)算法,使用此算法得到的最終比對(duì)結(jié)果通常遠(yuǎn)離全局最優(yōu)。最近還出現(xiàn)了粒子群算法(PSO),此算法也是一個(gè)局部最優(yōu)算法。

為了克服Baum-Welch算法和PSO算法的缺點(diǎn),我們使用量子粒子群優(yōu)化算法(QPSO)[9-10]來(lái)訓(xùn)練隱馬爾可夫模型,不僅參數(shù)個(gè)數(shù)少,隨機(jī)性強(qiáng),并且能覆蓋所有解空間,保證算法的全局收斂。

1 剖面隱馬爾可夫模型的拓?fù)浣Y(jié)構(gòu)

圖1 用于多序列比對(duì)的隱馬爾可夫模型Fig.1 An example of a simple HMM of length 3 for MSA

當(dāng)使用圖1所示的隱馬爾可夫模型進(jìn)行多序列比對(duì)時(shí),每條序列從開(kāi)始到結(jié)束通過(guò)這些狀態(tài)穿越模型,在這些待比對(duì)序列中進(jìn)行空位字符‘-’的插入和刪除操作,得到一個(gè)多序列比對(duì)結(jié)果的矩陣A=(aij)m×n,其中aij∈alph_set∪{-}。矩陣A中的每一列為一個(gè)位點(diǎn)上的比對(duì),矩陣A的第i行對(duì)應(yīng)參與比對(duì)的第i個(gè)序列,序列中非空字符的先后順序在比對(duì)中保持不變。

2 基于剖面隱馬爾可夫模型和QPSO的多序列比對(duì)

用量子粒子群優(yōu)化算法訓(xùn)練剖面隱馬爾可夫模型時(shí),每一個(gè)粒子代表一個(gè)隱馬爾可夫模型,通過(guò)不斷的更新粒子的位置來(lái)優(yōu)化隱馬爾可夫模型。在訓(xùn)練中保持模型的長(zhǎng)度不變,僅僅優(yōu)化模型的參數(shù):轉(zhuǎn)移概率和符號(hào)發(fā)出概率。對(duì)圖1所示的隱馬爾可夫模型的拓?fù)浣Y(jié)構(gòu),我們?nèi)〈葘?duì)序列的平均值m為模型的長(zhǎng)度,不考慮初始狀態(tài)和結(jié)束狀態(tài),則模型的狀態(tài)數(shù)為3m+1,狀態(tài)轉(zhuǎn)移概率參數(shù)為3(3m+1)個(gè);設(shè)字符集大小為|A|,共有(2m+1) |A|個(gè)符號(hào)發(fā)出概率。所以每個(gè)粒子是維數(shù)為9m+3+(2m+1)|A|的一個(gè)實(shí)數(shù)編碼串。所以DN A模型的參數(shù)個(gè)數(shù)是17m+7,蛋白質(zhì)模型的參數(shù)個(gè)數(shù)是49m+23。根據(jù)轉(zhuǎn)移概率和符號(hào)發(fā)出概率的性質(zhì),在對(duì)粒子對(duì)應(yīng)的隱馬爾可夫模型模型進(jìn)行評(píng)價(jià)前,需要先對(duì)隱馬爾可夫模型中的狀態(tài)轉(zhuǎn)移概率和符號(hào)發(fā)出概率進(jìn)行歸一化,以滿足3m+1個(gè)轉(zhuǎn)移概率歸一化約束方程和2m+1個(gè)符號(hào)發(fā)出概率歸一化方程。

根據(jù)粒子群算法訓(xùn)練的結(jié)果,得到全局最優(yōu)的粒子對(duì)應(yīng)的隱馬爾可夫模型,接下來(lái)用此模型使用Viterbi算法進(jìn)行序列的比對(duì),得到最優(yōu)的比對(duì)結(jié)果,并用基于SP(sum of pairs)打分系統(tǒng)的目標(biāo)函數(shù)評(píng)估比對(duì)的結(jié)果。

在整個(gè)算法的過(guò)程中,根據(jù)剖面隱馬爾可夫模型的拓?fù)浣Y(jié)構(gòu),所使用的轉(zhuǎn)移概率的形式見(jiàn)表1所示。

表1 轉(zhuǎn)移概率Tab.1 The transition probability

2.1 模型的訓(xùn)練問(wèn)題

2.1.1 量子粒子群優(yōu)化算法(QPSO) PSO[11]是基于種群的進(jìn)化搜索技術(shù),但是所有基本和改進(jìn)的PSO算法不能保證算法的全局收斂[12]。因?yàn)镻SO的進(jìn)化方程式使所有粒子在一個(gè)有限的樣本空間中搜索。根據(jù)粒子群的基本收斂性質(zhì),受量子物理基本理論的啟發(fā),Sun等人提出了QPSO(Quantum-behaved Particle Swarm Optimization[9-10]算法是對(duì)整個(gè)PSO算法進(jìn)化搜索策略的改變,進(jìn)化方程中不需要速度向量,形式更簡(jiǎn)單,參數(shù)更少且更容易控制。

在QPSO算法中,粒子按照下面3個(gè)公式進(jìn)行更新:

2.1.2 評(píng)估訓(xùn)練算法的質(zhì)量 在粒子群優(yōu)化算法中,需要評(píng)估每個(gè)粒子所代表的模型的質(zhì)量,使用的評(píng)估函數(shù)為:

這里O={O1,O2,…,OM}是給定的待比對(duì)序列的集合,序列的個(gè)數(shù)為M個(gè)。lS是序列OS的長(zhǎng)度。Log_odd的值越大,說(shuō)明使用量子粒子群優(yōu)化算法訓(xùn)練得到的隱馬爾可夫模型,模型的參數(shù)是最優(yōu)的,模型的穩(wěn)定性和可靠性都較好。

2.2 模型的聯(lián)配問(wèn)題

2.2.1 序列聯(lián)配的過(guò)程 給定3條待比對(duì)序列:O (1):A GGCT;O(2):GAACTGTA;O(3): AGCCTTA。按照下面的兩個(gè)步驟可以得到序列比對(duì)的結(jié)果:

1)根據(jù)Viterbi算法和圖1所示的隱馬爾可夫模型拓?fù)浣Y(jié)構(gòu)找出每條序列所對(duì)應(yīng)的狀態(tài)序列,如圖2所示,每條序列的氨基酸堿基對(duì)應(yīng)一個(gè)匹配狀態(tài)或一個(gè)插入狀態(tài)。

圖2 每條序列所對(duì)應(yīng)的狀態(tài)序列Fig.2 Each sequence corresponding to the state sequence

由圖3生成的狀態(tài)序列,我們可以進(jìn)行空位字符‘—’的插入操作,見(jiàn)圖3。所有的序列中與匹配狀態(tài)M相對(duì)應(yīng)的氨基酸堿基是比對(duì)的,這些氨基酸堿基位于同一列;與插入狀態(tài)Ik相對(duì)應(yīng)的氨基酸堿基位于Mk-1或DK-1之后,Mk或DK之前。

圖3 空位字符‘-’的插入Fig.3 Insert the gap characters’-’

由步驟(1)和步驟(2),得到聯(lián)配后的序列為圖4所示。

圖4 聯(lián)配后的序列Fig.4 Aligned sequences

2.2.2 評(píng)估比對(duì)序列的質(zhì)量在利用Viterrbi算法獲得的路徑進(jìn)行比對(duì)后,需要通過(guò)基于SP(sum of pairs)打分系統(tǒng)的目標(biāo)函數(shù)[13]對(duì)比對(duì)結(jié)果進(jìn)行評(píng)估。我們使用下面標(biāo)準(zhǔn)的sum-of-pairs打分函數(shù):

這里,li是已比對(duì)的序列,D為距離矩陣.

為了避免在比對(duì)過(guò)程中空位的積聚,我們從SOP分?jǐn)?shù)中推演出仿射幾何學(xué)的空位代價(jià),對(duì)于比對(duì)結(jié)果中一條序列的空位代價(jià)按照下面的公式進(jìn)行計(jì)算:

這里GOP表示第一個(gè)開(kāi)口空位的固定罰分, GEP表示對(duì)于擴(kuò)展的空位的罰分,n為一條序列中空位的個(gè)數(shù)。對(duì)于已比對(duì)的每條序列的空位都要計(jì)算相應(yīng)的空位代價(jià)。多序列比對(duì)結(jié)果的SOP的分值減去空位代價(jià)的總和,即為SOP的分值。

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

3.1 實(shí)驗(yàn)數(shù)據(jù)

3.1.1 數(shù)據(jù)集1(模擬的核苷酸序列) 使用軟件Rose[8]軟件產(chǎn)生如下的核酸序列:“l(fā)ow-short”“l(fā)ow-long”,“high-short”,和“high-long”.

3.1.2 數(shù)據(jù)集2(Pfam數(shù)據(jù)庫(kù)中的3個(gè)蛋白質(zhì)家族) 3個(gè)蛋白質(zhì)家族分別為G5,CagY_M and Interferon,它們來(lái)自Pfam數(shù)據(jù)庫(kù)[14],見(jiàn)表2。為了避免隱馬爾可夫模型模型中的過(guò)擬合,分別把3個(gè)蛋白質(zhì)家族分成訓(xùn)練集和測(cè)試集。

表2 蛋白質(zhì)家族N:序列的個(gè)數(shù),LSEQ:數(shù)列的平均,T:訓(xùn)練集的大小Tab.2 Protein familiesN:the number of sequences LSEQ:the average of seriesT:The size of training sets

3.2 實(shí)驗(yàn)設(shè)置

實(shí)驗(yàn)中,我們使用均勻分布初始化所有的種群;分別使用了Baum-Welch(BW)、Particle Swarm Optimization(PSO)、Quantum-behavedParticle Swarm Optimization(QPSO)3種算法訓(xùn)練隱馬爾可夫模型,用Log_odd目標(biāo)函數(shù)來(lái)評(píng)估4種訓(xùn)練算法的性能,并且使用SOP打分系統(tǒng)的目標(biāo)函數(shù)來(lái)評(píng)估比對(duì)效果。每種算法分別迭代1 000次,共進(jìn)行了20次模擬實(shí)驗(yàn)。

對(duì)于PSO算法,使用如下的參數(shù):種群數(shù)=20;慣性權(quán)重(ω)從1.0到0.5線性減少;c1=c2=2.0最大速度(v→max=1.0

對(duì)于QPSO算法使用如下的參數(shù):種群大小= 20;收縮擴(kuò)張系數(shù),從1.0到0.5線性減少;并且,在SOP打分函數(shù)中,核酸和蛋白質(zhì)分別使用下面的打分矩陣、空位開(kāi)放和空位擴(kuò)充。

1)對(duì)于核酸序列,使用ClustalW 1.81的“swgap”置換分?jǐn)?shù)表,空位開(kāi)放和空位擴(kuò)充分別為15、7。

2)對(duì)于蛋白質(zhì)序列,使用BLOSUM62打分矩陣,空位開(kāi)放和空位擴(kuò)充分別為11、2[9]。

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

表3給出了核酸序列的結(jié)果,這里是以Logodd score作為目標(biāo)函數(shù),用來(lái)評(píng)估訓(xùn)練隱馬爾可夫模型的質(zhì)量,從表中可以看出,QPSO取得了最好的平均值,PSO算法的結(jié)果和BW算法的結(jié)果相當(dāng),在Low-long序列和High-long序列上的結(jié)果不如BW得出的結(jié)果好,可以說(shuō)明,PSO算法隨著序列個(gè)數(shù)的增加效果不如BW的效果好。

表4也給出了核酸序列的結(jié)果,這里是以SOP score作為目標(biāo)函數(shù),用來(lái)評(píng)估核酸序列的比對(duì)的效果。從表4可以看出,QPSO優(yōu)于其他所有的算法,其次為CW工具得出的比對(duì)結(jié)果優(yōu)于BW和PSO算法得出的結(jié)果,PSO算法也是一種局部最優(yōu)算法,對(duì)BW算法沒(méi)有多大的改進(jìn)。

表3 隱馬爾可夫模型得到的核酸序列的Log_odd scores平均值和方差Tab.3 HMM log-odds scores±standard error of Nucleotidesequences

表5～8概述了3個(gè)蛋白質(zhì)家族的隱馬爾可夫模型訓(xùn)練和序列比對(duì)的結(jié)果。表5和表6分別給出了訓(xùn)練集的最優(yōu)的Log-odd scores和SOP scores的平均值和方差,表7和表8分別給出了評(píng)估集的最優(yōu)的Log-odd scores和SOP scores的平均值和方差。從4個(gè)表中,可以看出,對(duì)于蛋白質(zhì)家族,無(wú)論是在訓(xùn)練集還是評(píng)估集上,QPSO算法在模型的訓(xùn)練上和比對(duì)的效果上都優(yōu)于BW和PSO算法。

表4 核酸序列比對(duì)結(jié)果的SOP平均值和方差Tab.4 SOP scores for the final alignments of Nucleotide sequences

表5 隱馬爾可夫模型得到的蛋白質(zhì)序列訓(xùn)練集的Log_odd scores平均值和方差Tab.5 HMM log-odds scores±standard error for the training sets of three protein families

表6 蛋白質(zhì)序列訓(xùn)練集的最優(yōu)的SOP scores的平均值和方差Tab.6 SOP scores for the final alignments of the training sets of three protein families

表7 隱馬爾可夫模型得到的蛋白質(zhì)序列測(cè)試集的Log_odd scores平均值和方差Tab.7 HMM log-odds scores±standard error for the validation sets of three protein families

表8 蛋白質(zhì)序列測(cè)試集的最優(yōu)的SOP scores的平均值和方差Tab.8 SOP scores for the final alignments of the validation sets of three protein families

圖5、6、7刻畫了以Log-odd score作為目標(biāo)函數(shù),算法分別運(yùn)行20次得到的核酸序列和蛋白質(zhì)序列的平均值的收斂過(guò)程。

圖5 Low-short核酸序列的Log-odd平均分?jǐn)?shù)Fig.5 Mean Log-odds scores for Low-short Nucleotide

圖6 Cag_Y蛋白質(zhì)序列的Log_odd平均分?jǐn)?shù)Fig.6 Mean Log-odds scores for training sets of CagY_ M protein family

圖7 Interferon蛋白質(zhì)序列的Log_odd平均分?jǐn)?shù)Fig.7 Mean Log-odds scores for validation sets of Interferon protein family

從圖中看出,QPSO性能最好,在迭代的最后得到的Log-odd score的分值也最好,除了圖9, PSO算法的性能優(yōu)于BW算法。從收斂速度上看, BW和PSO的收斂速度最快,其次為QPSO算法,但是QPSO算法在迭代進(jìn)行到一半的時(shí)候,Logodd score的值就已經(jīng)超過(guò)PSO算法和BW算法,并且在整個(gè)運(yùn)行過(guò)程中,QPSO都在不斷的提高自身的性能,而PSO算法和BW算法已經(jīng)早早地收斂了。

圖8 Low-short核酸序列的SOP平均分?jǐn)?shù)Fig.8 Mean SOP scores for Low-short Nucleotide

圖9 Cag_Y蛋白質(zhì)序列的SOP平均分?jǐn)?shù)Fig.9 Mean SOP scores for training sets of CagY_M protein family

圖10 Interferon蛋白質(zhì)序列的SOP平均分?jǐn)?shù)Fig.10 Mean SOP scores for validation sets of Interferon protein family

圖8、9、10表明了以SOP score作為目標(biāo)函數(shù),算法分別運(yùn)行20次得到的核酸序列和蛋白質(zhì)序列的平均值的收斂過(guò)程。收斂情況與上述情況相同。

4 結(jié) 語(yǔ)

提出了使用QPSO算法來(lái)訓(xùn)練剖面隱馬爾可夫模型的參數(shù),進(jìn)行多序列的比對(duì).從實(shí)驗(yàn)結(jié)果可以看出,以Log-odd score作為目標(biāo)函數(shù),QPSO算法相比BW算法和PSO算法,是一種非常有效的訓(xùn)練HMM模型模型的方法。并且從實(shí)驗(yàn)結(jié)果還可以得出,以SOP作為目標(biāo)函數(shù),QPSO比其他所有的方法而言,能產(chǎn)生較好的比對(duì)結(jié)果。這是因?yàn)镼PSO算法是一種全局收斂算法,并且需要調(diào)整的參數(shù)也比較少。

在計(jì)算時(shí)間上,QPSO消耗的時(shí)間和PSO消耗的時(shí)間相當(dāng),遠(yuǎn)遠(yuǎn)大于BW算法的運(yùn)行時(shí)間。平均來(lái)說(shuō),QPSO和PSO的運(yùn)行時(shí)間為6小時(shí),BW的運(yùn)行時(shí)間為5小時(shí).并且,隨著序列個(gè)數(shù)和序列長(zhǎng)度的增加,QPSO消耗的時(shí)間也隨著增大。

在進(jìn)一步的研究工作中,我們將進(jìn)一步改善QPSO算法,來(lái)提高HMM的性能,減少序列比對(duì)的時(shí)間。

[1]GUAN Wei-hong,X Z-Y,ZHU Ping.Nonlinear prediction analysis of properties in protein sequences(Ⅰ)[J].Journal of Food Science and Biotechnology,2008,27(1):71-75.

[2]GUAN Wei-hong,X.Z.-Y.,ZHU Ping.Nonlinear prediction analysis of properties in protein sequences(Ⅱ)[J].Journal of Food Science and Biotechnology,2008,27(2):103-105.

[3]Thompson J D,Higgins D G,Gibson T J.CLUSTAL W:improving the sensitivity of progressive multiple sequence alignment through sequence weighting,position-specific gap penalties and weight matrix choice[J].Nucl Acids Res,1994, 22(22):4673-4680.

[4]Feng D F,Doolittle R.Progressive sequence alignment as a prerequisitetto correct phylogenetic trees[J].Journal of Molecular Evolution,1987,25(4):351-360.

[5]Kim J,Pramanik S,Chung M J.Multiple sequence alignment using simulated annealing[J].Comput Appl Biosci,1994, 10(4):419-426.

[6]E Jung Lee,S F S,Chen-Chia Chuang,et al.Genetic algorithm with ant colony optimization(GA-ACO)for multiple sequence alignment[J].Applied Soft Computing,2008,8(1):15-17.

[7]Mamitsuka H.Finding the biologically optimal alignment of multiple sequences[J].Artif Intell Med,2005,35(1-2):9-18.

[8]Loytynoja A,Milinkovitch M C.A hidden Markov model for progressive multiple alignment[J].Bioinformatics,2003, 19(12):1505-1513.

[9]Jun S,Wenbo X,Bin F.In A global search strategy of quantum-behaved particle swarm optimization[J].Cybernetics and Intelligent Systems,2004,321:325-331.

[10]Jun S,Bin F,Wenbo X.In Particle swarm optimization with particles having quantum behavior[J].Evolutionary Computation,2004,321:325-331.

[11]Kennedy J,Eberhart R.In Particle swarm optimization[J].Neural Networks,1995,1994:1942-1948.

[12]Clerc M,Kennedy J.The particle swarm-explosion,stability,and convergence in a multidimensional complex space[J]. Evolutionary Computation,2002,6(1):58-73.

[13]Thompson J D,Plewniak F,Poch O.A comprehensive comparison of multiple sequence alignment programs[J].Nucleic Acids Res,1999,27(13):2682-2690.

[14]Sonnhammer E L,Eddy S R,Durbin R.Pfam:a comprehensive database of protein domain families based on seed alignments[J].Proteins,1997,28(3):405-420.

(責(zé)任編輯:李春麗)

Multiple Sequence Alignment Based On the Profile Hidden Markov Model

LI Cheng-yuan1, LONG Hai-xia2, SUN Jun1, XU Wen-bo＊1
(1.School of Information Technology,Jiangnan University,Wuxi 214122,China;2.School of Education,Jiangnan University,Wuxi 214122,China)

Multiple sequence alignment(MSA),known as NP-complete problem,is one of the basic problems in computational biology.At present Profile Hidden Markov Model(HMM)was widely used in multiple sequence alignment.This manuscript presented the quantum-behaved particle swarm optimization(QPSO)which was based on particle swarm optimization.The proposed algorithm was used to optimize the profileHMM.Furthermore,an integration algorithm based on the profile HMM and QPSO for the MSA was constructed.Then the approach was evaluated by a set of standard instances which are chosen from nucleotides sequences and the benchmark alignment database,name as BAliBASE.Finally our results are compared with other algorithms.The result shown that the proposed algorithm not only finds out the perfect profile HMM,but also obtains the optimal alignment of multiple sequence.

multiple sequence alignment,profile hidden markov model,quantum-behaved particle swarm optimization

Q 811.4

:A

1673-1689(2010)04-0634-07

2009-08-07

李成淵(1980-),男,江蘇無(wú)錫人,生物信息學(xué)博士研究生。Email:lichengyuaning@gmail.com

＊通信作者:須文波(1946-),男,江蘇無(wú)錫人,教授,博士生導(dǎo)師,主要從事生物信息學(xué)方面的研究。Email:xwb_sytu@hotmail.com