基于隐马尔科夫模型的基因预测算法
第34卷 第4期
2008年11月
大连海事大学学报
JournalofDalianMaritimeUniversity
Vol.34 No.4Nov., 2008
文章编号:100627736(2008)0420041204
基于隐马尔科夫模型的基因预测算法
马宝山,朱义胜
(大连海事大学信息科学技术学院,辽宁大连 116026)
Ξ
摘要:为更好地从DNA序列中识别蛋白质编码区,提高计算效率,将隐马尔科夫模型与前向算法相结合,提出一种生物基因的预测算法.理论证明,该隐马尔科夫模型与经过EM算法优化后的模型具有相同的参数,能够降低计算量,提高计算效率.利用该方法对DNA序列F56F11.4a的外显子进行识别的仿真结果表明,该算法是有效的.
关键词:基因识别;隐马尔科夫模型;前向算法;EM算法中图分类号:TP18 文献标志码:A
主要的课题之一
YOON等
[2]
[1]
.
建立了一种csHMM,可以看作是对
隐马尔科夫模型的一种推广,它所包含的发生概率和过渡概率可随“上下文”而改变,这对研究某些具有相关特性的模型非常有效,RNA二级结构[2]
[YOON又提出了基
[]
,改进的算法.罗泽举等
[4]
Gene2predictiononMA,sheng
(CollegeofInfScienceandTechnology,DalianMaritimeUniversity,Dalian116026,China)
提出一种基于隐马
[5]
.实验表明,其性能优于支持向量机方法.朱红梅等提出一种延时隐马尔科夫模型.实验表明,该模型优于标准隐马尔科夫模型.通常,用于基因预测的隐马尔科夫模型一般需先对参数进行优化后再进行预测,且优化算法需要大量的运算.为此,本文建立基因外显子的隐马尔科夫模型,并结合前向算法,研究一种用于蛋白质编码区识别的新方法,提高了计算效率.
Abstract:PredictionalgorithmofbiologicalgeneswasdevelopedcombiningtheforwardalgorithmwithhiddenMarkovmodel(HMM)toidentifyproteincodingregionsfromDNAsequencesandimprovingcalculationefficiency.TheoreticalproofshowsthattheHMMhasthesameparameterswiththemodelopti2mizedbytheexpectationmaximization(EM)algorithm,whichcanreducethecalculationcapacityandimprovethecalculationefficiency.SimulationsontheDNAsequenceF56F11.4aalsoshowthattheproposedalgorithmiseffectiveandefficient.Keywords:geneidentification;hiddenMarkovmodel;forward
algorithm;expectationmaximization(EM)algo2rithm
1 隐马尔科夫模型
隐马尔科夫模型是由马尔科夫链扩展而来的一种随机过程过程.
1.1 隐马尔科夫模型的定义
[6]
,可以理解为一个双随机过程,分别为
系统状态变化的随机过程和由状态决定输出的随机
隐马尔科夫模型定义为λ={S,O,A,B,π},
0 引 言
真核生物基因的编码序列在DNA分子上不是
连续排列的,而是被不编码的序列隔开,其中外显子为编码序列,内含子不编码.因此,从真核生物基因组DNA序列中预测蛋白质编码区,是生物信息学最
Ξ收稿日期:2008203207.
其中:状态集合S={s1,s2,s3,…,sN},N为状态总数;观测事件集合O={o1,o2,o3,…,oM},M为观测事件总数;状态转移概率矩阵A={aij},aij表示由状态si转移到状态sj的转移概率,记为aij=
P(qt+1=sj|qt=si),1≤i≤N,1≤j≤N,qt
基金项目:重大基础研究前期研究专项资助项目(2005CCA02200),国家自然科学基金资助项目(60671061).作者简介:马宝山(1977-),男,吉林德惠人,博士研究生,讲师,E2mail:[email protected];
朱义胜(1945-),男,山东德平人,教授,博士生导师,E2mail:[email protected].
大连海事大学学报 第34卷 42
表示在t时刻模型所处的状态,并且有aij≥0,
Nj=1
～
T
t
T
t
bjk=[
k
∑a
ij
=1;在状态sj时产生观测事件vk的概率矩阵
t=1
o=v
γ(j)]Πγ(j)][∑∑
t=1
k
(3)
B={bjk},bjk=P(vk|qt=sj),1≤j≤N,1≤k≤M,vk∈O;初始状态分布矩阵π={πi},其中
N
3 前向算法
定义前向变量:α=P(o1o2…ot,qt=t(i))表示模型λ下,在t时刻观测事件为ot、si|λ状态
πi=P(q1=si),1≤i≤N,并且有πi≥0,=1.
1.2 隐马尔科夫模型
i=1
π∑
i
为si的概率.则在t+1时刻,
α)t+1(j)=P(o1o2…ot+1,qt+1=sj|λ
N
图1为基因外显子的一种三状态隐马尔科夫模
型.模型中3个状态分别代表三联密码子的3个状态,其中每个状态只能单向运动,即1至2,2至3,3又回到1,不能反向运动或发生跳跃.因此,该模型的转移概率为a12=a23=a31=1,其余转移概率为
0.
=[
i=1
α(i)a∑
i
ij
]bj(ot+1)
(4)
则前向算法为:
(1)初始:α1(i)=πibi(o1),1≤i≤N(2)递归:
N
(j(i)t
N
]bj(ot+1)
-1,1≤j≤N
图1)= (3)终止:P(O|λ
i=1
α(i)∑
T
2 EM算法
EM算法可以重新估计隐马尔科夫模型的参数
A、B和π,经常用来对隐马尔科夫模型进行训练,以
4 基于隐马尔科夫模型的EM算法
由图1所示隐马尔科夫模型,可以得到模型的概率转移矩阵A={aij},a12=a23=a31=1,其余转移概率为0;参数初始状态π={π0}={1,0,0}.
将初始条件代入EM算法中得到
T-1
T-1
t
增强隐马尔科夫模型的预测性能
[7]
.EM算法利用最
)优化的方法,使模型生成观测序列的概率P(O|λ
最大.算法的思想是先得到一组参数值,然后通过迭代计算,直到算法收敛至一个局部最优值.
定义t时刻状态为si以及t+1时刻状态为sj的概率为
ξλ)t(i,j)=P(qt=si,qt+1=sj|O,
=
)P(qsqsOλ
)P(O|λ
=
αij=
～
T-1
ξ(i,j)∑
t=1
=
∑P(q
t
=si,qt+1=sj)
t
T-1
γ(i)∑
t
t=1
∑P(q
=si)
()()
P(q1=si)+P(q2=si)+…+P(qT-1=si)
为方便分析,假设时间长度为3的倍数,由隐马尔科夫模型及初始条件可知P(qt=si,qt+1=sj)与P(qt=si)在t时刻总是同时为0或1,所以上式
(1)
αi(i)aijbj(ot)βt(j)=
)P(O|λ
αt(i)aijbj(ot+1)βt+1(j)=NN
αt(i)aijbj(ot+1)βt+1(j)∑∑
i=1j=1
简化为
αij=
～
其中:αt(i)和βt+1(j)分别为前向变量和后向变量,
t时刻状态为si的概率为
N
(TΠ3)×P(qt=si,qt=sj)
=
(TΠ3)×P(qt=si)
P(qt+1=sj|qt=si)=aij
(5)
同理
t
γt(i)=
则EM算法为
πi=γαij=[1(i),
～
～
j=1
ξ(i,j)∑
T-1
t
t
(2)
～
TT
t
T-1
bjk=
t=1o=v
k
γ(j)∑
k
T
=
t
t=1
o=v
k
∑P(q
k
t
=sj)
T
t=1
ξ(i,j)]Πγ(i)][∑∑
t=1
t=1
γ(j)∑
t=1
∑P(q
t
=sj)
第4期 马宝山,等:基于隐马尔科夫模型的基因预测算法 43
T
=
t=1
∑
P(qt=sj,ok=vk)
T
5.2 仿真步骤
(1)根据编码区密码子使用的偏好性,建立三状
t=1
∑P(q
t
=sj)
态隐马尔科夫模型,如图1所示.模型包含3个状态1、2和3,分别代表外显子的三个状态.每个状态对
P(q1=sj,ok=vk)+…+P(qT=sj,ok=vk)=
P(q1=sj)+P(q2=sj)+…+P(qT=sj)
==
(T)P(qsov)
(TΠ3)×P(qt=sj)
P(qt=sj,ot=vk)
P(qt=sj)
应的观测事件为4种,分别为四种碱基A、T、C和G.起始状态π={π0}={1,0,0}.状态转移矩阵A={aij},a12=a23=a31=1,i=1,2,3,j=1,2,3.
每个状态下取所有观测事件的概率分布B={bjk}.取已标注的DNA序列F56F11.4a中第1个外显子作
(6)
=P(ot=vk|qt=sj)=bjk 训练序列,长度为111,其中每个3联密码子中的3个碱基分别对应状态1、2和3,统计出每个状态对应的4种碱基的个数,除以该状态出现的总数得到每个状态对应观测事件的条件概率.假设状态sj=1,2,3,碱基k=A,T,C,G,bjk=nj(k)ΠNj,其中
nj)sjk出现的总次数,Nj
实际上,即使T不为3的倍数,上面的结论也是成立
的.因此,基于文中隐马尔科夫模型的EM算法简化为
～～～π=π(7)0 aij=aij bjk=bjk 可见,经过EM算法优化后的系统参数与图1所示隐马尔科夫模型的参数是一致的.文中状态总
数N=3,观测事件总数M=4.由式(7)～
sj,Nj=NΠ3,N为外显子.
(2)使用EM算法对隐马尔科夫模型的参数进行优化,其中A和π初值同上,B的初始值可以是任意值,但要为正数,并且每行元素之和为1.经过EM算法优化计算后,模型的参数值B与隐马尔科夫模型的参数B相同,即
0.2644
B=B=
～
～
bjkN×T=12×TEM～
～
都需要重新计算aijjk.由式(3)可知,计算aij需要N×N×[(T-1)+N(T-1)]次运算,计算bjk需要N×M×(N×T+N×T)次运算,优化一次模型参数共需要计算N×N×(N+1)×(T-1)+2×N×N×M×T次.代入数值为36(T-1)+72T=108T-36.故文中提出算法的计算量比EM算法优化过程的计算量大大降低了.
～
0.12620.25370.3235
0.25910.28990.1221
0.
0.11810.0.33830.2671
训练用的数据是已标注的DNA序列F56F11.4a的第1个外显子,与序列的其他4个外显子同属一
5 仿真及分析
5.1 仿真数据
个DNA序列,基因结构上具有很高的相似性,因此使用同源序列进行训练后,会使得算法对基因的预测更加准确.
(3)根据获得的模型参数,结合前向算法计算预测序列的生成概率.概率大于设定阈值的,可认为是编码序列,否则认为是非编码序列.5.3 仿真结果
本文计算使用的数据是从美国国家生物信息中心网站www.ncbi.com下载的基因序列F56F11.4a(在染色体中排列的碱基号从7949～14625,编号:AF099922),该序列共有5个外显子,如表1所示.
表1 基因序列F56F11.4a的外显子预测结构
Exon#12345
RelativetoItselfStart[***********]5
End[***********]05
RelativetoF56F11Start[***********]14275
End[***********]14625
通过EM算法优化后的模型参数与隐马尔科夫模型参数是相同的,并且预测的外显子结果也是相同的,如图2所示.图中5个黑色方框区域为预测的外显子位置.对照表1中的外显子实际位置,预测的外显子位置是正确的.
6 结 语
本文基于隐马尔科夫模型与前向算法,提出了一种基因外显子的预测算法.通过EM算法优化后
大连海事大学学报 第34卷 44
quenceswithcomplexsymbolcorrelationsusingprofilecont2
ext2sensitiveHMMSandpre2screeningfilters[C].ΠΠICASSP2007,Hawaii,USA:IEEEPress,2007,1:3452348.[4]罗泽举,李艳会,宋丽红,等.基于隐马尔可夫模型的DNA序列识别[J].华南理工大学学报:自然科学版,2007,35(8):1232126.
LUOZe2ju,LIYan2hui,SONGLi2hong,etal.Recogni2tionofDNAsequencesbasedonhiddenMarkovmodels[J].
图2 预测的基因外显子结果JournalofSouthChinaUniversityofTechnology:NaturalScienceEdition,2007,35(8):1232126.(inChinese)[5]朱红梅,王家 ,赵燕南,等.延时HMM在基因剪接供体
的模型参数与该隐马尔科夫模型参数相同,但该算法的计算量大大降低,计算效率明显提高.对已标注的DNA序列的预测结果表明,该方法是有效的.参考文献(References):
[1]张春霆.生物信息学的现状与展望[J].中国科技研究与
位点识别中的应用[J].计算机工程,2007,27(5):123.
ZHUHong2mei,WANGJia2xin,ZHAOYan2nan,etal.Applicationoftime2delayHMMingenesplicedonorsitesprediction[J].ComputerEngineering,2007,27(5):123.(inChinese)
[6]LRintroductiontohidden
].,SpeechandSignal,1986,3(1):4216.
[M,OPPENHEIMAV,WEINSTEINE.Maxi2mumlikelihoodnoisecancellationusingtheEMalgorithm[J].IEEETransactionsonAcoustics,SpeechandSignalProcessing,1989,37(2):2042216.
进展,2000,22(6):17220.
ZHANGChun2ting.Thecurrentstatusandtheprospectofbioinformatics[J].WorldSci2techR&D,2000,22(6):17220.(inChinese)
[2]YOONBJ,VAIDYANATHAN2tificationandofburied].ProcessingMagazine,2007,1)64274.
[3]YOONBJ,VAIDYANATHANPP.Fastsearchofse2
(上接第40页)
参考文献(References):
[1]JINGLi2ping,NGMK,HUANGJZ.Anentropy
weightingk2meansalgorithmforsubspaceclusteringofhigh2dimensionalsparsedata[J].IEEETransactionsonKnowledgeandDataEngineering,2007,19(8):102621041.
[2]HUANGJZ,NGMK,RONGH,etal.Automatedvari2ableweightingink2meanstypeclustering[J].IEEETransPatternAnalysisandMachineIntelligence,2005,27(5):1212.
[3]FRIGUIANDH,NASRAOUIO.Unsupervisedlearningofprototypesandattributeweights[J].PatternRecognition,2004,37(3):5672581.
[4]CHANY,CHINGW,NGMK,etal.Anoptimizationalgorithmforclusteringusingweighteddissimilaritymea2sures[J].PatternRecognition,2004,37(5):9432952.[5]DOMENICONIC.Locallyadaptivetechniquesforpattern
classification[D].NewYork:GeorgeMasonUniversity,2002.
[6]DOMENICONIC,PAPADOPOULOSD,GUNOPULOS
D,etal.Subspaceclusteringofhighdimensionaldata[C]ΠΠProceedingsoftheFourthSIAMInternationalConferenceonDataMining.Florida,USA:[s.n.].2004:5172521.[7]杨善林,李永森,胡笑旋,等.K2means算法中的k值优化
问题研究[J].系统工程理论与实践,2006(2):972101.
YANGShan2lin,LIYong2sen,HUXiao2xuan,etal.Op2timizationstudyonkvalueofk2meansalgorithm[J].Sys2temsEngineering-Theory&Practice,2006(2):972101.(inChinese)
[8]胡清河,张 爽,汪定伟.神经网络在项目评估中的应用[J].东北大学学报:自然科学版,2007,28(2):1692171.HUQing2he,ZHANGShuang,WANGDing2wei.Applica2tionofneuralnetworktoprojectevaluation[J].JournalofNortheasternUniversity:NaturalScience,2007,28(2):1692171.(inChinese)
[9]宋擒豹,沈钧毅.Web页面和客户群体的模糊聚类算法[J].小型微型计算机系统,2001,22(2):2292232.SONGQin2bao,SHENJun2yi.Fuzzyclusteringalgorithmsforwebpagesandcustomersegments[J].Mini2microSys2tems,2001,22(2):2292232.(inChinese)