山西醫(yī)科大學(xué)衛(wèi)生統(tǒng)計教研室（030001） 曹紅艷 馬 靖 李 治 張巖波

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

山西醫(yī)科大學(xué)衛(wèi)生統(tǒng)計教研室（030001） 曹紅艷 馬 靖 李 治 張巖波△

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

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

隨著2003年人類基因組測序的完成，快速、準(zhǔn)確的基因注釋對進(jìn)一步識別基因，解釋生命的起源和進(jìn)化等具有重要的意義［1］。基因注釋包括識別出基因序列中的啟動子、編碼區(qū)、調(diào)控區(qū)等區(qū)域以及其他一些未被發(fā)現(xiàn)的功能片段，其關(guān)鍵問題是找出基因組中所有的基因，即基因識別過程［2］。其中，編碼區(qū)域識別，即識別DNA序列中的編碼蛋白質(zhì)部分的序列，顯得尤為重要?；趯嶒炑芯康幕蜃R別耗時長、速度慢，遠(yuǎn)滯后于基因測序的速度，因此，需要尋求機器學(xué)習(xí)的方法用于基因識別［1］。

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

原理與方法

1.隱馬爾可夫模型定義

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

圖1 序列識別的HMM內(nèi)部示意圖

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

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

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

學(xué)習(xí)問題：根據(jù)給定的觀察值序列，調(diào)整模型參數(shù)，使其產(chǎn)生觀察值序列的概率P（0｜λ）最大。常用算法為Baum-Welch算法。

實例分析

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

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

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

表1 隱馬爾可夫模型初始參數(shù)

表2 隱馬爾可夫模型參數(shù)

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

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

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

表3 編碼序列識別結(jié)果

表4 非編碼序列識別結(jié)果

表5 編碼序列和非編碼序列識別結(jié)果比較

3.識別結(jié)果的評價

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

其中TP表示編碼序列中被正確識別為編碼序列的數(shù)目，F(xiàn)N表示編碼序列中未被正確識別為編碼序列數(shù)目，TN表示非編碼序列中被正確識別為非編碼序列，F(xiàn)P表示非編碼訓(xùn)練中未被正確識別出來的非編碼序列。

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

討 論

HMM不限制輸入序列的長度，不需要指導(dǎo)者，有形成模塊或?qū)哟谓Y(jié)構(gòu)進(jìn)而形成完整的模型識別系統(tǒng)的優(yōu)點，非常適用于有狀態(tài)轉(zhuǎn)移、線性序列的數(shù)據(jù)，同時，HMM運算速度比支持向量機快，對于大樣本的生物信息數(shù)據(jù)分析有很大的優(yōu)勢［3］，故在基因識別、序列比對等生物信息領(lǐng)域應(yīng)用廣泛。

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

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

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

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

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

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

5.潘海燕，丁元林，胡利人，等.隱Markov模型及其在慢性病流行病學(xué)研究中的應(yīng)用.中國衛(wèi)生統(tǒng)計，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的轉(zhuǎn)錄因子結(jié)合位點識別方法研究［碩士學(xué)位論文］.沈陽：東北師范大學(xué)，2009.

（責(zé)任編輯：劉壯）

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

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

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