基于POA有向無環圖及隱馬爾科夫模型的優化PHMM算法

王靖會 孫啟明

摘 要：隨著基因組計劃的實施，大量的基因組數據被測序，如何從海量數據中提取出有用的生物學信息，成為生物信息學研究中的重點。序列比對是生物信息學研究的基本方法和主要手段，目前已經有很多比對算法，但如何提高算法敏感性及準確性，仍舊是一個難題。基于這種情況，本文基于POA有向無環圖算法及隱馬爾科夫模型，創新性的提出一種新的算法PHMM算法，這種算法打破傳統的線性序列比對，采用了一種基于圖論的漸進比對算法，對參數選擇具有一定先進性。

關鍵詞：生物信息學；序列比對；POA有向無環圖；隱馬爾科夫模型；PHMM算法

中圖分類號：Q811.2 文獻標識碼：A DOI：10.11974/nyyjs.20160532013

1 引言

隨著人類基因組計劃的實施和發展，基因組和蛋白質組研究的也不斷深化，伴隨而來產生大量的基因序列數據。生物信息學[1]在生物數據分析和處理中的應用越來越廣泛，這使得序列分析[2]成為計算機在生物學中研究的熱點，并使現代生物學研究方法發生了深刻的改變。生物信息學的研究重點主要分2部分：基因組學方面、蛋白質組學方面。簡單的說就是通過分析核酸和蛋白質的序列，得到其表達結構，找到功能及進化的信息。從而闡明生物所攜帶的生物信息及生命意義。

2 多序列比對算法優化

2.1 多序列比對（MSA）問題描述

對于給定的N條序列，S1，S2，…SN，Sk為N條序列中的1條，Sk∈[S1，SN]， Sk=sk1sk2…skn ， k=1，2，…N，n>0.skn為N條序列中的任意1條上的1個殘基，若該序列為DNA分子序列，則skn屬于集合{A，T，C，G}，若該序列為RNA分子序列，則skn屬于集合{A，U，G，C}，若該序列為蛋白質序列，則skn為20種氨基酸中的1種，它的集合中包含20種氨基酸。在生物長期進化發展過程中，會出現少數殘基缺失，插入或變異的現象，因此在多序列比對過程中可能會出現序列中斷（break）現象，既間隙（gap）問題。當出現間隙時，可以用“—”來表示。當用Σ‘表示序列的集合時，Σ=Σ∪{—}。則多序列比對可以按如下定義：

1個多重序列比對A是1個二維字符矩陣，即A={skn‘}，k=1，2，…N，n>0。

其中skn=skn或“—”，并滿足以下條件：

序列的數目等于矩陣的行數；

如果移去每行中的“—”將得到原來的序列；

每1列中不允許所有元素同時為“—”。

從定義可以看出，多序列的比對[3-5]可以模擬生物序列的進化過程，當序列中的殘基經過插入，刪除或替代，則產生了變異，再經過多代的遺傳與積累，最終產生較大差異。

2.2 POA算法

POA算法是1種使用偏序圖（partially ordered graph）的算法，這種算法不同于其他多序列比對算法，摒棄了傳統的線性比對，在POA比對算法中，每個字母代表1個節點，在相鄰的節點間用有向邊表示，1條序列就是1個有向無環的節點序列。對于MSA格式的序列，首先將每條序列都用POA算法畫圖，然后將相同的字母對齊合并，對于不能合并的節點，稱為分離節點，分離節點單獨畫出；對于有多條入邊或多條出邊的節點，稱為匯合節點；與匯合節點相連的多條邊叫做分支。

POA算法計算相似矩陣時首先是用聚合最近鄰聚類算法，根據計算得分構建進化樹。對于匯合節點，它可以有多個前驅，所以POA算法擴展了動態規劃算法，對于2條序列S和S，設其得分為W（u ， v），u和v分別是S和S的節點，為了計算W（u ， v），需要考慮所有可能到達這2個節點的分支，則有：

其中，p是序列S中u節點的前分支，q是序列S中v節點的前分支。Δ是空格罰分，W（u ， v）是u和v的得分。

2.3 隱馬爾科夫模型

隱馬爾科夫模型（Hidden Markov Model ，HMM）[6-7]是由馬爾科夫模型發展而來。馬爾科夫模型是馬爾科夫鏈的模型化過程，這是一種狀態空間中從一種狀態到另一種狀態的隨機過程。如果一個將來狀態僅依賴于現在狀態，而不依賴于過去的狀態，這個過程具有馬爾科夫性，可稱為馬爾科夫過程。用公式表示，為：

隱馬爾科夫模型是由2個隨機變量序列組成的，一條是觀察不到的隱馬爾科夫鏈，（Yn ，n≥0）；另外一條是可以觀察到的隨機序列（Xn ，n≥0），{ Yn≥0}稱為馬爾科夫鏈，{ Xn ≥0}稱為其觀察鏈。馬爾科夫鏈與觀察鏈之間通過一組概率相聯系。

2.4 基于POA算法的HMM模型優化算法PHMM

2.4.1 POA模型的引入

PHMM算法是一種建立于圖論比對與隱馬爾科夫模型結合的一種比對算法，基于圖論比對算法打破傳統的線性比對算法，根據圖中匯合節點計算序列間的距離，然后構建指導樹，按指導樹順序進行漸進比對。

對待比對序列進行POA畫圖，設A為待比對序列的集合，共有N條序列，則A={A1，A2，，…AN}，每條序列長為Li，即含每條序列含有i個殘基，每個殘疾對應自己的位置i，則有L1，L2，…LI。對于集合A中的序列，要找到調和序列。調和序列是保守區域最多的序列，用調和序列與剩余的其他序列進行相似性計算。對于輸入的多個序列，會存在很多重復片段，所以構造POA圖的第1步是去環，采用了深度優先搜索算法找到環，然后根據其他信息去掉序列中的環。去環的方法最簡單的是切斷環中的一條邊，根據環上的信息，選取關聯信息最少的邊進行切斷。經過去環過程后，可以得到一個有向無環序列。將待比對序列全部去環處理后，下一步就是找到調和序列。可以通過對節點的重疊次數進行加權計算，這樣就可以找到最大權值路徑，即要找到的調和序列。

2.4.2 PHMM模型的引入

對于已找到的調和序列AK，作為HMM模型中的初始狀態概率π，待比對序列作為固定的狀態序列，每條序列與調和序列的相似性作為觀察序列的概率，則有：

待比對序列中的任意一條序列Ai與調和序列比對產生的相似概率Q：

accuracu（Ai，Ak ）=（2條序列Ai 與AK 中的正確匹配數目）/（AK 序列數目）

在給定的待比對序列（狀態序列）A={A1，A2，，…AN}下，產生的2條序列比對的相似概率（觀察序列）Q的概率為：

構建進化樹有2類方法：使用特征數據（characters）、使用距離數據（distances）。使用數據特征構建進化樹一般是針對有限的不同狀態的特征，而距離數據一般是用于序列間的差異的衡量。一旦確認了相似性，就可以轉換成距離數據。

計算2個序列間的距離D，可以用匹配的殘基數目m，比上序列的總數目，通過計算出2個序列間的距離后，建立距離矩陣，距離矩陣是所有待比對序列的總結。根據距離矩陣構建進化樹。PHMM算法采用的是類似于UPGMA（Unweighted pair group method using arithmetic averages）的方法構建進化樹。PHMM算法對于待比對序列結合A中的序列，每條序列與調和序列都可以進行期望精確度計算E=Ei（Ai，AK），通過對期望精確度的計算，期望精確度越高，則優先加入比對。

3 結論

隨著實驗方法的增多和檢測手段的發展，生物基因序列數據被大量測序，如何解決這些生物數據，從數據中提取有用信息成為生物信息學研究中的重要任務。而多序列比對作為生物信息學研究中的重要手段，如何研究新的比對方法，優化算法也成為生物信息學中的主要內容。

本文針對這一問題，基于POA無環有向圖的拓撲結構和統計學算法中的隱馬爾科夫模型，研究了一種新的算法PHMM算法。PHMM算法屬于漸進比對算法的一種，對于漸進比對算法參數選擇問題[8]上有一定先進性，可以提高比對質量的準確性。

參考文獻

[1]陳潤生.生物信息學[J].生物物理學報，1999（15），1.

[2]唐玉榮.生物信息學中的序列比對算法[J].計算機工程與應用，2003（29）：5-7.

[3]霍紅衛.序列比較問題的分治法[J].西安電子科技大學學報，1998，25（3）：345-348.

[4]郭衛斌，施保昌，王能超.多重生物序列對準及其算法綜述[J].高科技通訊，2001（6）：96-102.

[5]Yuzhen Ye and Adam Godzik . Multiple flexible structure alignment using partial order graphs[J].Bioinformatics，2005，21（10）：2362-2369.

[6] Ocak，H.，Loparo，K.A. new bearing fault detection and diagnosis based on hidden Markov modeling of vibration signals[J].Acoustics，Speech，and Signal Processing.2001.Proceedings.2001 IEEE International Conference on，2001（5）：3141-3144.

[7]方紹武，戴蓓倩.一種離散隱Markov模型參數的全局優化算法[J].電路與系統學報，2000，30（6）：659-665.

[8]高雨清，陳永彬. 隱Markov模型參數估計的一種新方法[J].自動化學報.1991，17（1）：56-62.