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

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

基于剖面隱馬氏模型的多序列比對

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

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

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

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

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

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

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

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

1 剖面隱馬爾可夫模型的拓撲結(jié)構(gòu)

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

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

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

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

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

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

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

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

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

在QPSO算法中,粒子按照下面3個公式進行更新:

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

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

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

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

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

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

由圖3生成的狀態(tài)序列,我們可以進行空位字符‘—’的插入操作,見圖3。所有的序列中與匹配狀態(tài)M相對應(yīng)的氨基酸堿基是比對的,這些氨基酸堿基位于同一列;與插入狀態(tài)Ik相對應(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 評估比對序列的質(zhì)量在利用Viterrbi算法獲得的路徑進行比對后,需要通過基于SP(sum of pairs)打分系統(tǒng)的目標(biāo)函數(shù)[13]對比對結(jié)果進行評估。我們使用下面標(biāo)準的sum-of-pairs打分函數(shù):

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

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

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

3 實驗結(jié)果

3.1 實驗數(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ù)庫中的3個蛋白質(zhì)家族) 3個蛋白質(zhì)家族分別為G5,CagY_M and Interferon,它們來自Pfam數(shù)據(jù)庫[14],見表2。為了避免隱馬爾可夫模型模型中的過擬合,分別把3個蛋白質(zhì)家族分成訓(xùn)練集和測試集。

表2 蛋白質(zhì)家族N:序列的個數(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è)置

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

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

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

1)對于核酸序列,使用ClustalW 1.81的“swgap”置換分數(shù)表,空位開放和空位擴充分別為15、7。

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

3.3 實驗結(jié)果

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

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

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

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

表4 核酸序列比對結(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ì)序列測試集的Log_odd scores平均值和方差Tab.7 HMM log-odds scores±standard error for the validation sets of three protein families

表8 蛋白質(zhì)序列測試集的最優(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ù),算法分別運行20次得到的核酸序列和蛋白質(zhì)序列的平均值的收斂過程。

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

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

圖7 Interferon蛋白質(zhì)序列的Log_odd平均分數(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算法在迭代進行到一半的時候,Logodd score的值就已經(jīng)超過PSO算法和BW算法,并且在整個運行過程中,QPSO都在不斷的提高自身的性能,而PSO算法和BW算法已經(jīng)早早地收斂了。

圖8 Low-short核酸序列的SOP平均分數(shù)Fig.8 Mean SOP scores for Low-short Nucleotide

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

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

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

4 結(jié) 語

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

在計算時間上,QPSO消耗的時間和PSO消耗的時間相當(dāng),遠遠大于BW算法的運行時間。平均來說,QPSO和PSO的運行時間為6小時,BW的運行時間為5小時.并且,隨著序列個數(shù)和序列長度的增加,QPSO消耗的時間也隨著增大。

在進一步的研究工作中,我們將進一步改善QPSO算法,來提高HMM的性能,減少序列比對的時間。

[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-),男,江蘇無錫人,生物信息學(xué)博士研究生。Email:lichengyuaning@gmail.com

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