一个基于隐马尔可夫模型和生物知识修正.doc

精品论文一个基于隐马尔可夫模型和生物知识修正的 cpg 岛识别系统徐瑜，兰曼5（华东师范大学信息科学技术学院， 上海 200241）摘要： cpg 岛的存在能识别某些基因的启动子，而且 cpg 岛的异常甲基化多与人类肿瘤的 发生有关，因此 cpg 岛的识别在生物基因组测序中很重要。本文实现了一个基于隐马尔可夫模型（hmm）和后期生物知识修正的 cpg 岛识别系统，在 embl 的 dna 序列数据集上10进行系统性能测试的结果显示，该系统对于 cpg 岛有较好的识别能力，同时又比较精确地 定位 cpg 岛的位置。与其它常用 cpg 岛识别工具的对比实验结果表明，该系统的识别准确性不亚于其它软件。此外，该系统的优点是一经训练，可以用于自动识别。关键词：隐马尔可夫模型；机器学习；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引言生物学上，cpg岛（cpg islands）是指dna上包含大量相邻的胞嘧啶（c）、鸟嘌呤（g），35以及使两者相连的磷酸酯键（p）的一个区域。通常，cpg岛的长度为几百到几千个碱基对（nucleotides，单位bp）。在人的基因组中，如果双碱基对cg出现，则c通常被甲基化。并 且，甲基化的c很快会突变成t，因此基因组中cpg岛非常少。然而，在基因的起始位置（启 动子），因为功能的保守性，其序列很少突变，在哺乳类基因中的启动子上含有约40%的cpg 岛，人类基因中含有约70%的cpg岛，因此，cpg岛的存在能识别某些基因的启动子，并可40作为限制酶的辨识位置。从已知的dna序列统计来看，几乎所有的管家基因（house-keeping genes）及约40%的组织特异性基因（tissue-specific genes）的附近(通常是在5上游区域)都含基金项目：教育部博士点基金(20090076120029)作者简介：徐瑜（1988-）， 女，硕士研究生，自然语言处理，文本挖掘通信联系人：兰曼（1974-）， 女，副教授，自然语言处理，数据挖掘，生物文本挖掘等. e-mail:- 2 -有非甲基化的cpg岛，其序列可能包括基因转录的启动子的第一外显子1。因此，cpg岛的 识别在生物基因组测序工作中是非常重要的。此外，人类肿瘤的发生与多种肿瘤相关基因的cpg岛的甲基化水平有关，cpg岛甲基化可以调控基因转录的效率，使基因转录失活，是细45胞内dna转录状态的重要表观遗传学标记，因此，启动子区域的cpg岛的异常甲基化也是识 别癌症的重要标志之一。针对cpg岛的定义，1987年，gardiner-garden和formmer2认为长度在200bp以上，g+c 含量大于50，并且实际cpg含量与期望cpg含量的比值(obscpgexpc g)大于0.6的区域 即为cpg岛。2002年，daiya takai 和 peter a. jones3的研究认为，长度在500bp以上，g+c50含量大于55，并且实际cpg含量与期望cpg含量的比值(obscpgexpc g)大于0.65的区域 更可能和基因的5区域有关。因此，根据研究的不同的目的，生物学研究人员可以采用不同 的参数进行cpg岛预测。dna序列中的cpg岛既可以通过生物学实验的方法识别，也可以通过计算机来识别可能 的cpg岛。由于dna序列的多样性，如果借助于生物信息和计算机工具首先识别出潜在的55cpg岛，就可以大大减少生物实验的成本和盲目性，极大地提高生物实验的时间效率和针对 性。近年来，利用计算机开发cpg岛识别工具，吸引了许多生物学领域和传统计算机领域研 究人员的兴趣，他们已经开发出许多有用的工具来帮助生物学家提高工作效率。根据这些工 具所采用的方法，我们将它们大致分为两类。第一类方法是采用传统的滑动窗口（sliding window），这些工具包括cpgis3，cpgprod4和cpgie5等。滑动窗口法从头至尾对dna60序列进行扫描（每次移动1n位），分别计算窗口内序列是cpg岛 (用“+”表示) 和不是cpg 岛 (用“-”表示) 的概率值（通常采用log-odds ratio），如果高于设定的域值，则认为窗口内 的序列是cpg岛，否则认为不是。采用滑动窗口法的不足之处在于：（1）窗口的大小很难 确定，通常窗口的大小直接影响cpg岛的长度和数量，如果窗口过大，一些长度较短的cpg 岛会被合并在一起，从而被预测成一个大的cpg岛，降低识别的准确率；（2）窗口只能向65一个方向移动，对于之前错误的判断结果，不能在后面得到发现和修改，从而影响识别的精 度；（3）系统运行时间较长。第二类方法是采用聚集（clustering）的技术，这些工具包括 cpgcluster6和cpgif7等。这类工具的基本思想是基于cpg岛定义中的密度或物理距离等生 物学属性，不断迭代地扫描序列，并将符合这些属性的相邻cpg岛合并为大的cpg岛簇。聚 集方法的提出是为了解决滑动窗口方法的不足，然而它带有如下的缺陷：（1）识别的结果70依赖于扫描序列的顺序，而不同的扫描顺序则会识别出不同的cpg岛，因此系统稳定性不足；（2）识别的敏感度较低，倾向于识别长度较长的cpg岛序列。 与上述两类方法不同，在本文，我们尝试使用隐马尔可夫模型（hmm）建立一个自动识别cpg岛的系统，并结合生物知识对预测结果进行后期修正，从而达到较为理想的识别预 测结果。隐马尔可夫模型是一种机器学习的方法，它广泛应用于序列分析和模式识别的研究75中，例如语音识别、手写识别，命名实体识别等。这类应用的一个重要特点是，数据的序列 性。这种序列可以被认为是后一状态依赖于前一状态的时间序列，即具有马尔可夫性。正是 因为考虑到dna序列（碱基对）也可以被看做是具有马尔可夫性的有序数据，我们在本文 中利用hmm模型实现了一个自动从人类dna序列中识别cpg岛的系统，即cpg-discover系 统。80本文中，我们首先介绍结合hmm模型和生物知识修正的cpg-discover系统的框架结构。 然后在欧洲分子生物学实验室（embl）发布的人类基因中cpg岛的数据集上测试本系统的 性能，并与其它常用的cpg岛识别工具（基于滑动窗口和聚集技术）进行对比实验，深入分析和讨论实验结果。最后我们总结这一工作和未来的方向。1cpg-discover 系统框架85cpg-discover 系统的框架结构主要包括三个子系统：（1）建立 hmm 模型的参数系统；（2）基于 hmm 模型建立 cpg 岛识别系统；（2）基于生物知识的后期识别修正系统。下 面我们分别详细介绍这三个子系统的工作。1.1建立 hmm 模型的参数系统这个子系统完成对 hmm 模型的求解过程，即实现参数初始化、参数逼近和参数投票这90三个步骤，建立起 hmm 模型的参数系统。在参数初始化阶段，针对 4 种碱基（a，c，t， g）分别在 cpg 岛内外的状态，我们标记此 hmm 模型包括 8 个状态，分别是 a+, c+, g+, t+, a, c, g,t，标号中“+”和“”分别表示此核苷酸在和不在 cpg 岛上。这 8 个状态之间的 转移方式如图 1 所示。这 8 个状态间的状态转移概率可以从有标识的训练序列中得到，后面 会详细叙述。本文中，我们基于欧洲分子生物学实验室（embl）核苷序列数据库进行实验95建立 hmm 模型，因此，在以下三个参数求解的过程中，我们以 id 为 bc008880 的人类 dna序列为例来进行说明。c+g+a+t+a-t- c-g- 11 -1001.1.1参数初始化图 1hmm 模型中 8 个状态之间的转移图示fig. 1 the transform of 8 states in hmm105110这个步骤完成 hmm 模型中状态转移概率矩阵（a）、发射概率矩阵（b）和状态概率 矩阵（）的初始化，作为下一步参数逼近的输入种子。三个概率矩阵的初始化步骤如下：（1）初始化状态转移概率矩阵 a：对于每个训练样本序列，我们采用基于长度为 2 的 滑动窗口的方法来初始化参数，即通过计算从第一个状态转移到第二个状态的出现频数，来 计算转移概率。（2）初始化发射概率矩阵 b：假设一个状态到一个符号的发射概率都为 0 或 1。（3）初始化状态概率矩阵 ：每个状态的初始化概率由序列中各个状态的出现频数计 算。参数完成初始化后，我们得到初始化后的各个概率矩阵。表 1 显示了人类 dna 序列bc008880 经过初始化之后的三个概率矩阵的分布结果。1.1.2参数逼近在第二步中，我们采用 baum-welch 算法来重新估计基于训练序列的概率矩阵。上一步 的三个初始化概率矩阵作为 baum-welch 算法的输入，输出期望的转移和发射概率矩阵，再115将得到的期望概率代替旧的参数，这个过程不断迭代进行直至达到收敛。在每次迭代过程中，模型概率的准确度都得到提高，直到一个极限概率，整个迭代过程最终收敛于一个局部最优 解。参数估计后，对每个训练序列，我们都得到三个逼近后的概率矩阵。表 2 显示了人类 dna 序列 bc008880 经过 baum-welch 近似之后的三个概率矩阵的分布结果。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 经过初始化之后的三个概率矩阵的结果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 经过 baum-welch 逼近之后的三个概率矩阵的结果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参数投票至此，对每个 dna 训练序列，我们都得到三个逼近后的概率矩阵。为了整合各个不同 长度的训练序列的概率矩阵，我们采取投票（vote）平均技术，就是说，对每个概率矩阵， 根据每个训练序列的长度对其进行归一化（normalization），每个 dna 训练序列的权重计算公式为： w = ni，其中，n 表示第 i 个训练序列中核苷酸的数目，n 表示训练集中的核i n i135苷酸总数。采用投票平均后，我们得到最终的 hmm 模型的参数。表 3 显示了使用 1400 条 “bc”打头的人类 dna 序列作为 hmm 系统的训练数据，经过加权平均（投票）之后的 hmm 模型的三个概率矩阵的分布结果（关于数据集见后面 2.1 节）。1401451.2基于 hmm 模型建立 cpg 岛识别系统经过求解和加权平均之后得到的 hmm 模型，可以通过 viterbi 解码算法来寻找能够产 生给定符号序列的最大似然状态序列，即对每个测试 dna 序列中的每个碱基对（atcg）， viterbi 解码算法能够输出每个碱基对对应的状态为”+”或者”-”（其中“+”标识碱基对在 cpg 岛内，”-”表示碱基对在 cpg 岛之外），并使得标识后的 dna 序列输出状态序列为最优解。表 3 使用 1400 条以“bc”开头的人类 dna 序列作为训练数据集，经过加权平均（投票）之后的 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基于生物知识的后期识别修正系统由viterbi 算法求解出的dna序列输出状态是基于训练模型计算而得，系统输出的状态 序列并不一定满足生物学上对cpg岛的定义。为了进一步提高结果的精度，我们结合cpg岛 的生物定义，如cpg岛中c和g的概率、相邻cpg片段间距以及cpg片段的长度阈值等，对 hmm模型识别dna序列产生的误差进行后期的知识修正。我们对cpg岛作如下的限定：至 少包含140个碱基对的长度，区域内gc所占含量超过60%，且cpg的观察值预测值比例必 须高于0.6（根据研究的不同目的，可以在软件系统中自行修改这些限定的值）。这个后期 识别修正系统包含下面三个模块：（1）相邻cpg岛的合并模块：由于计算误差，系统经常把cpg岛中少量的核苷酸（通 常长度在1bp到15bp之间）错误地识别为非cpg岛上的核苷酸。为了解决这个问题，我们对 相邻cpg岛设定最小距离限定。经过多次实验，这个最小距离的阈值设定为20，即，如果识 别出的两个相邻cpg岛间的距离小于20，系统就修改这些间隔中的非cpg岛核苷酸状态，使 之状态转变为cpg岛上核苷酸，即合并这两个相邻cpg岛为一条长cpg岛。合并后的cpg岛 需要进一步经过后面两个模块的修正。（2）密度检测过滤模块：生物学上，cpg岛通常含有较高的g/c（或c/g）含量(密度高 于60%)，如果前一步的错误合并或者系统识别的误差，可能产生不满足c/g含量要求的cpg 岛序列。因此，对识别出的不能满足最低密度要求的cpg岛序列，即区域内gc所占含量至 少达到60%，这个密度检测过滤模块将改变这些cpg岛中所有核苷酸状态为非cpg岛状态。（3）长度检测过滤模块：除了（2）中的密度约束模块外，这个系统还进行了长度约束， 即根据cpg岛的长度定义，去除长度少于140个碱基对的cpg岛序列，即，对识别出的长度175180小于140bp的cpg岛将被改变为非cpg岛状态。经过以上三个后期修正模块之后的输出序列状态，才成为cpg-discover系统的最终输入结果。表4列出了人类dna序列bc009465在修正前后的对比结果，其中491255分别表示cpg岛序列的起始位置为第49个碱基对和结束位置为第1255个碱基对，这6个评估指标的含义在后面2.2节中有详述。从表4可以看到，基于生物知识修正后的识别结果除了correlation coefficient指标没有变化外，其它的指标都有提高，特别是specificity，从原来的33.76%增至了51.68%（提高60.15%）。原来预测位置在1284和2961509的两条cpg岛因为相隔只有12个碱基对（bp），小于最小间隔阈值20bp，因而被修正合并为一条长的cpg岛；而原来预测位置在15841695的cpg岛序列长度只有120bp，在后处理模块中由于长度过短，不满足 生物学上cpg岛的长度要求，因而被过滤除去。表 4 人类 dna 序列 bc009465 基于生物知识的修正前后识别结果的对比tab. 4 the comparison of domain knowledge revision (bc009456)评测指标未修正的结果修正后的结果真实cpg岛的始末端位置491255491255系统识别的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 系统的性能做一个准确而全面的分析和比较，我们从欧洲分子生 物学实验室(embl)的核苷序列数据库（http:/www.ebi.ac.uk/embl/）中下载了关于人类基因 的 cpg 岛数据库 emb173hum（ftp:/ftp.ebi.ac.uk/pub/databases/cpgisle），并在这个数据集上 测试了本系统的性能。这个数据库是由欧洲生物信息中心负责维护，与日本基因数据库（ddbj）和美国国立生物技术信息中心（ncbi）相互合作，他们之间每天都要交换最新的 核苷序列数据。此外，为了比较 cpg-discover 系统与其它常用的 cpg 岛识别工具（基于滑 动窗口和聚集技术）的性能，我们还在这个数据集上选择部分 dna 序列进行了对比实验。2.1数据来源emb173hum 数据库包含了当前 embl 数据库发布的已标定 cpg 岛位置的人类 dna 序 列，共有 233,004 条 dna 序列和 142,325 个 cpg 岛，其中只有 61,051 条 dna 序列中包含 有一个或者若干个 cpg 岛，平均每条 dna 序列含有 2.33 个 cpg 岛。我们下载了这个数据 库中以“ab”标记开头的 1,955 条 dna 序列和以“bc”标记开头的 6,124 条 dna 序列，除去 数据库中的空序列和有错误信息的序列之后，我们得到近 8,000 条 dna 序列。其中，以“ab” 标记开头的 dna 序列中平均每条 dna 序列含有 1.70 个 cpg 岛，以“bc”标记开头的 dna 序列中平均每条 dna 序列含有 1.34 个 cpg 岛。为了使实验结果真实合理有效，我们从以 “ab”和“bc”标记开头的 dna 序列中分别选择了 1400 条训练数据和 200 条测试数据（共有2800 条训练数据和 400 条测试数据），选择的原则是这个数据子集里面的 cpg 岛的分布分 别与原来“ab”和“bc”标记开头的 dna 序列中 cpg 岛的分布相同，目的是使我们的实验数 据尽可能符合原来 cpg 岛在 dna 序列中的真实数据分布概率。2052.2评价手段为了对系统进行准确客观的评价，对每一条测试 dna 序列，我们使用以下 6 种评估指 标对 cpg-discover 系统识别出的 cpg 岛位置和真实数据库中标定的 cpg 岛位置进行对比。 这 6 种评估指标都在生物信息学领域有广泛的应用，它们的表示公式如下：(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)其中各个参数的定义如下：lntp：是 cpg 岛内的碱基且被正确的预测为 cpg 岛内碱基的碱基个数。lnfn: 是 cpg 岛内的碱基且被错误的预测为非 cpg 岛内碱基的碱基个数。lnfp: 不是 cpg 岛内的碱基且被错误的预测为 cpg 岛内碱基的碱基个数。lntn: 不是 cpg 岛内的碱基被正确的预测为非 cpg 岛内碱基的碱基个数。2.3实验结果与分析在下面的实验过程中，根据 cpg 岛的定义，参考其它软件给出的可选参数，我们选择 了如下的参数作为我们的筛选参数：长度在 140bp 以上，g+c 含量大于 60。为了研究训练数据集的大小对系统测试结果的影响，我们分别把“ab”和“bc”开头的dna 序列训练数据集分成 7 组不同大小的训练数据集，即分别取200, 400, 600, 800, 1000,1200, 1400条训练数据来建立 hmm 系统，然后使用预先留出的 200 条测试 dna 序列来评 估系统的性能参考。图 2 和图 3 分别显示对以“ab”和“bc”开头的不同训练数据集对 6 个系 统性能评估指标的结果。230235图 2 ab 序列的测试结果fig.2 the test result of nda sequence(ab)图 3 bc 序列的测试结果fig.3 the test result of nda sequence(bc)从图 2 和图 3 显示的结果可以看出，对每一种性能评价指标来说，随着训练数据数量的增加，系统的性能指标的均值总体是增加趋势，而当训练数据集包含大约 1000 条数据的时240245250255260265270275候，系统的各个性能指标趋于稳定。这个结果表明，当系统在接受 1000 条训练数据之后，即可以输出相对稳定的测试结果。此外，从图 2 和图 3 的对比结果还可以看出，“bc”开头的 dna 序列数据测试结果相对 于“ab”开头的数据的测试结果要更加稳定。通过对 dna 序列数据的分析，我们认为可能的 原因如下：l首先，从训练数据来看，“ab”开头的数据中，每条 dna 序列中的 cpg 岛的数量波动 范围集中在为 17，而以“bc”开头的数据中，每条 dna 序列中的 cpg 岛的数量波动范 围集中在 14。l其次，“ab”开头的 dna 序列平均长度为 5924bp, “bc”开头的 dna 序列平均长度为1985bp。因此，相对“ab”开头的 dna 序列，“bc”开头的数据本身长度短，cpg 岛数量波动小，比“ab”开头的数据更加规范和集中，因而系统的波动性也更小，表现出较为平稳的结果趋势。 从上述实验结果可以看出，cpg-discover 系统的 sensitivity 性能很高，最高达到了 97%，即使选择不同训练数据集大小，系统的平均 sensitivity 值也超过 90%，尤其是“bc”开头的较短 dna 序列数据。从 sensitivity 的定义来看，cpg 岛内 90%的碱基都能被正确识别标定 出来，也证明 cpg-discover 系统能达到较为满意的从人类 dna 序列中自动识别 cpg 岛的 作用。此外，在最好的情况下，cpg-discover 系统的 accuracy 可以达到近 80%的正确率， 而且绝大部分的 accuracy 也超过了 50%。2.4与其它工具的对比为了进一步全面地评估系统的性能，我们随机地从测试样本中选择若干条 dna 序列， 将 cpg-discover 系统与前面提到过的基于滑动窗口和聚集技术的工具进行了对比实验。我 们选取了目前较为有名的 cpg 岛识别工具或者系统，例如基于滑动窗口技术的 cpgprod4 和 cpgie5，基于聚集技术的工具 cpgcluster6和 cpgif7来进行比较，这些工具软件涵 盖了当前广泛采用的算法，并且在 cpg 岛挖掘方面公认有较好的性能。表 5 列出了部分 dna 序列分别使用这 5 个系统工具预测的 cpg 岛位置和真实 cpg 岛位置的对比结果此外，为了更全面比较 cpg-discover 与其它同类软件的结果，我们对其他软件的预测 结果使用 2.2 节中的 5 种评价指标进行衡量比较。表 6 以 ab046787 和 bc011652 这两条 dna 序列为例，列出了各个系统对这两条 dna 序列进行识别的 5 个性能评价指标的结果。从表 5 和表 6 的结果可以看出，当相邻的 cpg 岛的间隔较小的时，以上五个软件都不 能准确确定 cpg 岛的起始边界位置，cpgif 和 cpgpord 倾向于将间隔较小的两个 cpg 岛 合并成一个 cpg 岛，而 cpg-discover 则不能将其间隔准确定位（如 bc008956、bc017346）， 而 cpgie 和 cpgcluster 出现了较多无法预测出 cpg 岛的情况。cpgcluster 较其他软件而言有很好的 specificity 和 positive predictive value，这表明 cpgcluster 识别出来的 cpg 岛极可能是真正的 cpg 岛。但是 cpgcluster 的 sensitivity 却与 其他软件有很大的差距，这说明对于很多是 cpg 岛的区域，cpgcluster 很可能不能识别出 来而将其判成了非 cpg 岛的区域。原因是 cpgcluster 的识别结果与输入序列有关，如果待 测序列的非岛区域较大，则会导致其 cpg 岛区域的 p-value 偏低，因此造成了 cpgcluster 将 cpg 岛的区域判断成非 cpg 岛的区域的情况。cpgprod 的 specificity 相对于其他软件而 言比较低。而 cpgie、cpgif 和 cpg-discover 的各个评价参数都相对较好。这个结果非常 有意思，因为这三个工具正是分别采用了滑动窗口，聚集技术和 hmm 模型，这也解释了为 什么目前采用不同技术的 cpg 岛识别系统同时都在得到广泛应用中。280表 5 cpg-discover 系统与其它四个软件工具的比较结果。选用参数长度在 140bp 以上，g+c 含量大于60 ，并且实际 cpg 含量与期望 cpg 含量的比值(obscpgexpc g)大于 06 ，其中 cpgcluster 工具的 参数取默认值为 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序列标号bc011973bc008956bc014043ab037841ab046787ab007883长度 (bp)161812333079451659906246真实 cpg 岛的起始和结束位置49507483764961024804100414751989239526104839847264190408cpgcluster135191813900611725217576639cpgif125327104979210361324271522525552464811cpgprod1883113087352628180317492031200cpgie16071522140556485870cpg-discover17051507577938130242873179696810683038137512551780长度 (bp)161812333079451659906246285表 6 cpg-discover 与其它同类软件使用 5 个性能评价指标的对比结果tab. 6 the performance evaluation of cpg-discover and other tools.ab046787bc011652系统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 系统在预测 cpg 岛的方面有很好的表现。首先，在 specificity 和 sensitivity 方面都比较优秀，这表明 cpg-discover 系统既可以识别出绝 大部分 cpg 岛，同时又有很好的识别精度，可以比较精确的定位 cpg 岛的位置。其次， cpg-discover 系统是基于训练数据的自动识别系统，也就是说，给定标记的 cpg 岛训练数 据，系统一经训练建立，可以应用到各种 cpg 岛的自动识别中，大大减少人工干预和劳动。3结论本文中，我们实现了一个基于隐马尔科夫模型（hmm ）的 cpg 岛自动识别系 统cpg-discover，在欧洲分子生物学实验室(embl)的核苷序列数据库发布的人类 dna 序列数据集上