一個基于隱馬爾可夫模型和生物知識修正.doc

精品論文一個基于隱馬爾可夫模型和生物知識修正的 cpg 島識別系統(tǒng)徐瑜，蘭曼5（華東師范大學(xué)信息科學(xué)技術(shù)學(xué)院， 上海 200241）摘要： cpg 島的存在能識別某些基因的啟動子，而且 cpg 島的異常甲基化多與人類腫瘤的 發(fā)生有關(guān)，因此 cpg 島的識別在生物基因組測序中很重要。本文實現(xiàn)了一個基于隱馬爾可夫模型（hmm）和后期生物知識修正的 cpg 島識別系統(tǒng)，在 embl 的 dna 序列數(shù)據(jù)集上10進行系統(tǒng)性能測試的結(jié)果顯示，該系統(tǒng)對于 cpg 島有較好的識別能力，同時又比較精確地 定位 cpg 島的位置。與其它常用 cpg 島識別工具的對比實驗結(jié)果表明，該系統(tǒng)的識別準(zhǔn)確性不亞于其它軟件。此外，該系統(tǒng)的優(yōu)點是一經(jīng)訓(xùn)練，可以用于自動識別。關(guān)鍵詞：隱馬爾可夫模型；機器學(xué)習(xí)；cpg 島識別；知識修正中圖分類號：請查閱中國圖書館分類法15a cpg island identification system based on hmm anddomain knowledge revisionxu yu, lan man(school of information science and technology, east china normal university, shanghai20200241)abstract: cpg islands are often associated with the promoters of most housekeeping genes and many tissue-specific genes. finding candidate regions for aberrant dna methylation contributesto the understanding of the epigenetic causes of cancer. therefore, identification of cpg islands isvery important in genome mapping projects. in this work, we implemented an automatic cpg25islands identification system, i.e. cpg-discover, based on hidden markov model (hmm) and post-processing revision including biologic domain knowledge. the experimental results on the data set of human dna sequences from emnl showed that this system has good identification performance. moreover, in comparison with other popular tools, cpg-discover shows comparable effectiveness. in addition, the significant advantage is it can be applied to other data sets once built30on an annotated training data set.key words: hidden markov model; machine learning; cpg islands identification; knowledge modification0引言生物學(xué)上，cpg島（cpg islands）是指dna上包含大量相鄰的胞嘧啶（c）、鳥嘌呤（g），35以及使兩者相連的磷酸酯鍵（p）的一個區(qū)域。通常，cpg島的長度為幾百到幾千個堿基對（nucleotides，單位bp）。在人的基因組中，如果雙堿基對cg出現(xiàn)，則c通常被甲基化。并 且，甲基化的c很快會突變成t，因此基因組中cpg島非常少。然而，在基因的起始位置（啟 動子），因為功能的保守性，其序列很少突變，在哺乳類基因中的啟動子上含有約40%的cpg 島，人類基因中含有約70%的cpg島，因此，cpg島的存在能識別某些基因的啟動子，并可40作為限制酶的辨識位置。從已知的dna序列統(tǒng)計來看，幾乎所有的管家基因（house-keeping genes）及約40%的組織特異性基因（tissue-specific genes）的附近(通常是在5上游區(qū)域)都含基金項目：教育部博士點基金(20090076120029)作者簡介：徐瑜（1988-）， 女，碩士研究生，自然語言處理，文本挖掘通信聯(lián)系人：蘭曼（1974-）， 女，副教授，自然語言處理，數(shù)據(jù)挖掘，生物文本挖掘等. e-mail:- 2 -有非甲基化的cpg島，其序列可能包括基因轉(zhuǎn)錄的啟動子的第一外顯子1。因此，cpg島的 識別在生物基因組測序工作中是非常重要的。此外，人類腫瘤的發(fā)生與多種腫瘤相關(guān)基因的cpg島的甲基化水平有關(guān)，cpg島甲基化可以調(diào)控基因轉(zhuǎn)錄的效率，使基因轉(zhuǎn)錄失活，是細(xì)45胞內(nèi)dna轉(zhuǎn)錄狀態(tài)的重要表觀遺傳學(xué)標(biāo)記，因此，啟動子區(qū)域的cpg島的異常甲基化也是識 別癌癥的重要標(biāo)志之一。針對cpg島的定義，1987年，gardiner-garden和formmer2認(rèn)為長度在200bp以上，g+c 含量大于50，并且實際cpg含量與期望cpg含量的比值(obscpgexpc g)大于0.6的區(qū)域 即為cpg島。2002年，daiya takai 和 peter a. jones3的研究認(rèn)為，長度在500bp以上，g+c50含量大于55，并且實際cpg含量與期望cpg含量的比值(obscpgexpc g)大于0.65的區(qū)域 更可能和基因的5區(qū)域有關(guān)。因此，根據(jù)研究的不同的目的，生物學(xué)研究人員可以采用不同 的參數(shù)進行cpg島預(yù)測。dna序列中的cpg島既可以通過生物學(xué)實驗的方法識別，也可以通過計算機來識別可能 的cpg島。由于dna序列的多樣性，如果借助于生物信息和計算機工具首先識別出潛在的55cpg島，就可以大大減少生物實驗的成本和盲目性，極大地提高生物實驗的時間效率和針對 性。近年來，利用計算機開發(fā)cpg島識別工具，吸引了許多生物學(xué)領(lǐng)域和傳統(tǒng)計算機領(lǐng)域研 究人員的興趣，他們已經(jīng)開發(fā)出許多有用的工具來幫助生物學(xué)家提高工作效率。根據(jù)這些工 具所采用的方法，我們將它們大致分為兩類。第一類方法是采用傳統(tǒng)的滑動窗口（sliding window），這些工具包括cpgis3，cpgprod4和cpgie5等?；瑒哟翱诜◤念^至尾對dna60序列進行掃描（每次移動1n位），分別計算窗口內(nèi)序列是cpg島 (用“+”表示) 和不是cpg 島 (用“-”表示) 的概率值（通常采用log-odds ratio），如果高于設(shè)定的域值，則認(rèn)為窗口內(nèi) 的序列是cpg島，否則認(rèn)為不是。采用滑動窗口法的不足之處在于：（1）窗口的大小很難 確定，通常窗口的大小直接影響cpg島的長度和數(shù)量，如果窗口過大，一些長度較短的cpg 島會被合并在一起，從而被預(yù)測成一個大的cpg島，降低識別的準(zhǔn)確率；（2）窗口只能向65一個方向移動，對于之前錯誤的判斷結(jié)果，不能在后面得到發(fā)現(xiàn)和修改，從而影響識別的精 度；（3）系統(tǒng)運行時間較長。第二類方法是采用聚集（clustering）的技術(shù)，這些工具包括 cpgcluster6和cpgif7等。這類工具的基本思想是基于cpg島定義中的密度或物理距離等生 物學(xué)屬性，不斷迭代地掃描序列，并將符合這些屬性的相鄰cpg島合并為大的cpg島簇。聚 集方法的提出是為了解決滑動窗口方法的不足，然而它帶有如下的缺陷：（1）識別的結(jié)果70依賴于掃描序列的順序，而不同的掃描順序則會識別出不同的cpg島，因此系統(tǒng)穩(wěn)定性不足；（2）識別的敏感度較低，傾向于識別長度較長的cpg島序列。 與上述兩類方法不同，在本文，我們嘗試使用隱馬爾可夫模型（hmm）建立一個自動識別cpg島的系統(tǒng)，并結(jié)合生物知識對預(yù)測結(jié)果進行后期修正，從而達到較為理想的識別預(yù) 測結(jié)果。隱馬爾可夫模型是一種機器學(xué)習(xí)的方法，它廣泛應(yīng)用于序列分析和模式識別的研究75中，例如語音識別、手寫識別，命名實體識別等。這類應(yīng)用的一個重要特點是，數(shù)據(jù)的序列 性。這種序列可以被認(rèn)為是后一狀態(tài)依賴于前一狀態(tài)的時間序列，即具有馬爾可夫性。正是 因為考慮到dna序列（堿基對）也可以被看做是具有馬爾可夫性的有序數(shù)據(jù)，我們在本文 中利用hmm模型實現(xiàn)了一個自動從人類dna序列中識別cpg島的系統(tǒng)，即cpg-discover系 統(tǒng)。80本文中，我們首先介紹結(jié)合hmm模型和生物知識修正的cpg-discover系統(tǒng)的框架結(jié)構(gòu)。 然后在歐洲分子生物學(xué)實驗室（embl）發(fā)布的人類基因中cpg島的數(shù)據(jù)集上測試本系統(tǒng)的 性能，并與其它常用的cpg島識別工具（基于滑動窗口和聚集技術(shù)）進行對比實驗，深入分析和討論實驗結(jié)果。最后我們總結(jié)這一工作和未來的方向。1cpg-discover 系統(tǒng)框架85cpg-discover 系統(tǒng)的框架結(jié)構(gòu)主要包括三個子系統(tǒng)：（1）建立 hmm 模型的參數(shù)系統(tǒng)；（2）基于 hmm 模型建立 cpg 島識別系統(tǒng)；（2）基于生物知識的后期識別修正系統(tǒng)。下 面我們分別詳細(xì)介紹這三個子系統(tǒng)的工作。1.1建立 hmm 模型的參數(shù)系統(tǒng)這個子系統(tǒng)完成對 hmm 模型的求解過程，即實現(xiàn)參數(shù)初始化、參數(shù)逼近和參數(shù)投票這90三個步驟，建立起 hmm 模型的參數(shù)系統(tǒng)。在參數(shù)初始化階段，針對 4 種堿基（a，c，t， g）分別在 cpg 島內(nèi)外的狀態(tài)，我們標(biāo)記此 hmm 模型包括 8 個狀態(tài)，分別是 a+, c+, g+, t+, a, c, g,t，標(biāo)號中“+”和“”分別表示此核苷酸在和不在 cpg 島上。這 8 個狀態(tài)之間的 轉(zhuǎn)移方式如圖 1 所示。這 8 個狀態(tài)間的狀態(tài)轉(zhuǎn)移概率可以從有標(biāo)識的訓(xùn)練序列中得到，后面 會詳細(xì)敘述。本文中，我們基于歐洲分子生物學(xué)實驗室（embl）核苷序列數(shù)據(jù)庫進行實驗95建立 hmm 模型，因此，在以下三個參數(shù)求解的過程中，我們以 id 為 bc008880 的人類 dna序列為例來進行說明。c+g+a+t+a-t- c-g- 11 -1001.1.1參數(shù)初始化圖 1hmm 模型中 8 個狀態(tài)之間的轉(zhuǎn)移圖示fig. 1 the transform of 8 states in hmm105110這個步驟完成 hmm 模型中狀態(tài)轉(zhuǎn)移概率矩陣（a）、發(fā)射概率矩陣（b）和狀態(tài)概率 矩陣（）的初始化，作為下一步參數(shù)逼近的輸入種子。三個概率矩陣的初始化步驟如下：（1）初始化狀態(tài)轉(zhuǎn)移概率矩陣 a：對于每個訓(xùn)練樣本序列，我們采用基于長度為 2 的 滑動窗口的方法來初始化參數(shù)，即通過計算從第一個狀態(tài)轉(zhuǎn)移到第二個狀態(tài)的出現(xiàn)頻數(shù)，來 計算轉(zhuǎn)移概率。（2）初始化發(fā)射概率矩陣 b：假設(shè)一個狀態(tài)到一個符號的發(fā)射概率都為 0 或 1。（3）初始化狀態(tài)概率矩陣 ：每個狀態(tài)的初始化概率由序列中各個狀態(tài)的出現(xiàn)頻數(shù)計 算。參數(shù)完成初始化后，我們得到初始化后的各個概率矩陣。表 1 顯示了人類 dna 序列bc008880 經(jīng)過初始化之后的三個概率矩陣的分布結(jié)果。1.1.2參數(shù)逼近在第二步中，我們采用 baum-welch 算法來重新估計基于訓(xùn)練序列的概率矩陣。上一步 的三個初始化概率矩陣作為 baum-welch 算法的輸入，輸出期望的轉(zhuǎn)移和發(fā)射概率矩陣，再115將得到的期望概率代替舊的參數(shù)，這個過程不斷迭代進行直至達到收斂。在每次迭代過程中，模型概率的準(zhǔn)確度都得到提高，直到一個極限概率，整個迭代過程最終收斂于一個局部最優(yōu) 解。參數(shù)估計后，對每個訓(xùn)練序列，我們都得到三個逼近后的概率矩陣。表 2 顯示了人類 dna 序列 bc008880 經(jīng)過 baum-welch 近似之后的三個概率矩陣的分布結(jié)果。120a+c+t+g+a-c-t-g-0.19150.23450.51060.063800000.18520.23460.34570.22220.01230000.12930.32760.33620.206900000.13730.23530.50980.1176000000000.39360.13650.23290.236900000.40830.1750.10.3167000.005400.32970.15680.25950.248600000.19180.16440.29680.347a：a+ c+ t+ g+ a- c- t- g-:表 1 人類 dna 序列 bc008880 經(jīng)過初始化之后的三個概率矩陣的結(jié)果tab. 1 the result of initialization of dna sequence (bc008880)actg10000100001000011000010000100001b： a+ c+ t+ g+ a-c- t- g-0.04388 0.0756 0.0476 0.1092 0.2334 0.1120 0.2045 0.1718125a+c+t+g+a-c-t-g-0.00610.04830.00580.94060.00120.00160.00170.00200.56020.02560.12330.29120.00110.00150.00160.00230.00570.48620.16540.33490.00170.00460.00570.00280.01780.23530.50980.11760.00120.00140.00210.00350.00110.00110.00110.00110.87760.11810.00360.00330.00120.00110.00110.00110.01670.26820.70460.01300.00110.00120.00110.00120.01660.28190.01690.68710.00110.00110.00110.00110.00280.22230.29170.4856a： a+ c+ t+ g+ a-c- t-g-表 2 人類 dna 序列 bc008880 經(jīng)過 baum-welch 逼近之后的三個概率矩陣的結(jié)果tab. 2 the result of initialization of dna sequence (bc008880)b：actga+0.43010.04100.31770.2140c+0.05590.59650.00800.3426t+0.00530.24020.71210.0454g+0.13010.14830.00830.7163a-0.92480.03080.0040.0450c-0.35340.54340.02170.0845t-0.55900.00970.42400.0103g-0.06590.01960.40140.5162：0.001113 0.999155 0.001543 0.001164 0.001001 0.001019 0.001002 0.0010031301.1.3參數(shù)投票至此，對每個 dna 訓(xùn)練序列，我們都得到三個逼近后的概率矩陣。為了整合各個不同 長度的訓(xùn)練序列的概率矩陣，我們采取投票（vote）平均技術(shù)，就是說，對每個概率矩陣， 根據(jù)每個訓(xùn)練序列的長度對其進行歸一化（normalization），每個 dna 訓(xùn)練序列的權(quán)重計算公式為： w = ni，其中，n 表示第 i 個訓(xùn)練序列中核苷酸的數(shù)目，n 表示訓(xùn)練集中的核i n i135苷酸總數(shù)。采用投票平均后，我們得到最終的 hmm 模型的參數(shù)。表 3 顯示了使用 1400 條 “bc”打頭的人類 dna 序列作為 hmm 系統(tǒng)的訓(xùn)練數(shù)據(jù)，經(jīng)過加權(quán)平均（投票）之后的 hmm 模型的三個概率矩陣的分布結(jié)果（關(guān)于數(shù)據(jù)集見后面 2.1 節(jié)）。1401451.2基于 hmm 模型建立 cpg 島識別系統(tǒng)經(jīng)過求解和加權(quán)平均之后得到的 hmm 模型，可以通過 viterbi 解碼算法來尋找能夠產(chǎn) 生給定符號序列的最大似然狀態(tài)序列，即對每個測試 dna 序列中的每個堿基對（atcg）， viterbi 解碼算法能夠輸出每個堿基對對應(yīng)的狀態(tài)為”+”或者”-”（其中“+”標(biāo)識堿基對在 cpg 島內(nèi)，”-”表示堿基對在 cpg 島之外），并使得標(biāo)識后的 dna 序列輸出狀態(tài)序列為最優(yōu)解。表 3 使用 1400 條以“bc”開頭的人類 dna 序列作為訓(xùn)練數(shù)據(jù)集，經(jīng)過加權(quán)平均（投票）之后的 hmm模型的三個概率矩陣分布tab. 3 the result of majority vote, using 1400 dna-sequencea：a+c+t+g+a-c-t-g-a+0.238770.17100.14740.31790.00880.03950.03450.0492c+0.22390.26740.30300.08320.02120.05930.03230.0165t+0.051780.26000.20220.37280.00660.03770.03200.0437g+0.188420.30830.14300.22190.02210.03620.03050.0565a-0.006390.01700.01230.02340.55030.09970.09690.201c-0.066440.06970.04020.01890.17580.24300.30320.0898t-0.015650.03810.03130.03670.05420.21380.23460.3825g-0.057660.03640.03300.05960.17050.26060.15760.2317b：actga+0.80710.06550.06470.0657c+0.06990.83340.04080.0589t+0.10640.07040.75500.0712g+0.04570.06380.05280.8406a-0.87800.03400.04130.0497c-0.08120.81540.04100.0654t-0.13210.06330.73610.0715g-0.05260.05790.05540.8371150155160165170：0.05035 0.1919 0.03460 0.2662 0.04592 0.1488 0.03663 0.23251.3基于生物知識的后期識別修正系統(tǒng)由viterbi 算法求解出的dna序列輸出狀態(tài)是基于訓(xùn)練模型計算而得，系統(tǒng)輸出的狀態(tài) 序列并不一定滿足生物學(xué)上對cpg島的定義。為了進一步提高結(jié)果的精度，我們結(jié)合cpg島 的生物定義，如cpg島中c和g的概率、相鄰cpg片段間距以及cpg片段的長度閾值等，對 hmm模型識別dna序列產(chǎn)生的誤差進行后期的知識修正。我們對cpg島作如下的限定：至 少包含140個堿基對的長度，區(qū)域內(nèi)gc所占含量超過60%，且cpg的觀察值預(yù)測值比例必 須高于0.6（根據(jù)研究的不同目的，可以在軟件系統(tǒng)中自行修改這些限定的值）。這個后期 識別修正系統(tǒng)包含下面三個模塊：（1）相鄰cpg島的合并模塊：由于計算誤差，系統(tǒng)經(jīng)常把cpg島中少量的核苷酸（通 常長度在1bp到15bp之間）錯誤地識別為非cpg島上的核苷酸。為了解決這個問題，我們對 相鄰cpg島設(shè)定最小距離限定。經(jīng)過多次實驗，這個最小距離的閾值設(shè)定為20，即，如果識 別出的兩個相鄰cpg島間的距離小于20，系統(tǒng)就修改這些間隔中的非cpg島核苷酸狀態(tài)，使 之狀態(tài)轉(zhuǎn)變?yōu)閏pg島上核苷酸，即合并這兩個相鄰cpg島為一條長cpg島。合并后的cpg島 需要進一步經(jīng)過后面兩個模塊的修正。（2）密度檢測過濾模塊：生物學(xué)上，cpg島通常含有較高的g/c（或c/g）含量(密度高 于60%)，如果前一步的錯誤合并或者系統(tǒng)識別的誤差，可能產(chǎn)生不滿足c/g含量要求的cpg 島序列。因此，對識別出的不能滿足最低密度要求的cpg島序列，即區(qū)域內(nèi)gc所占含量至 少達到60%，這個密度檢測過濾模塊將改變這些cpg島中所有核苷酸狀態(tài)為非cpg島狀態(tài)。（3）長度檢測過濾模塊：除了（2）中的密度約束模塊外，這個系統(tǒng)還進行了長度約束， 即根據(jù)cpg島的長度定義，去除長度少于140個堿基對的cpg島序列，即，對識別出的長度175180小于140bp的cpg島將被改變?yōu)榉莄pg島狀態(tài)。經(jīng)過以上三個后期修正模塊之后的輸出序列狀態(tài)，才成為cpg-discover系統(tǒng)的最終輸入結(jié)果。表4列出了人類dna序列bc009465在修正前后的對比結(jié)果，其中491255分別表示cpg島序列的起始位置為第49個堿基對和結(jié)束位置為第1255個堿基對，這6個評估指標(biāo)的含義在后面2.2節(jié)中有詳述。從表4可以看到，基于生物知識修正后的識別結(jié)果除了correlation coefficient指標(biāo)沒有變化外，其它的指標(biāo)都有提高，特別是specificity，從原來的33.76%增至了51.68%（提高60.15%）。原來預(yù)測位置在1284和2961509的兩條cpg島因為相隔只有12個堿基對（bp），小于最小間隔閾值20bp，因而被修正合并為一條長的cpg島；而原來預(yù)測位置在15841695的cpg島序列長度只有120bp，在后處理模塊中由于長度過短，不滿足 生物學(xué)上cpg島的長度要求，因而被過濾除去。表 4 人類 dna 序列 bc009465 基于生物知識的修正前后識別結(jié)果的對比tab. 4 the comparison of domain knowledge revision (bc009456)評測指標(biāo)未修正的結(jié)果修正后的結(jié)果真實cpg島的始末端位置491255491255系統(tǒng)識別的cpg島的始末端位置128429615091584169511509accuracy76.80%83.51%sensitivity99.09%100%specificity33.76%51.68%positive predictive value74.29%79.99%performance coefficient73.78%79.99%correlation coefficient47.72%47.72%1851901952002實驗部分為了對 cpg-discover 系統(tǒng)的性能做一個準(zhǔn)確而全面的分析和比較，我們從歐洲分子生 物學(xué)實驗室(embl)的核苷序列數(shù)據(jù)庫（http:/www.ebi.ac.uk/embl/）中下載了關(guān)于人類基因 的 cpg 島數(shù)據(jù)庫 emb173hum（ftp:/ftp.ebi.ac.uk/pub/databases/cpgisle），并在這個數(shù)據(jù)集上 測試了本系統(tǒng)的性能。這個數(shù)據(jù)庫是由歐洲生物信息中心負(fù)責(zé)維護，與日本基因數(shù)據(jù)庫（ddbj）和美國國立生物技術(shù)信息中心（ncbi）相互合作，他們之間每天都要交換最新的 核苷序列數(shù)據(jù)。此外，為了比較 cpg-discover 系統(tǒng)與其它常用的 cpg 島識別工具（基于滑 動窗口和聚集技術(shù)）的性能，我們還在這個數(shù)據(jù)集上選擇部分 dna 序列進行了對比實驗。2.1數(shù)據(jù)來源emb173hum 數(shù)據(jù)庫包含了當(dāng)前 embl 數(shù)據(jù)庫發(fā)布的已標(biāo)定 cpg 島位置的人類 dna 序 列，共有 233,004 條 dna 序列和 142,325 個 cpg 島，其中只有 61,051 條 dna 序列中包含 有一個或者若干個 cpg 島，平均每條 dna 序列含有 2.33 個 cpg 島。我們下載了這個數(shù)據(jù) 庫中以“ab”標(biāo)記開頭的 1,955 條 dna 序列和以“bc”標(biāo)記開頭的 6,124 條 dna 序列，除去 數(shù)據(jù)庫中的空序列和有錯誤信息的序列之后，我們得到近 8,000 條 dna 序列。其中，以“ab” 標(biāo)記開頭的 dna 序列中平均每條 dna 序列含有 1.70 個 cpg 島，以“bc”標(biāo)記開頭的 dna 序列中平均每條 dna 序列含有 1.34 個 cpg 島。為了使實驗結(jié)果真實合理有效，我們從以 “ab”和“bc”標(biāo)記開頭的 dna 序列中分別選擇了 1400 條訓(xùn)練數(shù)據(jù)和 200 條測試數(shù)據(jù)（共有2800 條訓(xùn)練數(shù)據(jù)和 400 條測試數(shù)據(jù)），選擇的原則是這個數(shù)據(jù)子集里面的 cpg 島的分布分 別與原來“ab”和“bc”標(biāo)記開頭的 dna 序列中 cpg 島的分布相同，目的是使我們的實驗數(shù) 據(jù)盡可能符合原來 cpg 島在 dna 序列中的真實數(shù)據(jù)分布概率。2052.2評價手段為了對系統(tǒng)進行準(zhǔn)確客觀的評價，對每一條測試 dna 序列，我們使用以下 6 種評估指 標(biāo)對 cpg-discover 系統(tǒng)識別出的 cpg 島位置和真實數(shù)據(jù)庫中標(biāo)定的 cpg 島位置進行對比。 這 6 種評估指標(biāo)都在生物信息學(xué)領(lǐng)域有廣泛的應(yīng)用，它們的表示公式如下：(ntp+ntn)（1） accuracy (ac)：(ntp+nfn+nfp+ntn)ntp210（2） sensitivity (sn)：（3） sepcificity (sp)：ntp+nfn ntn ntn+nfpntp（4） positive predictive value (ppv)：（5） performance coefficient (pc):ntp+nfpntp（6） correlation coefficient (cc)：ntp+nfp+nfn215220225ntp*ntn-nfn*nfp(ntp+nfn) * (ntn+nfp) * (ntp+nfp) * (ntn+nfn)其中各個參數(shù)的定義如下：lntp：是 cpg 島內(nèi)的堿基且被正確的預(yù)測為 cpg 島內(nèi)堿基的堿基個數(shù)。lnfn: 是 cpg 島內(nèi)的堿基且被錯誤的預(yù)測為非 cpg 島內(nèi)堿基的堿基個數(shù)。lnfp: 不是 cpg 島內(nèi)的堿基且被錯誤的預(yù)測為 cpg 島內(nèi)堿基的堿基個數(shù)。lntn: 不是 cpg 島內(nèi)的堿基被正確的預(yù)測為非 cpg 島內(nèi)堿基的堿基個數(shù)。2.3實驗結(jié)果與分析在下面的實驗過程中，根據(jù) cpg 島的定義，參考其它軟件給出的可選參數(shù)，我們選擇 了如下的參數(shù)作為我們的篩選參數(shù)：長度在 140bp 以上，g+c 含量大于 60。為了研究訓(xùn)練數(shù)據(jù)集的大小對系統(tǒng)測試結(jié)果的影響，我們分別把“ab”和“bc”開頭的dna 序列訓(xùn)練數(shù)據(jù)集分成 7 組不同大小的訓(xùn)練數(shù)據(jù)集，即分別取200, 400, 600, 800, 1000,1200, 1400條訓(xùn)練數(shù)據(jù)來建立 hmm 系統(tǒng)，然后使用預(yù)先留出的 200 條測試 dna 序列來評 估系統(tǒng)的性能參考。圖 2 和圖 3 分別顯示對以“ab”和“bc”開頭的不同訓(xùn)練數(shù)據(jù)集對 6 個系 統(tǒng)性能評估指標(biāo)的結(jié)果。230235圖 2 ab 序列的測試結(jié)果fig.2 the test result of nda sequence(ab)圖 3 bc 序列的測試結(jié)果fig.3 the test result of nda sequence(bc)從圖 2 和圖 3 顯示的結(jié)果可以看出，對每一種性能評價指標(biāo)來說，隨著訓(xùn)練數(shù)據(jù)數(shù)量的增加，系統(tǒng)的性能指標(biāo)的均值總體是增加趨勢，而當(dāng)訓(xùn)練數(shù)據(jù)集包含大約 1000 條數(shù)據(jù)的時240245250255260265270275候，系統(tǒng)的各個性能指標(biāo)趨于穩(wěn)定。這個結(jié)果表明，當(dāng)系統(tǒng)在接受 1000 條訓(xùn)練數(shù)據(jù)之后，即可以輸出相對穩(wěn)定的測試結(jié)果。此外，從圖 2 和圖 3 的對比結(jié)果還可以看出，“bc”開頭的 dna 序列數(shù)據(jù)測試結(jié)果相對 于“ab”開頭的數(shù)據(jù)的測試結(jié)果要更加穩(wěn)定。通過對 dna 序列數(shù)據(jù)的分析，我們認(rèn)為可能的 原因如下：l首先，從訓(xùn)練數(shù)據(jù)來看，“ab”開頭的數(shù)據(jù)中，每條 dna 序列中的 cpg 島的數(shù)量波動 范圍集中在為 17，而以“bc”開頭的數(shù)據(jù)中，每條 dna 序列中的 cpg 島的數(shù)量波動范 圍集中在 14。l其次，“ab”開頭的 dna 序列平均長度為 5924bp, “bc”開頭的 dna 序列平均長度為1985bp。因此，相對“ab”開頭的 dna 序列，“bc”開頭的數(shù)據(jù)本身長度短，cpg 島數(shù)量波動小，比“ab”開頭的數(shù)據(jù)更加規(guī)范和集中，因而系統(tǒng)的波動性也更小，表現(xiàn)出較為平穩(wěn)的結(jié)果趨勢。 從上述實驗結(jié)果可以看出，cpg-discover 系統(tǒng)的 sensitivity 性能很高，最高達到了 97%，即使選擇不同訓(xùn)練數(shù)據(jù)集大小，系統(tǒng)的平均 sensitivity 值也超過 90%，尤其是“bc”開頭的較短 dna 序列數(shù)據(jù)。從 sensitivity 的定義來看，cpg 島內(nèi) 90%的堿基都能被正確識別標(biāo)定 出來，也證明 cpg-discover 系統(tǒng)能達到較為滿意的從人類 dna 序列中自動識別 cpg 島的 作用。此外，在最好的情況下，cpg-discover 系統(tǒng)的 accuracy 可以達到近 80%的正確率， 而且絕大部分的 accuracy 也超過了 50%。2.4與其它工具的對比為了進一步全面地評估系統(tǒng)的性能，我們隨機地從測試樣本中選擇若干條 dna 序列， 將 cpg-discover 系統(tǒng)與前面提到過的基于滑動窗口和聚集技術(shù)的工具進行了對比實驗。我 們選取了目前較為有名的 cpg 島識別工具或者系統(tǒng)，例如基于滑動窗口技術(shù)的 cpgprod4 和 cpgie5，基于聚集技術(shù)的工具 cpgcluster6和 cpgif7來進行比較，這些工具軟件涵 蓋了當(dāng)前廣泛采用的算法，并且在 cpg 島挖掘方面公認(rèn)有較好的性能。表 5 列出了部分 dna 序列分別使用這 5 個系統(tǒng)工具預(yù)測的 cpg 島位置和真實 cpg 島位置的對比結(jié)果此外，為了更全面比較 cpg-discover 與其它同類軟件的結(jié)果，我們對其他軟件的預(yù)測 結(jié)果使用 2.2 節(jié)中的 5 種評價指標(biāo)進行衡量比較。表 6 以 ab046787 和 bc011652 這兩條 dna 序列為例，列出了各個系統(tǒng)對這兩條 dna 序列進行識別的 5 個性能評價指標(biāo)的結(jié)果。從表 5 和表 6 的結(jié)果可以看出，當(dāng)相鄰的 cpg 島的間隔較小的時，以上五個軟件都不 能準(zhǔn)確確定 cpg 島的起始邊界位置，cpgif 和 cpgpord 傾向于將間隔較小的兩個 cpg 島 合并成一個 cpg 島，而 cpg-discover 則不能將其間隔準(zhǔn)確定位（如 bc008956、bc017346）， 而 cpgie 和 cpgcluster 出現(xiàn)了較多無法預(yù)測出 cpg 島的情況。cpgcluster 較其他軟件而言有很好的 specificity 和 positive predictive value，這表明 cpgcluster 識別出來的 cpg 島極可能是真正的 cpg 島。但是 cpgcluster 的 sensitivity 卻與 其他軟件有很大的差距，這說明對于很多是 cpg 島的區(qū)域，cpgcluster 很可能不能識別出 來而將其判成了非 cpg 島的區(qū)域。原因是 cpgcluster 的識別結(jié)果與輸入序列有關(guān)，如果待 測序列的非島區(qū)域較大，則會導(dǎo)致其 cpg 島區(qū)域的 p-value 偏低，因此造成了 cpgcluster 將 cpg 島的區(qū)域判斷成非 cpg 島的區(qū)域的情況。cpgprod 的 specificity 相對于其他軟件而 言比較低。而 cpgie、cpgif 和 cpg-discover 的各個評價參數(shù)都相對較好。這個結(jié)果非常 有意思，因為這三個工具正是分別采用了滑動窗口，聚集技術(shù)和 hmm 模型，這也解釋了為 什么目前采用不同技術(shù)的 cpg 島識別系統(tǒng)同時都在得到廣泛應(yīng)用中。280表 5 cpg-discover 系統(tǒng)與其它四個軟件工具的比較結(jié)果。選用參數(shù)長度在 140bp 以上，g+c 含量大于60 ，并且實際 cpg 含量與期望 cpg 含量的比值(obscpgexpc g)大于 06 ，其中 cpgcluster 工具的 參數(shù)取默認(rèn)值為 distance threshold 50 p-value threshold 1e-5。tab. 5 the comparison of cpg-discover and other tools. the length of over 140bp. g+c 60%. obscpg/expcpg 0.6. the default parameter of cpgcluster is distance threshold 50 p-value threshold 1e-5序列標(biāo)號bc011973bc008956bc014043ab037841ab046787ab007883長度 (bp)161812333079451659906246真實 cpg 島的起始和結(jié)束位置49507483764961024804100414751989239526104839847264190408cpgcluster135191813900611725217576639cpgif125327104979210361324271522525552464811cpgprod1883113087352628180317492031200cpgie16071522140556485870cpg-discover17051507577938130242873179696810683038137512551780長度 (bp)161812333079451659906246285表 6 cpg-discover 與其它同類軟件使用 5 個性能評價指標(biāo)的對比結(jié)果tab. 6 the performance evaluation of cpg-discover and other tools.ab046787bc011652系統(tǒng)snspppvpcccsnspppvpccccpgcluster0.7840.9930.8030.6580.7880.297110.2970.486cpgif00.962000.0380.9870.9760.9400.9290.953cpgprod10.9130.3020.3020.52710.6920.5510.5510.620cpgie10.9210.3230.3230.54710.9160.8180.8180.869cpg-discover0.9590.9920.8200.7920.88510.9320.8480.8480.892290295300從上述分析比較可以看出，cpg-discover 系統(tǒng)在預(yù)測 cpg 島的方面有很好的表現(xiàn)。首先，在 specificity 和 sensitivity 方面都比較優(yōu)秀，這表明 cpg-discover 系統(tǒng)既可以識別出絕 大部分 cpg 島，同時又有很好的識別精度，可以比較精確的定位 cpg 島的位置。其次， cpg-discover 系統(tǒng)是基于訓(xùn)練數(shù)據(jù)的自動識別系統(tǒng)，也就是說，給定標(biāo)記的 cpg 島訓(xùn)練數(shù) 據(jù)，系統(tǒng)一經(jīng)訓(xùn)練建立，可以應(yīng)用到各種 cpg 島的自動識別中，大大減少人工干預(yù)和勞動。3結(jié)論本文中，我們實現(xiàn)了一個基于隱馬爾科夫模型（hmm ）的 cpg 島自動識別系 統(tǒng)cpg-discover，在歐洲分子生物學(xué)實驗室(embl)的核苷序列數(shù)據(jù)庫發(fā)布的人類 dna 序列數(shù)據(jù)集上