基于隱馬爾可夫模型對原核生物編碼序列的識別*

山西醫科大學衛生統計教研室（030001） 曹紅艷 馬 靖 李 治 張巖波

基于隱馬爾可夫模型對原核生物編碼序列的識別*

山西醫科大學衛生統計教研室（030001） 曹紅艷 馬 靖 李 治 張巖波△

目的探討隱馬爾可夫模型在大腸桿菌編碼序列識別中的應用，為生物信息挖掘、致病位點研究提供方法參考。方法對大腸桿菌訓練集數據進行訓練建模，并對測試序列進行識別，用特異度、靈敏度以及精確度三個指標進行評價。結果利用本試驗的方法識別編碼序列的靈敏度為73.33%，特異度為67.78%，精確度為70.56%。結論隱馬爾可夫模型能很好地模擬離散狀態間的轉換，適用于識別有狀態轉移、線性序列的數據。

隱馬爾可夫模型 編碼區序列識別 大腸桿菌

隨著2003年人類基因組測序的完成，快速、準確的基因注釋對進一步識別基因，解釋生命的起源和進化等具有重要的意義［1］。基因注釋包括識別出基因序列中的啟動子、編碼區、調控區等區域以及其他一些未被發現的功能片段，其關鍵問題是找出基因組中所有的基因，即基因識別過程［2］。其中，編碼區域識別，即識別DNA序列中的編碼蛋白質部分的序列，顯得尤為重要。基于實驗研究的基因識別耗時長、速度慢，遠滯后于基因測序的速度，因此，需要尋求機器學習的方法用于基因識別［1］。

隱馬爾可夫模型（hidden markov model，HMM）是當前機器學習的研究熱點［3］，由于HMM結構和基因結構非常吻合，且HMM基于堅實的數學理論體系，能保證精確的分析［4］，另外，HMM的計算處理時間遠小于偽支持向量機，對于海量的生物信息數據分析有很大的優勢，因此，HMM廣泛應用于生物信息領域［3］。本文采用HMM對原核生物編碼序列進行識別，以期為生物信息挖掘、致病位點研究提供方法參考。

原理與方法

1.隱馬爾可夫模型定義

隱馬爾可夫模型是由馬爾可夫鏈發展而來的一種隨機模型，由兩個隨機過程｛xk，yk，k＝1，2，3，…｝組成，其中｛xk｝是由離散隱狀態組成的狀態序列，稱為路徑或狀態鏈，描述狀態之間的轉換；｛yk｝是由可觀察字符組成的觀察序列，稱之為觀測鏈，用來描述狀態與觀察值的對應關系。對編碼序列進行識別時，序列的狀態集合為｛編碼區Y，非編碼區N｝，而序列是由四個堿基組成，故觀察集合為｛A，T，C，G｝。作為觀察者，不能直接看到狀態，而是通過一個隨機過程去感知狀態的存在及其特性，因而稱為“隱”馬爾可夫模型［5］。圖1表示了DNA序列識別的隱馬爾可夫模型內部關系。

圖1 序列識別的HMM內部示意圖

2.隱馬爾可夫模型的三個基本問題

評估問題：根據給定模型求某個觀察值序列發生的概率P（0｜λ），用來評估模型和給定觀察輸出序列的匹配程度，從而達到在一系列候選對象中選取最佳的匹配對象。常用算法為向前／向后算法，該算法可以有效地減少計算量［6－7］。

解碼問題：在給定觀察序列和模型的前提下，求產生觀察值序列的最有可能的狀態序列，常用算法為Viterbi算法。

學習問題：根據給定的觀察值序列，調整模型參數，使其產生觀察值序列的概率P（0｜λ）最大。常用算法為Baum-Welch算法。

實例分析

數據來源于美國國家生物信息技術中心（NCBI），從共享資源中下載到已標識出編碼區和非編碼區的大腸桿菌全基因組序列。

1.隱馬爾可夫模型的建立

采用Baum-Welch算法，利用Rstudio軟件中的HMM包，對大腸桿菌編碼區和非編碼區序列分別進行訓練建模：選取三分之二的大腸桿菌編碼區序列共2750條建立HMM-gene模型，同樣針對2330條大腸桿菌非編碼序列建立HMM-nogene模型。迭代次數設置為100次，delta默認為1E-9，偽數默認為0。序列的狀態集合為｛編碼區Y，非編碼區N｝，觀察集合為｛A，T，C，G｝。初始的轉移概率矩陣A和發射概率B的定義如表1，訓練后建立的模型結果如表2。

表1 隱馬爾可夫模型初始參數

表2 隱馬爾可夫模型參數

2.隱馬爾可夫模型對DNA序列的識別結果

根據建立好的HMM-gene模型以及HMM-nogene模型，可以直接判斷某位置上核苷酸的屬性，但由于基因刪除、插入等多種基因突變的出現，單一分析某一位點的屬性顯然是不夠的，需要對一定長度序列的性質進行判斷。分別統計在兩模型下，被識別為編碼狀態的核苷酸與識別為非編碼狀態的核苷酸，即Y／N的比例，分別記為（1）和（2）。以兩者的差作為特征指標，若差值大于0，則認為特征傾向于HMM-gene模型，判斷序列為編碼序列，若差值小于0，則識別為非編碼序列，差值為0時，尚不能判斷。該方法避免了復雜的計算，同時也避免了下溢現象，即由于HMM在計算概率時利用了條件概率思想，因此多次計算結果會出現無限趨于0的現象。

從余下的1／3序列中隨機選取編碼和非編碼序列各180條作為測試數據集。表3，表4分別為利用HMM-nogene及HMM-gene模型對編碼和非編碼序列各180條的識別結果，當全部識別為Y時，為了避免分母為0，將N記為1。表5為編碼和非編碼序列識別結果比較，認為兩者識別結果差別有統計學意義（χ2＝10.840，P＝0.028），其中，非編碼序列中相持現象較多，可能與其結構有一定的相關性。

表3 編碼序列識別結果

表4 非編碼序列識別結果

表5 編碼序列和非編碼序列識別結果比較

3.識別結果的評價

識別結果的評價采用敏感度Sn（sensitivity）、特異度Tn（specificity）及精確度Ac（accuracy）指標［8］。靈敏度，又稱之為真陽性率，反映了正確識別出編碼序列的能力。特異度，又稱真陰性率，反映了正確識別出非編碼序列的能力。精確度定義為特異度和靈敏度的均數。三者范圍均在0～1之間，值越大，評價效果越理想。公式如下：

其中TP表示編碼序列中被正確識別為編碼序列的數目，FN表示編碼序列中未被正確識別為編碼序列數目，TN表示非編碼序列中被正確識別為非編碼序列，FP表示非編碼訓練中未被正確識別出來的非編碼序列。

表6 對HMM模型識別序列的評價

討 論

HMM不限制輸入序列的長度，不需要指導者，有形成模塊或層次結構進而形成完整的模型識別系統的優點，非常適用于有狀態轉移、線性序列的數據，同時，HMM運算速度比支持向量機快，對于大樣本的生物信息數據分析有很大的優勢［3］，故在基因識別、序列比對等生物信息領域應用廣泛。

但是，本文中HMM對原核生物編碼序列的識別精確度不夠理想，小于同類研究［3］，究其原因如下：（1）對訓練序列作預處理時，僅剔除了序列長度小于80bp或大于20000bp的序列，未采取其他措施控制訓練序列的質量，進一步可選取等長度、較短序列進行訓練，以研究訓練序列長度對識別結果的影響；（2）未考慮原核生物的重疊信息，進一步應對兩條基因的重疊序列進行處理，以增大識別精度；（3）一階隱馬爾可夫模型有三個條件獨立基本假設，即狀態獨立性、觀察獨立性、狀態與具體時間無關，事實上，在此刻出現的觀測輸出概率不僅依賴于系統前一狀態，也很可能同之前的系統狀態有關，因此進一步可在HMM的基礎上進行堿基相關性的計算［9］。

總之，本文通過對原核生物序列進行訓練后建立HMM，對DNA序列進行了識別，用特異度、靈敏度以及精確度三個指標進行了評價，為生物序列識別提供了方法學參考。今后可進一步研究HMM在序列比對中的應用，以及研究HMM對真核生物序列以及蛋白質序列的識別。

1.Goel N，Singh S，Aseri TC.A Review of Soft Computing Techniques for Gene Prediction.ISRN Genomics，2013，2013.

2.房穎.基于統計的基因識別算法研究［碩士論文］.長春：吉林大學，2007.

3.羅澤舉，李艷會，宋麗紅，等.基于隱馬爾可夫模型的DNA序列識別.華南理工大學報（自然科學版）.2007，135（8）：123-126.

4.Stormo G.Gene-finding approaches for eukaryotes.Genome research，2000，10（4）：394-397.

5.潘海燕，丁元林，胡利人，等.隱Markov模型及其在慢性病流行病學研究中的應用.中國衛生統計，2009，26（1）：38-40.

6.Baum LE，Egon JA.An inequality with applications to statistical estimation for probabilistic functions of a Markov process and to amodel for ecology.Bull Amer，Meterol，SOC，1967，73：360-363.

7.Baum LE，Sell GR.Grow th functions for transformations on manifolds. Pac，J，Math，1968，27（2）：211-227.

8.Tompa M，Li N，Bailey TL，et al.Assessing computational tools for the discovery of transcription factor binding sites.Nat Biotechnol，2005，1（23）：137-144.

9.張爽.基于HMM的轉錄因子結合位點識別方法研究［碩士學位論文］.沈陽：東北師范大學，2009.

（責任編輯：劉壯）

The Coding Sequence Recognition of Prokaryote Based on Hidden Markov Model

Cao Hongyan，Ma Jing，Li Zhi，et al.（School of Public Health，Shanxi Medical University（030001），Taiyuan）

ObjectiveTo explore the identification of Escherichia coli coding sequence with Hidden Markov Model，so as to provide methods for the research of mining biological information and pathogenic loci.MethodsWe train the data set of Escherichia coli to model and identify the test set，and then evaluate the results using specificity，sensitivity and accuracy.ResultsThe specificity is 67.78%，the sensitivity is 73.33%and the accuracy is 70.56%based on the method of the paper.ConclusionHidden Markov Model can simulate the transformation of the discrete state very well，applied to identify the data of transformation state and linear sequence.

Hidden markov model；Coding region recognition；Escherichia coli

*：國家自然科學基金資助項目（31071156）

△通信作者：張巖波，Email：yanbozh＠126.com