基于條件隨機場進行蛋白質二級結構預測
2009-01-01 00:00:00羅亮邵澤輝
計算機應用研究 2009年3期
（華中科技大學 控制科學與工程系， 武漢 430074）
摘 要：使用一種新的概率圖模型——條件隨機場對蛋白質二級結構進行預測，并給出了模型的構建、訓練以及解碼的算法。應用這一模型對一個典型的蛋白質數據集CB513的二級結構進行了預測，并將預測結果與其他方法進行比較，預測準確度有明顯的提高。
關鍵詞：二級結構預測;條件隨機場;概率圖模型
中圖分類號：O224；TP301 文獻標志碼：A
文章編號：10013695(2009)03083202
Protein secondary structure predict based on conditional random fields
LUO Liang， SHAO Zehui
(Dept. of Control Science Engineering， Huazhong University of Science Technology， Wuhan 430074， China)
Abstract:Conditional random fields(CRFs) is a new probabilistic graph model and this work introduce conditional random fields in protein secondary structure prediction. This paper gave the method of constructing the model and the algorithms to train and decode the model and used the model to predict the second structure of a famous protein dataset(CB513). Finallycompared the results with some other methods.
Key words：prediction of protein secondary structure; conditional random fields(CRFs); probabilistic graph model
蛋白質是生命科學研究的重要對象，研究蛋白質的功能需要了解它們的結構，特別是空間結構，因為結構決定功能。X射線晶體衍射分析法和多維核磁共振技術是測定蛋白質空間結構的兩種主要實驗方法，然而實驗方法不僅耗資耗時，還受實驗條件的限制，而且用實驗方法測定結構的速度和人類的測序速度之間還存在很大的差距。而蛋白質的空間折疊結構取決于構成該蛋白質的氨基酸序列[1]， 即蛋白質的空間折疊結構的全部信息都隱藏在氨基酸序列中。因此蛋白質結構預測問題是當前國際上最具有理論研究和實際應用價值的問題之一。目前，解決蛋白質二級結構預測方法非常多，也取得了比較好的預測結果，如神經網絡[2， 3]、支持向量機[4~8]、隱馬爾可夫模型（HMM）[9~12]、半馬爾可夫模型（SSMM）[13]、最大熵模型（MEM）[14]等。
隱馬爾可夫模型一個最大的缺點就是由于其輸出獨立性假設，導致其不能考慮上下文的特征，限制了特征的選擇。而最大熵模型則解決了這一問題，可以任意地選擇特征，但由于其在每一節點都要進行歸一化，只能找到局部的最優值，同時也帶來了標記偏見的問題（label bias），即凡是訓練集中未出現的情況全都忽略掉。條件隨機場(CRFs)則能很好地解決了這一問題。條件隨機場是一種新的概率圖模型[15]，它并不在每一個節點進行歸一化，而是所有特征進行全局歸一化，具有表達元素長距離依賴性和交疊性特征的能力， 能方便地在模型中包含領域知識，因此可以求得全局的最優值。Liu等人[14]應用條件隨機場對蛋白質折疊進行預測，取得了比較好的預測結果。本文將條件隨機場用于三狀態二級結構預測，也提高了預測的精度。
1 材料和方法
1.1 線性條件隨機場
條件隨機場是一種以給定的輸入節點值為條件來預測輸出節點值概率的無向圖模型。條件隨機場是一個典型的判別式模型，其聯合概率可以寫成若干勢函數聯乘的形式，其中每個最大團都是兩個相鄰變量以一條邊連接構成的子圖。圖1是一種最簡單且最重要的CRFs，稱為線鏈CRFs。線鏈CRFs假設在各個輸出節點之間存在一階馬爾可夫獨立性，其輸出節點被無向邊連接成一條線性鏈。
在圖1中，S1，S2，…，Sn為CRF的N個狀態，O1，O2，…，On為CRFs的N個觀測字符。觀測字符從字符集{V1，…，VM}中取出，對于蛋白質結構預測問題，M=20，即字符集為包含了20個氨基酸殘基符號的字符集。
CRFs模型由N個狀態組成，這樣，在一個輸入序列給定的情況下，對于參數為λ={1，2，…，K}的線性鏈CRFs， 其狀態序列的條件概率由式（1）得到：
P（S｜Q)= exp [∑Tt=1∑Kk=1λkfk(st-1，st，o，t)]/Z0(1)
其中Z0為歸一化因子，由式(2)得到：
Z0=∑s exp [∑Tt=1∑Kk=1λkfk(st-1，st，o，t)](2)
λifi(S，O)在線性鏈條件隨機場模型中包括定義在點上的特征和定義在邊上的特征，一般都取0、1特征。
1.2 基于序列模體的條件隨機場模型
使用CRFs模型前，先將待預測以及用于訓練的蛋白質序列轉換為表示進化信息的序列模體。本文使用序列對比程序PSIBLAST產生序列模體，程序的下載地址為ftp://ftp.ncbi.nih.gov/blast。給定一個長度為L的蛋白質序列，使用PSIBLAST程序與NCBI數據庫中的蛋白質序列進行對比，可以得到一個20×L的位置得分矩陣S（PSSM）：
S=(S1(1)，S1(2)，…，S1(M))(S2(1)，S2(2)，…，S2(M)) （SL（1），SL(2)，…，SL(M))(3)
每個向量st的分量為正數，并且其和為常數S，其值獨立于向量的位置t。
按照DSSP[16]的方法，蛋白質的二級結構分為八種（H，G，I，E，B，T，S，_），本文按照Rost等人[17]的分類方法把八種二階結構簡化為三種（H，E，C）。其中，H、G歸入H；E、B歸入E；剩下的歸入C。按照這種方法，長度為L的蛋白質序列的二級結構可以用列向量O=(o1，o2，…，oL)T表示。
1.3 RCFs的參數訓練與解碼
與HMM、MEMM類似，在進行二級結構預測之前需要首先對RCFs的參數進行訓練。參數訓練是從訓練數據集D學習每個特征的權重參數，求解向量Λ=λ1，λ2，…，λK的過程，這個過程通過最大對數似然估計實現。
將訓練數據集表示為D={O(i)，S(i)}Ni=1。其中：O(i)={o(i)1，o(i)2，…，oiL}是第i個輸入數據序列，S(i)={s(i)1，s(i)2，…，siL}是相應的輸出數據序列。在訓練數據集D下對數似然為
L=∑Ni=1 log P(S(i)｜O(i))(4)
將式(1)(2)帶入并進行高斯先驗調整，式(4)變為
L=∑Ni=1∑Tt=1∑Kk=1λkfk[s(i)t-1，sit，o(i)，t]-∑Ni=1 log Z0(i)-∑k(k2/2σ2)(5)
其中最后一項是用于進行調整的特征參數的高斯先驗值，σ2表示先驗方差。
Lafferty建立CRFs時給出了一種訓練CRFs的算法，然而這種算法對于參數較多的RCFs來說非常慢。2003年Wallach[18]給出了一種LBFGS算法對CRFs的參數進行訓練。LBFGS算法通過充分利用以前的梯度和修正值來近似曲率值的二階方法，因而使用LBFGS算法進行CRFs訓練只要求提供似然函數的一階導數，很大程度上提高了參數訓練的速度。
對于訓練好的CRFs模型，給定待預測的序列，應用Viterbi算法求出產生此序列的最佳狀態路徑，即解碼問題，其算法復雜度為O(NFLC)。對于本文使用的線性鏈RCFs，C=2。根據路徑中各狀態的類型即可確定序列中每個殘基的預測結構。
S=arg max(S/O)=
arg max ∑nj=1log {exp[∑ki=1λifi(sc1，sc2，O)]}/Z(o)(6)
1.4 CB513數據集
由于膜蛋白從組成和結構上都與球蛋白有很大區別，本文研究的是球形蛋白的二級結構，本文使用的數據從不包含膜蛋白的數據集CB513中進行選取。這個數據集可以從http://www.compbio.dundee.ac.uk下載。在使用PSIBLAST進行序列對比時，刪除了殘基數小于40的序列。
2 結果與討論
為了評價預測的準確性， 使用Q3、CH、CC、CE及SOV五個指標。Q3是測量每一個殘基預測的執行情況(不是1就是0)，當x=y，δ(x，y)=1;當x≠y，δ(x，y)=0。總體上蛋白質結構預測的執行情況可以用Q3的平均值Q3表示。其中x表示正確的結構；y表示預測的結構。
Q3=δ[x，y](7)
Q3=(∑Q3)/N(8)
其中：CH表示H結構中殘基預測的準確度，其定義為
CH=（PHnH-uHoH)/([nH+uH][nH+oH]
[pH+uH][pH+oH])(9)
其中：pH為正確預測的數目，nH為正確負預測的數目，oH是過度預測（假陽性）數目，而uH是不定預測殘基（已丟失）的數目。該系數值越接近1，說明該方法越能成功預測螺旋殘基。類似地，CC、CE顯示對于C和E兩種結構預測的準確度。
基于片斷重疊的測度SOV考慮了β折疊和α螺旋的連續片斷特性，更符合蛋白質二級結構的真實情況，其定義為式(10)，更加詳細的說明可以參考文獻[19]。
SOV=[(1/N) ∑1∈{H，E，C}∑S(i)[minov(s1+s2)+δ(s1，s2)]/minov(s1，s2)×len(s1)]×100(10)
五個指標從不同方面反映預測的準確度。本文將篩選過的CB513分為八組，每組約為50條鏈，每組分別作為測試集，剩下的七組作為訓練樣本集進行預測。預測結果如表1所示。
本文進一步將Q3和SOV的平均值與其他文獻中使用支持向量機[6]和最大熵模型[8]得到的結果進行比較，比較結果如表2所示。從表中可以看到，兩個指標均得到了提高，特別是Q3提高了約1%。
表2 與其他預測方法的比較
方法Q3SOV
支持向量機76.673.5
最大熵模型76.976.1
條件隨機場77.773.6
3 結束語
本文成功地將條件隨機場應用于蛋白質三狀態二級結構預測，并取得了比較好的預測結果，類似的方法可以擴展到對蛋白質八狀態的二級結構預測。本文發現，對于E的預測準確度低于對于C和H的預測，并且SOV的提高并不明顯，需要進一步研究改進蛋白質結構圖模型的建立，提高預測準確度。對于RCFs的訓練的算法速度有待進一步的提高。
