一個(gè)基于隱馬爾可夫模型和生物知識(shí)修正.doc

精品論文一個(gè)基于隱馬爾可夫模型和生物知識(shí)修正的 cpg 島識(shí)別系統(tǒng)徐瑜，蘭曼5（華東師范大學(xué)信息科學(xué)技術(shù)學(xué)院， 上海 200241）摘要： cpg 島的存在能識(shí)別某些基因的啟動(dòng)子，而且 cpg 島的異常甲基化多與人類腫瘤的 發(fā)生有關(guān)，因此 cpg 島的識(shí)別在生物基因組測序中很重要。本文實(shí)現(xiàn)了一個(gè)基于隱馬爾可夫模型（hmm）和后期生物知識(shí)修正的 cpg 島識(shí)別系統(tǒng)，在 embl 的 dna 序列數(shù)據(jù)集上10進(jìn)行系統(tǒng)性能測試的結(jié)果顯示，該系統(tǒng)對(duì)于 cpg 島有較好的識(shí)別能力，同時(shí)又比較精確地 定位 cpg 島的位置。與其它常用 cpg 島識(shí)別工具的對(duì)比實(shí)驗(yàn)結(jié)果表明，該系統(tǒng)的識(shí)別準(zhǔn)確性不亞于其它軟件。此外，該系統(tǒng)的優(yōu)點(diǎn)是一經(jīng)訓(xùn)練，可以用于自動(dòng)識(shí)別。關(guān)鍵詞：隱馬爾可夫模型；機(jī)器學(xué)習(xí)；cpg 島識(shí)別；知識(shí)修正中圖分類號(hào)：請(qǐng)查閱中國圖書館分類法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）的一個(gè)區(qū)域。通常，cpg島的長度為幾百到幾千個(gè)堿基對(duì)（nucleotides，單位bp）。在人的基因組中，如果雙堿基對(duì)cg出現(xiàn)，則c通常被甲基化。并 且，甲基化的c很快會(huì)突變成t，因此基因組中cpg島非常少。然而，在基因的起始位置（啟 動(dòng)子），因?yàn)楣δ艿谋Ｊ匦裕湫蛄泻苌偻蛔?，在哺乳類基因中的啟?dòng)子上含有約40%的cpg 島，人類基因中含有約70%的cpg島，因此，cpg島的存在能識(shí)別某些基因的啟動(dòng)子，并可40作為限制酶的辨識(shí)位置。從已知的dna序列統(tǒng)計(jì)來看，幾乎所有的管家基因（house-keeping genes）及約40%的組織特異性基因（tissue-specific genes）的附近(通常是在5上游區(qū)域)都含基金項(xiàng)目：教育部博士點(diǎn)基金(20090076120029)作者簡介：徐瑜（1988-）， 女，碩士研究生，自然語言處理，文本挖掘通信聯(lián)系人：蘭曼（1974-）， 女，副教授，自然語言處理，數(shù)據(jù)挖掘，生物文本挖掘等. e-mail:- 2 -有非甲基化的cpg島，其序列可能包括基因轉(zhuǎn)錄的啟動(dòng)子的第一外顯子1。因此，cpg島的 識(shí)別在生物基因組測序工作中是非常重要的。此外，人類腫瘤的發(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)記，因此，啟動(dòng)子區(qū)域的cpg島的異常甲基化也是識(shí) 別癌癥的重要標(biāo)志之一。針對(duì)cpg島的定義，1987年，gardiner-garden和formmer2認(rèn)為長度在200bp以上，g+c 含量大于50，并且實(shí)際cpg含量與期望cpg含量的比值(obscpgexpc g)大于0.6的區(qū)域 即為cpg島。2002年，daiya takai 和 peter a. jones3的研究認(rèn)為，長度在500bp以上，g+c50含量大于55，并且實(shí)際cpg含量與期望cpg含量的比值(obscpgexpc g)大于0.65的區(qū)域 更可能和基因的5區(qū)域有關(guān)。因此，根據(jù)研究的不同的目的，生物學(xué)研究人員可以采用不同 的參數(shù)進(jìn)行cpg島預(yù)測。dna序列中的cpg島既可以通過生物學(xué)實(shí)驗(yàn)的方法識(shí)別，也可以通過計(jì)算機(jī)來識(shí)別可能 的cpg島。由于dna序列的多樣性，如果借助于生物信息和計(jì)算機(jī)工具首先識(shí)別出潛在的55cpg島，就可以大大減少生物實(shí)驗(yàn)的成本和盲目性，極大地提高生物實(shí)驗(yàn)的時(shí)間效率和針對(duì) 性。近年來，利用計(jì)算機(jī)開發(fā)cpg島識(shí)別工具，吸引了許多生物學(xué)領(lǐng)域和傳統(tǒng)計(jì)算機(jī)領(lǐng)域研 究人員的興趣，他們已經(jīng)開發(fā)出許多有用的工具來幫助生物學(xué)家提高工作效率。根據(jù)這些工 具所采用的方法，我們將它們大致分為兩類。第一類方法是采用傳統(tǒng)的滑動(dòng)窗口（sliding window），這些工具包括cpgis3，cpgprod4和cpgie5等?；瑒?dòng)窗口法從頭至尾對(duì)dna60序列進(jìn)行掃描（每次移動(dòng)1n位），分別計(jì)算窗口內(nèi)序列是cpg島 (用“+”表示) 和不是cpg 島 (用“-”表示) 的概率值（通常采用log-odds ratio），如果高于設(shè)定的域值，則認(rèn)為窗口內(nèi) 的序列是cpg島，否則認(rèn)為不是。采用滑動(dòng)窗口法的不足之處在于：（1）窗口的大小很難 確定，通常窗口的大小直接影響cpg島的長度和數(shù)量，如果窗口過大，一些長度較短的cpg 島會(huì)被合并在一起，從而被預(yù)測成一個(gè)大的cpg島，降低識(shí)別的準(zhǔn)確率；（2）窗口只能向65一個(gè)方向移動(dòng)，對(duì)于之前錯(cuò)誤的判斷結(jié)果，不能在后面得到發(fā)現(xiàn)和修改，從而影響識(shí)別的精 度；（3）系統(tǒng)運(yùn)行時(shí)間較長。第二類方法是采用聚集（clustering）的技術(shù)，這些工具包括 cpgcluster6和cpgif7等。這類工具的基本思想是基于cpg島定義中的密度或物理距離等生 物學(xué)屬性，不斷迭代地掃描序列，并將符合這些屬性的相鄰cpg島合并為大的cpg島簇。聚 集方法的提出是為了解決滑動(dòng)窗口方法的不足，然而它帶有如下的缺陷：（1）識(shí)別的結(jié)果70依賴于掃描序列的順序，而不同的掃描順序則會(huì)識(shí)別出不同的cpg島，因此系統(tǒng)穩(wěn)定性不足；（2）識(shí)別的敏感度較低，傾向于識(shí)別長度較長的cpg島序列。 與上述兩類方法不同，在本文，我們嘗試使用隱馬爾可夫模型（hmm）建立一個(gè)自動(dòng)識(shí)別cpg島的系統(tǒng)，并結(jié)合生物知識(shí)對(duì)預(yù)測結(jié)果進(jìn)行后期修正，從而達(dá)到較為理想的識(shí)別預(yù) 測結(jié)果。隱馬爾可夫模型是一種機(jī)器學(xué)習(xí)的方法，它廣泛應(yīng)用于序列分析和模式識(shí)別的研究75中，例如語音識(shí)別、手寫識(shí)別，命名實(shí)體識(shí)別等。這類應(yīng)用的一個(gè)重要特點(diǎn)是，數(shù)據(jù)的序列 性。這種序列可以被認(rèn)為是后一狀態(tài)依賴于前一狀態(tài)的時(shí)間序列，即具有馬爾可夫性。正是 因?yàn)榭紤]到dna序列（堿基對(duì)）也可以被看做是具有馬爾可夫性的有序數(shù)據(jù)，我們在本文 中利用hmm模型實(shí)現(xiàn)了一個(gè)自動(dòng)從人類dna序列中識(shí)別cpg島的系統(tǒng)，即cpg-discover系 統(tǒng)。80本文中，我們首先介紹結(jié)合hmm模型和生物知識(shí)修正的cpg-discover系統(tǒng)的框架結(jié)構(gòu)。 然后在歐洲分子生物學(xué)實(shí)驗(yàn)室（embl）發(fā)布的人類基因中cpg島的數(shù)據(jù)集上測試本系統(tǒng)的 性能，并與其它常用的cpg島識(shí)別工具（基于滑動(dòng)窗口和聚集技術(shù)）進(jìn)行對(duì)比實(shí)驗(yàn)，深入分析和討論實(shí)驗(yàn)結(jié)果。最后我們總結(jié)這一工作和未來的方向。1cpg-discover 系統(tǒng)框架85cpg-discover 系統(tǒng)的框架結(jié)構(gòu)主要包括三個(gè)子系統(tǒng)：（1）建立 hmm 模型的參數(shù)系統(tǒng)；（2）基于 hmm 模型建立 cpg 島識(shí)別系統(tǒng)；（2）基于生物知識(shí)的后期識(shí)別修正系統(tǒng)。下 面我們分別詳細(xì)介紹這三個(gè)子系統(tǒng)的工作。1.1建立 hmm 模型的參數(shù)系統(tǒng)這個(gè)子系統(tǒng)完成對(duì) hmm 模型的求解過程，即實(shí)現(xiàn)參數(shù)初始化、參數(shù)逼近和參數(shù)投票這90三個(gè)步驟，建立起 hmm 模型的參數(shù)系統(tǒng)。在參數(shù)初始化階段，針對(duì) 4 種堿基（a，c，t， g）分別在 cpg 島內(nèi)外的狀態(tài)，我們標(biāo)記此 hmm 模型包括 8 個(gè)狀態(tài)，分別是 a+, c+, g+, t+, a, c, g,t，標(biāo)號(hào)中“+”和“”分別表示此核苷酸在和不在 cpg 島上。這 8 個(gè)狀態(tài)之間的 轉(zhuǎn)移方式如圖 1 所示。這 8 個(gè)狀態(tài)間的狀態(tài)轉(zhuǎn)移概率可以從有標(biāo)識(shí)的訓(xùn)練序列中得到，后面 會(huì)詳細(xì)敘述。本文中，我們基于歐洲分子生物學(xué)實(shí)驗(yàn)室（embl）核苷序列數(shù)據(jù)庫進(jìn)行實(shí)驗(yàn)95建立 hmm 模型，因此，在以下三個(gè)參數(shù)求解的過程中，我們以 id 為 bc008880 的人類 dna序列為例來進(jìn)行說明。c+g+a+t+a-t- c-g- 11 -1001.1.1參數(shù)初始化圖 1hmm 模型中 8 個(gè)狀態(tài)之間的轉(zhuǎn)移圖示fig. 1 the transform of 8 states in hmm105110這個(gè)步驟完成 hmm 模型中狀態(tài)轉(zhuǎn)移概率矩陣（a）、發(fā)射概率矩陣（b）和狀態(tài)概率 矩陣（）的初始化，作為下一步參數(shù)逼近的輸入種子。三個(gè)概率矩陣的初始化步驟如下：（1）初始化狀態(tài)轉(zhuǎn)移概率矩陣 a：對(duì)于每個(gè)訓(xùn)練樣本序列，我們采用基于長度為 2 的 滑動(dòng)窗口的方法來初始化參數(shù)，即通過計(jì)算從第一個(gè)狀態(tài)轉(zhuǎn)移到第二個(gè)狀態(tài)的出現(xiàn)頻數(shù)，來 計(jì)算轉(zhuǎn)移概率。（2）初始化發(fā)射概率矩陣 b：假設(shè)一個(gè)狀態(tài)到一個(gè)符號(hào)的發(fā)射概率都為 0 或 1。（3）初始化狀態(tài)概率矩陣 ：每個(gè)狀態(tài)的初始化概率由序列中各個(gè)狀態(tài)的出現(xiàn)頻數(shù)計(jì) 算。參數(shù)完成初始化后，我們得到初始化后的各個(gè)概率矩陣。表 1 顯示了人類 dna 序列bc008880 經(jīng)過初始化之后的三個(gè)概率矩陣的分布結(jié)果。1.1.2參數(shù)逼近在第二步中，我們采用 baum-welch 算法來重新估計(jì)基于訓(xùn)練序列的概率矩陣。上一步 的三個(gè)初始化概率矩陣作為 baum-welch 算法的輸入，輸出期望的轉(zhuǎn)移和發(fā)射概率矩陣，再115將得到的期望概率代替舊的參數(shù)，這個(gè)過程不斷迭代進(jìn)行直至達(dá)到收斂。在每次迭代過程中，模型概率的準(zhǔn)確度都得到提高，直到一個(gè)極限概率，整個(gè)迭代過程最終收斂于一個(gè)局部最優(yōu) 解。參數(shù)估計(jì)后，對(duì)每個(gè)訓(xùn)練序列，我們都得到三個(gè)逼近后的概率矩陣。表 2 顯示了人類 dna 序列 bc008880 經(jīng)過 baum-welch 近似之后的三個(gè)概率矩陣的分布結(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)過初始化之后的三個(gè)概率矩陣的結(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 逼近之后的三個(gè)概率矩陣的結(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ù)投票至此，對(duì)每個(gè) dna 訓(xùn)練序列，我們都得到三個(gè)逼近后的概率矩陣。為了整合各個(gè)不同 長度的訓(xùn)練序列的概率矩陣，我們采取投票（vote）平均技術(shù)，就是說，對(duì)每個(gè)概率矩陣， 根據(jù)每個(gè)訓(xùn)練序列的長度對(duì)其進(jìn)行歸一化（normalization），每個(gè) dna 訓(xùn)練序列的權(quán)重計(jì)算公式為： w = ni，其中，n 表示第 i 個(gè)訓(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 模型的三個(gè)概率矩陣的分布結(jié)果（關(guān)于數(shù)據(jù)集見后面 2.1 節(jié)）。1401451.2基于 hmm 模型建立 cpg 島識(shí)別系統(tǒng)經(jīng)過求解和加權(quán)平均之后得到的 hmm 模型，可以通過 viterbi 解碼算法來尋找能夠產(chǎn) 生給定符號(hào)序列的最大似然狀態(tài)序列，即對(duì)每個(gè)測試 dna 序列中的每個(gè)堿基對(duì)（atcg）， viterbi 解碼算法能夠輸出每個(gè)堿基對(duì)對(duì)應(yīng)的狀態(tài)為”+”或者”-”（其中“+”標(biāo)識(shí)堿基對(duì)在 cpg 島內(nèi)，”-”表示堿基對(duì)在 cpg 島之外），并使得標(biāo)識(shí)后的 dna 序列輸出狀態(tài)序列為最優(yōu)解。表 3 使用 1400 條以“bc”開頭的人類 dna 序列作為訓(xùn)練數(shù)據(jù)集，經(jīng)過加權(quán)平均（投票）之后的 hmm模型的三個(gè)概率矩陣分布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基于生物知識(shí)的后期識(shí)別修正系統(tǒng)由viterbi 算法求解出的dna序列輸出狀態(tài)是基于訓(xùn)練模型計(jì)算而得，系統(tǒng)輸出的狀態(tài) 序列并不一定滿足生物學(xué)上對(duì)cpg島的定義。為了進(jìn)一步提高結(jié)果的精度，我們結(jié)合cpg島 的生物定義，如cpg島中c和g的概率、相鄰cpg片段間距以及cpg片段的長度閾值等，對(duì) hmm模型識(shí)別dna序列產(chǎn)生的誤差進(jìn)行后期的知識(shí)修正。我們對(duì)cpg島作如下的限定：至 少包含140個(gè)堿基對(duì)的長度，區(qū)域內(nèi)gc所占含量超過60%，且cpg的觀察值預(yù)測值比例必 須高于0.6（根據(jù)研究的不同目的，可以在軟件系統(tǒng)中自行修改這些限定的值）。這個(gè)后期 識(shí)別修正系統(tǒng)包含下面三個(gè)模塊：（1）相鄰cpg島的合并模塊：由于計(jì)算誤差，系統(tǒng)經(jīng)常把cpg島中少量的核苷酸（通 常長度在1bp到15bp之間）錯(cuò)誤地識(shí)別為非cpg島上的核苷酸。為了解決這個(gè)問題，我們對(duì) 相鄰cpg島設(shè)定最小距離限定。經(jīng)過多次實(shí)驗(yàn)，這個(gè)最小距離的閾值設(shè)定為20，即，如果識(shí) 別出的兩個(gè)相鄰cpg島間的距離小于20，系統(tǒng)就修改這些間隔中的非cpg島核苷酸狀態(tài)，使 之狀態(tài)轉(zhuǎn)變?yōu)閏pg島上核苷酸，即合并這兩個(gè)相鄰cpg島為一條長cpg島。合并后的cpg島 需要進(jìn)一步經(jīng)過后面兩個(gè)模塊的修正。（2）密度檢測過濾模塊：生物學(xué)上，cpg島通常含有較高的g/c（或c/g）含量(密度高 于60%)，如果前一步的錯(cuò)誤合并或者系統(tǒng)識(shí)別的誤差，可能產(chǎn)生不滿足c/g含量要求的cpg 島序列。因此，對(duì)識(shí)別出的不能滿足最低密度要求的cpg島序列，即區(qū)域內(nèi)gc所占含量至 少達(dá)到60%，這個(gè)密度檢測過濾模塊將改變這些cpg島中所有核苷酸狀態(tài)為非cpg島狀態(tài)。（3）長度檢測過濾模塊：除了（2）中的密度約束模塊外，這個(gè)系統(tǒng)還進(jìn)行了長度約束， 即根據(jù)cpg島的長度定義，去除長度少于140個(gè)堿基對(duì)的cpg島序列，即，對(duì)識(shí)別出的長度175180小于140bp的cpg島將被改變?yōu)榉莄pg島狀態(tài)。經(jīng)過以上三個(gè)后期修正模塊之后的輸出序列狀態(tài)，才成為cpg-discover系統(tǒng)的最終輸入結(jié)果。表4列出了人類dna序列bc009465在修正前后的對(duì)比結(jié)果，其中491255分別表示cpg島序列的起始位置為第49個(gè)堿基對(duì)和結(jié)束位置為第1255個(gè)堿基對(duì)，這6個(gè)評(píng)估指標(biāo)的含義在后面2.2節(jié)中有詳述。從表4可以看到，基于生物知識(shí)修正后的識(shí)別結(jié)果除了correlation coefficient指標(biāo)沒有變化外，其它的指標(biāo)都有提高，特別是specificity，從原來的33.76%增至了51.68%（提高60.15%）。原來預(yù)測位置在1284和2961509的兩條cpg島因?yàn)橄喔糁挥?2個(gè)堿基對(duì)（bp），小于最小間隔閾值20bp，因而被修正合并為一條長的cpg島；而原來預(yù)測位置在15841695的cpg島序列長度只有120bp，在后處理模塊中由于長度過短，不滿足 生物學(xué)上cpg島的長度要求，因而被過濾除去。表 4 人類 dna 序列 bc009465 基于生物知識(shí)的修正前后識(shí)別結(jié)果的對(duì)比tab. 4 the comparison of domain knowledge revision (bc009456)評(píng)測指標(biāo)未修正的結(jié)果修正后的結(jié)果真實(shí)cpg島的始末端位置491255491255系統(tǒng)識(shí)別的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實(shí)驗(yàn)部分為了對(duì) cpg-discover 系統(tǒng)的性能做一個(gè)準(zhǔn)確而全面的分析和比較，我們從歐洲分子生 物學(xué)實(shí)驗(yàn)室(embl)的核苷序列數(shù)據(jù)庫（http:/www.ebi.ac.uk/embl/）中下載了關(guān)于人類基因 的 cpg 島數(shù)據(jù)庫 emb173hum（ftp:/ftp.ebi.ac.uk/pub/databases/cpgisle），并在這個(gè)數(shù)據(jù)集上 測試了本系統(tǒng)的性能。這個(gè)數(shù)據(jù)庫是由歐洲生物信息中心負(fù)責(zé)維護(hù)，與日本基因數(shù)據(jù)庫（ddbj）和美國國立生物技術(shù)信息中心（ncbi）相互合作，他們之間每天都要交換最新的 核苷序列數(shù)據(jù)。此外，為了比較 cpg-discover 系統(tǒng)與其它常用的 cpg 島識(shí)別工具（基于滑 動(dòng)窗口和聚集技術(shù)）的性能，我們還在這個(gè)數(shù)據(jù)集上選擇部分 dna 序列進(jìn)行了對(duì)比實(shí)驗(yàn)。2.1數(shù)據(jù)來源emb173hum 數(shù)據(jù)庫包含了當(dāng)前 embl 數(shù)據(jù)庫發(fā)布的已標(biāo)定 cpg 島位置的人類 dna 序 列，共有 233,004 條 dna 序列和 142,325 個(gè) cpg 島，其中只有 61,051 條 dna 序列中包含 有一個(gè)或者若干個(gè) cpg 島，平均每條 dna 序列含有 2.33 個(gè) cpg 島。我們下載了這個(gè)數(shù)據(jù) 庫中以“ab”標(biāo)記開頭的 1,955 條 dna 序列和以“bc”標(biāo)記開頭的 6,124 條 dna 序列，除去 數(shù)據(jù)庫中的空序列和有錯(cuò)誤信息的序列之后，我們得到近 8,000 條 dna 序列。其中，以“ab” 標(biāo)記開頭的 dna 序列中平均每條 dna 序列含有 1.70 個(gè) cpg 島，以“bc”標(biāo)記開頭的 dna 序列中平均每條 dna 序列含有 1.34 個(gè) cpg 島。為了使實(shí)驗(yàn)結(jié)果真實(shí)合理有效，我們從以 “ab”和“bc”標(biāo)記開頭的 dna 序列中分別選擇了 1400 條訓(xùn)練數(shù)據(jù)和 200 條測試數(shù)據(jù)（共有2800 條訓(xùn)練數(shù)據(jù)和 400 條測試數(shù)據(jù)），選擇的原則是這個(gè)數(shù)據(jù)子集里面的 cpg 島的分布分 別與原來“ab”和“bc”標(biāo)記開頭的 dna 序列中 cpg 島的分布相同，目的是使我們的實(shí)驗(yàn)數(shù) 據(jù)盡可能符合原來 cpg 島在 dna 序列中的真實(shí)數(shù)據(jù)分布概率。2052.2評(píng)價(jià)手段為了對(duì)系統(tǒng)進(jìn)行準(zhǔn)確客觀的評(píng)價(jià)，對(duì)每一條測試 dna 序列，我們使用以下 6 種評(píng)估指 標(biāo)對(duì) cpg-discover 系統(tǒng)識(shí)別出的 cpg 島位置和真實(shí)數(shù)據(jù)庫中標(biāo)定的 cpg 島位置進(jìn)行對(duì)比。 這 6 種評(píng)估指標(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)其中各個(gè)參數(shù)的定義如下：lntp：是 cpg 島內(nèi)的堿基且被正確的預(yù)測為 cpg 島內(nèi)堿基的堿基個(gè)數(shù)。lnfn: 是 cpg 島內(nèi)的堿基且被錯(cuò)誤的預(yù)測為非 cpg 島內(nèi)堿基的堿基個(gè)數(shù)。lnfp: 不是 cpg 島內(nèi)的堿基且被錯(cuò)誤的預(yù)測為 cpg 島內(nèi)堿基的堿基個(gè)數(shù)。lntn: 不是 cpg 島內(nèi)的堿基被正確的預(yù)測為非 cpg 島內(nèi)堿基的堿基個(gè)數(shù)。2.3實(shí)驗(yàn)結(jié)果與分析在下面的實(shí)驗(yàn)過程中，根據(jù) cpg 島的定義，參考其它軟件給出的可選參數(shù)，我們選擇 了如下的參數(shù)作為我們的篩選參數(shù)：長度在 140bp 以上，g+c 含量大于 60。為了研究訓(xùn)練數(shù)據(jù)集的大小對(duì)系統(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 序列來評(píng) 估系統(tǒng)的性能參考。圖 2 和圖 3 分別顯示對(duì)以“ab”和“bc”開頭的不同訓(xùn)練數(shù)據(jù)集對(duì) 6 個(gè)系 統(tǒng)性能評(pí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é)果可以看出，對(duì)每一種性能評(píng)價(jià)指標(biāo)來說，隨著訓(xùn)練數(shù)據(jù)數(shù)量的增加，系統(tǒng)的性能指標(biāo)的均值總體是增加趨勢，而當(dāng)訓(xùn)練數(shù)據(jù)集包含大約 1000 條數(shù)據(jù)的時(shí)240245250255260265270275候，系統(tǒng)的各個(gè)性能指標(biāo)趨于穩(wěn)定。這個(gè)結(jié)果表明，當(dāng)系統(tǒng)在接受 1000 條訓(xùn)練數(shù)據(jù)之后，即可以輸出相對(duì)穩(wěn)定的測試結(jié)果。此外，從圖 2 和圖 3 的對(duì)比結(jié)果還可以看出，“bc”開頭的 dna 序列數(shù)據(jù)測試結(jié)果相對(duì) 于“ab”開頭的數(shù)據(jù)的測試結(jié)果要更加穩(wěn)定。通過對(duì) dna 序列數(shù)據(jù)的分析，我們認(rèn)為可能的 原因如下：l首先，從訓(xùn)練數(shù)據(jù)來看，“ab”開頭的數(shù)據(jù)中，每條 dna 序列中的 cpg 島的數(shù)量波動(dòng) 范圍集中在為 17，而以“bc”開頭的數(shù)據(jù)中，每條 dna 序列中的 cpg 島的數(shù)量波動(dòng)范 圍集中在 14。l其次，“ab”開頭的 dna 序列平均長度為 5924bp, “bc”開頭的 dna 序列平均長度為1985bp。因此，相對(duì)“ab”開頭的 dna 序列，“bc”開頭的數(shù)據(jù)本身長度短，cpg 島數(shù)量波動(dòng)小，比“ab”開頭的數(shù)據(jù)更加規(guī)范和集中，因而系統(tǒng)的波動(dòng)性也更小，表現(xiàn)出較為平穩(wěn)的結(jié)果趨勢。 從上述實(shí)驗(yàn)結(jié)果可以看出，cpg-discover 系統(tǒng)的 sensitivity 性能很高，最高達(dá)到了 97%，即使選擇不同訓(xùn)練數(shù)據(jù)集大小，系統(tǒng)的平均 sensitivity 值也超過 90%，尤其是“bc”開頭的較短 dna 序列數(shù)據(jù)。從 sensitivity 的定義來看，cpg 島內(nèi) 90%的堿基都能被正確識(shí)別標(biāo)定 出來，也證明 cpg-discover 系統(tǒng)能達(dá)到較為滿意的從人類 dna 序列中自動(dòng)識(shí)別 cpg 島的 作用。此外，在最好的情況下，cpg-discover 系統(tǒng)的 accuracy 可以達(dá)到近 80%的正確率， 而且絕大部分的 accuracy 也超過了 50%。2.4與其它工具的對(duì)比為了進(jìn)一步全面地評(píng)估系統(tǒng)的性能，我們隨機(jī)地從測試樣本中選擇若干條 dna 序列， 將 cpg-discover 系統(tǒng)與前面提到過的基于滑動(dòng)窗口和聚集技術(shù)的工具進(jìn)行了對(duì)比實(shí)驗(yàn)。我 們選取了目前較為有名的 cpg 島識(shí)別工具或者系統(tǒng)，例如基于滑動(dòng)窗口技術(shù)的 cpgprod4 和 cpgie5，基于聚集技術(shù)的工具 cpgcluster6和 cpgif7來進(jìn)行比較，這些工具軟件涵 蓋了當(dāng)前廣泛采用的算法，并且在 cpg 島挖掘方面公認(rèn)有較好的性能。表 5 列出了部分 dna 序列分別使用這 5 個(gè)系統(tǒng)工具預(yù)測的 cpg 島位置和真實(shí) cpg 島位置的對(duì)比結(jié)果此外，為了更全面比較 cpg-discover 與其它同類軟件的結(jié)果，我們對(duì)其他軟件的預(yù)測 結(jié)果使用 2.2 節(jié)中的 5 種評(píng)價(jià)指標(biāo)進(jìn)行衡量比較。表 6 以 ab046787 和 bc011652 這兩條 dna 序列為例，列出了各個(gè)系統(tǒng)對(duì)這兩條 dna 序列進(jìn)行識(shí)別的 5 個(gè)性能評(píng)價(jià)指標(biāo)的結(jié)果。從表 5 和表 6 的結(jié)果可以看出，當(dāng)相鄰的 cpg 島的間隔較小的時(shí)，以上五個(gè)軟件都不 能準(zhǔn)確確定 cpg 島的起始邊界位置，cpgif 和 cpgpord 傾向于將間隔較小的兩個(gè) cpg 島 合并成一個(gè) cpg 島，而 cpg-discover 則不能將其間隔準(zhǔn)確定位（如 bc008956、bc017346）， 而 cpgie 和 cpgcluster 出現(xiàn)了較多無法預(yù)測出 cpg 島的情況。cpgcluster 較其他軟件而言有很好的 specificity 和 positive predictive value，這表明 cpgcluster 識(shí)別出來的 cpg 島極可能是真正的 cpg 島。但是 cpgcluster 的 sensitivity 卻與 其他軟件有很大的差距，這說明對(duì)于很多是 cpg 島的區(qū)域，cpgcluster 很可能不能識(shí)別出 來而將其判成了非 cpg 島的區(qū)域。原因是 cpgcluster 的識(shí)別結(jié)果與輸入序列有關(guān)，如果待 測序列的非島區(qū)域較大，則會(huì)導(dǎo)致其 cpg 島區(qū)域的 p-value 偏低，因此造成了 cpgcluster 將 cpg 島的區(qū)域判斷成非 cpg 島的區(qū)域的情況。cpgprod 的 specificity 相對(duì)于其他軟件而 言比較低。而 cpgie、cpgif 和 cpg-discover 的各個(gè)評(píng)價(jià)參數(shù)都相對(duì)較好。這個(gè)結(jié)果非常 有意思，因?yàn)檫@三個(gè)工具正是分別采用了滑動(dòng)窗口，聚集技術(shù)和 hmm 模型，這也解釋了為 什么目前采用不同技術(shù)的 cpg 島識(shí)別系統(tǒng)同時(shí)都在得到廣泛應(yīng)用中。280表 5 cpg-discover 系統(tǒng)與其它四個(gè)軟件工具的比較結(jié)果。選用參數(shù)長度在 140bp 以上，g+c 含量大于60 ，并且實(shí)際 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)號(hào)bc011973bc008956bc014043ab037841ab046787ab007883長度 (bp)161812333079451659906246真實(shí) cpg 島的起始和結(jié)束位置49507483764961024804100414751989239526104839847264190408cpgcluster135191813900611725217576639cpgif125327104979210361324271522525552464811cpgprod1883113087352628180317492031200cpgie16071522140556485870cpg-discover17051507577938130242873179696810683038137512551780長度 (bp)161812333079451659906246285表 6 cpg-discover 與其它同類軟件使用 5 個(gè)性能評(píng)價(jià)指標(biāo)的對(duì)比結(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)既可以識(shí)別出絕 大部分 cpg 島，同時(shí)又有很好的識(shí)別精度，可以比較精確的定位 cpg 島的位置。其次， cpg-discover 系統(tǒng)是基于訓(xùn)練數(shù)據(jù)的自動(dòng)識(shí)別系統(tǒng)，也就是說，給定標(biāo)記的 cpg 島訓(xùn)練數(shù) 據(jù)，系統(tǒng)一經(jīng)訓(xùn)練建立，可以應(yīng)用到各種 cpg 島的自動(dòng)識(shí)別中，大大減少人工干預(yù)和勞動(dòng)。3結(jié)論本文中，我們實(shí)現(xiàn)了一個(gè)基于隱馬爾科夫模型（hmm ）的 cpg 島自動(dòng)識(shí)別系 統(tǒng)cpg-discover，在歐洲分子生物學(xué)實(shí)驗(yàn)室(embl)的核苷序列數(shù)據(jù)庫發(fā)布的人類 dna 序列數(shù)據(jù)集上