基于隐马尔可夫模型的蛋白质序列筛选算法

张毅，梅挺

（成都医学院 人文信息管理学院，四川 成都 610083）

近年来，随着人们对蛋白质测序工作的快速发展，蛋白质数据库中的序列数量呈现指数级的增长速度，在这存储有海量的蛋白质数据库中，存在着大量的冗余蛋白质序列。虽然目前对冗余蛋白质序列尚未有很明确和统一的定义，但是普遍认为，在蛋白质数据库中，如果两条蛋白质序列的具有非常高的相似度，尤其是在整个序列中控制蛋白质功能的特征序列具有很强的相似性时，则认为这两条蛋白质序列是互为冗余的序列。造成这种现象的原因很多，一个典型的原因即针对某一同源的蛋白质序列进行的测序，并将测量的结果存入数据库中。

由于蛋白质数据库在医学研究、物种研究等方面发挥着非常重要的作用，利用蛋白质数据库中的信息，有助于人们发现新物种，寻找物种之间的生物关系，研究针对某些特殊病毒的抗生药物等。然而，如果在蛋白质数据库中存在大量的冗余序列，则可能导致对这些蛋白质数据分析的误差加大。比如在某一蛋白质序列簇中，如果冗余序列过多，可能会夸大这一序列簇的某些功能特征，从而对蛋白质序列间相互关系的研究产生误导[1]。

目前，针对蛋白质数据库冗余的问题，国内外有不少相关学者开展了研究，比较有代表性的有Hoblhm和Sander提出的CD-HIT去冗余算法，该算法的设计思想是在每个蛋白质序列簇中选取一个序列作为特征序列，然后再将该簇中的其余序列进行冗余检测，如果某条序列与该序列的匹配程度超过某一阈值，则该序列被视为冗余序列。该算法是一种非常经典的蛋白质去冗余算法，目前也有很多算法是基于这一思想进行变形实现。此外还有在2000年由Yona和Linial共同提出了蛋白质序列聚类Protomap算法，Enright和Ouzounis提出了蛋白质序列分级聚类Generage算法，2004年由Kawaji和Takenaka提出了基于图论的蛋白质序列分类算法。总体而言，这些蛋白质分类算法在分类的精确性和分类准确性两方面还有待提高。

1 模型的建立与参数估计

1.1 模型的建立

从蛋白质序列数据库中选择一个蛋白质序列作为研究对象，记为D=D1D2…Dn，该研究对象也被称为观测对象。一个数据库中的蛋白质序列也被视为隐马尔可夫链的初始状态的序列分布π，每一条蛋白质序列作为一个状态的迁移中的中间状态，观测结果是初始状态经过一定的迁移变化，且和一些随机过程共同作用的状态结果。状态的迁移概率P1和观测得到某一特定状态结果的概率P2在分析之前需要进行确立，确立之后为隐马尔可夫[2]可表示为 φ=（π，P1，P2）。

为了能够对蛋白质数据库中所有的序列进行匹配筛选，因此需要将蛋白质数据库中的序列出现的概率进行抽象和描述。然后才能在此基础上设定筛选的规则，并按照筛选规则确定的权值，筛选出最具有代表性的蛋白质序列。

假设在给定的隐马尔可夫模型中，已经有观测到某一序列的部分值为：{D1D2…Dk}，且在时刻k，序列的状态为 w1的概率[3]为：

由于该概率的初始值可根据蛋白质数据库中的常量计算得到，因此，该计算式可以采用递归的方式计算得到[4]：

该式经过递归化简后可得：

1.2 模型参数的估计

模型参数建立起来之后，将考虑如何选取合适的参数，使得对于一个给定的隐马尔可夫模型φ=（π，P1，P2），如何得到其特定子序列的概率最大值，即P（D|φ）为最大。

假定对一个给定的隐马尔可夫模型 φ=（π，P1，P2）和得到的观测序列 D={D1，D2，…，Dn}，在时刻 k 的状态为 wi，时刻 k+1 的状态为 wj，满足这种特征的概率[5]记为 P（i，j），则有：

所以，通过上面的概率计算式，可以得到在模型确定的条件下，对于一个给定的观测序列（即某一条蛋白质序列）下，k时刻的状态为wi的概率[6]：

其中 N≥i≥1， N≥d≥1。

2 同源蛋白质序列筛选过程

蛋白质序列筛选过程中最为关键的是找出每个蛋白质序列的关键信息，这些信息直接决定了该蛋白质序列的主要功能和特征。因此定位和筛选出每条蛋白质序列中的关键信息是进行同源蛋白质序列筛选的重要前提。

蛋白质数据库中每条序列记为D=D1D2-Dm-1Dm，隐马尔可夫链阶梯步长记为spl，在一条完整的序列中定义一个子片断记为d=d1d2-dr-1dr。则子片断d出现的概率[8]为：

由此，可定义若存在子片断d条件下，各个序列值的概率[9]：

按照这种条件概率计算式，可以进一步得到子片断d的详细计算式[10]：

再定义整个蛋白质序列中，每一种可能的子片断xd出现的概率[11]：

根据如上定义，可以得到任意一个子片断的重要度描述参数Weight[12]。

确定了任意一个子片断的权值后，再按照序列的模式匹配，定位序列之间的匹配位置和匹配程度。若两个对比的蛋白质序列长度不等，记 D1中有 d[1：x]特征序列，记 D2中有 d[1：y]特征序列。通过调整特征序列的长度，记录两个蛋白质序列中的最大匹配权值，该权值即为两个蛋白质序列的匹配程度描述参数。对于两个特征序列d[1：x]和d[1：y]的最大匹配权值计算式为[7]：

如此递归计算，最后得到最终结果。

3 测试结果

采用文中设计的基于隐马尔可夫模型的蛋白质序列筛选算法SWISS-PROT对蛋白质数据库中的数据进行了分类筛选测试，通过本文设计的算法，对蛋白质数据库中的蛋白质序列特征进行提取和匹配，当匹配到两条蛋白质序列的关键信息是一致的，则将这两条蛋白质序列标记为同源蛋白质序列，同时将得到的筛选结果与目前公认的同源序列结果进行对比，得出筛选的正确率。测试结果如表1所示。

表1 基于隐马尔可夫模型的蛋白质序列筛选算法测试结果Tab.1 Hidden Markov Model-based protein sequence selection algorithm results

测试结果表明，文中所设计的蛋白质筛选算法对蛋白质数据库中筛选出了大量的同源蛋白质序列，根据与目前公认的同源蛋白质序列进行结果对比表明，在筛选过程中筛选出了蛋白质数据库中的绝大多数的同源蛋白质序列，且筛选结果的正确率达到了95%以上，从筛选的精度而言，本设计的算法具有较高的筛选正确率。

4 结束语

从蛋白质数据库中对蛋白质序列进行分类和筛选有着非常现实的意义，尤其是面对指数级增长的蛋白质数据库，只有通过对蛋白质数据库进行去冗余的处理，筛选出真正代表每个蛋白质簇的特征序列，才能建立更有实用价值的蛋白质数据库。目前这一问题也是全球蛋白质序列研究的一个热点问题，本文通过引入隐马尔可夫模型对蛋白质序列进行分析筛选，探索了蛋白质数据库特征序列筛选的新方法，已进行的测试结果也表明该方面在筛选的精度上可以达到95%以上。

[1]张成岗，欧阳曙光，张绍文，等.基于PC/Linux的核酸序列分析系统的构建及其应用 [J].生物化学与生物物理进展，2001（2）：263-266.

ZHANG Cheng-gang， OU YANG Shu-guang， ZHANG Shaowen， et al.Based PC/Linux system， Construction and application of the nucleic acid sequence analysis system based on PC/Linux[J].Biochemistry and Biophysics，2001（2）：263-266.

[2]陈英，彭心昭，朴英杰.自噬基因APG5基因结构的生物信息学分析[J].遗传学报，2001，28（11）：1077-1084.

CHEN Ying， PEN Xin-zhao， PU Ying-jie.Bioinformatics analysis of autophagy gene APG5 gene structure[J].Genetics，2001,28(11)：1077-1084.

[3]齐建勋，肖奕.基于小波方法的蛋白质非规则二级结构预测[J].科学通报，2002（6）：425-430.

QI Jian-xun，XIAO Yi.Non-wavelet-based method of protein secondary structure prediction rules[J].Chinese Science Bulletin，2002（6）：425-430.

[4]任力锋，张波，刘辉.蛋白质序列信息的提取与蛋白质结构预测[J].北京生物医学工程，2005（3）：237-238.

REN Li-feng， ZHANG Bo， LIU Hui.Protein sequence information extraction and protein structure prediction[J].Beijing Biomedical Engineering，2005（3）：237-238.

[5]霍红卫，肖智伟.基于最大权值路径算法的DNA多序列比对方法[J].软件学报，2007，18（2）：185-195.

HUO Hong-wei，XIAO Zhi-wei.A multiplealignment approachforDNA sequencesbasedonthemaximum weighted path algorithms[J].Journal of Software，2007，18（2）：185-195.

[6]邹权，郭茂祖，王晓凯，等.基于关键字树的DNA多序列星比对算法[J].电子学报，2009，37（8）：1764-1850.

ZOU Quan， GUO Mao-zu， WANG Xiao-kai， etal.Keyword-based tree of the DNA sequence star more than the algorithm[J].Electronics Technology，2009，37（8）：1764-1850.

[7]王艳春，何东健.神经网络在蛋白质二级结构预测中的应用[J].安徽农业科学，2006（16）：4172-4174.

WANG Yan-chun；HE Dong-jian，Neural network in protein secondary structure prediction in two[J].Anhui Agricultural Sciences，2006（16）：4172-4174.

[8]阮晓钢，孙海军.编码方式对蛋白质二级结构预测精度的影响[J].北京工业大学学报，2005，31（3）：227-235.

RUAN Xiao-gang，SUN Hai-jun.Researchon encode influencing protein secondary structure prediction[J].Journal of Beijing University of Technology，2005，31（3）：227-235.

[9]刘帅，马志强，刘清雪，等.基于自适应免疫遗传算法的多序列比对[J].信息技术，2007（2）：15-17，111.

LIu Shuai， MA Zhi-qiang， LIU Qing-xue， et al.Adaptive immune genetic algorithm based on multiplesequence alignment[J].Information Technology，2007（2）：15-17，111.

[10]郭卫斌，施保昌，王能超.多重生物序列对准及其算法综述[J].高技术通讯，2001，11（6）：96-102.

GUO Wei-bin， SHI Bao-chang， WANG Neng-chao.Multiple biological sequence alignment and its algorithm[J].High Technology，2001，11（6）：96-102.

[11]关敏，辜华良，常雅萍，等.DNA核苷酸碱基序列分析软件的编写和应用[J]，白求恩医科大学学报，2001，27（5）：467-469.

GUAN Min， GU Hua-liang， CHANG Ya-ping， et al.DNA nucleotide base sequence analysis software and application[J].BethuneUniversityofMedicalSciences，2001，27（5）：467-469.

[12]杜世平.隐马尔可夫模型在生物信息学中的应用[J].大学数学，2004，20（5）：24-29.

DU Shi-ping. HMM in bioinformatics applications[J].University Mathematics，2004，20（5）：24-29.