生物信息学实验市公开课一等奖百校联赛特等奖课件

生物信息学试验试验2隐马尔科夫模型上海交通大学生命科学技术学院生物信息学与生物统计学系10/10/1第1页生物学中惯用统计模型 StructuredprobabilitymodelsMarkovmodelsHiddenmarkovmodelsArtificialNeuralNetwork(A.N.N)

10/10/2第2页IntroductionHiddenMarkovModels(HMMs)最早是在上个世纪60年代末70年代初提出来。进入80年代以后，逐步被利用在各个领域。10/10/3第3页IntroductionHiddenMarkovModels作为一个强有力统计学模型，主要被应用在一些连续行或时间延续性事件建模上语音识别系统。生物学中DNA/protein序列分析机器人控制。文本文件信息提取。10/10/4第4页HMM优点1，它数学结构非常丰富，适合用于各个领域研究。2，在很多领域中，已经证实它结果和实际符合相当好。10/10/5第5页ProbabilityReview10/10/6第6页独立事件概率构想我们做一连串试验，而每次试验所可能发生结果定为E1,E2,…En,…。（可能是有限也可能是无限）。每一个结果Ek，假如给定一个出现可能性pk（即概率），则某一特定样本之序列Ej1

Ej2…Ejn出现概率为p(Ej1

Ej2…Ejn)=pj1…Pjn。10/10/7第7页马尔科夫链普通及惯用统计中，彼此相互「独立」大约是最有用一个观念。用简单术语來说，相互「独立」就是彼此毫不相干，一点牵涉都沒有。不过实际生活中很多事件是相互关联[不是相互独立」也就是相互关联意思，不过要怎样相关呢？怎样在相关中作一些简单分类呢？马尔科夫链就是要描述在「相关」这个概念中最简单一个。但即使如此，相关马可夫链理论已经相当丰富了。在概率理论中，它几乎占了绝大部分。10/10/8第8页马尔科夫链在马尔科夫链中考虑最简单「相关」性。在在这种情况下，我们不能给任一个事件Ej一個概率pj但我们给一对事件(Ej,Ek)一個概率pjk，这个时候pjk解释是一个条件概率，就是假设在某次试验中Ej已经出现，而在下一次试验中Ek出现概率。除了pjk之外，还需要知道第一次试验中Ej出現機率aj。有了这些资料后，一個样本序列Ej0

Ej1…Ejn（也就是说第零次试验结果是Ej0，第一次一次是Ej1……第n次试验是Ejn）概率就很清楚是P(Ej0,Ej1,Ejn)=aj

pj0j1

pj1j2…pjn-1jn。10/10/9第9页隐马尔科夫模型不过在大多数情况下我们所观察到值并不是序列本身元素。即观察值不等于状态值。故我们引入隐马尔科夫模型。10/10/10第10页定义一个HMM是一个五元组：(ΩX,ΩO,A,B,π)其中：ΩX={q1,...qN}：状态有限集合ΩO={v1,...,vM}：观察值有限集合A={aij}，aij=p(Xt+1=qj|Xt=qi)：转移概率B={bik}，bik=p(Ot=vk|Xt=qi)：输出概率π={πi}，πi=p(X1=qi)：初始状态分布10/10/11第11页假设对于一个随机事件，有一个观察值序列：O1,...,OT该事件隐含着一个状态序列：X1,...,XT假设1：马尔可夫假设（状态组成一阶马尔可夫链）

p(Xi|Xi-1…X1)=p(Xi|Xi-1)假设2：不动性假设（状态与详细时间无关）p(Xi+1|Xi)=p(Xj+1|Xj)，对任意i,j成立假设3：输出独立性假设（输出仅与当前状态相关）

p(O1,...,OT|X1,...,XT)=Πp(Ot|Xt)10/10/12第12页马尔科夫链Vs隐马尔科夫模型Markovchainshaveentirelyobservablestates.Howevera“HiddenMarkovModel”isamodelofaMarkovSourcewhichadmitsanelementeachtimeslotdependinguponthestate.Thestatesarenotdirectlyobserved10/10/13第13页Problems令λ={A,B,π}为给定HMM参数，令σ=O1,...,OT为观察值序列，隐马尔可夫模型（HMM）三个基本问题：评定问题：对于给定模型，求某个观察值序列概率p(σ|λ)；forwardalgorithm解码问题：对于给定模型和观察值序列，求可能性最大状态序列；viterbialgorithm学习问题：对于给定一个观察值序列，调整参数λ，使得观察值出现概率p(σ|λ)最大。Forward-backwardalgorithm10/10/14第14页SolutionsEvaluationproblem：forwardalgorithm定义向前变量采取动态规划算法，复杂度O(N2T)Decodingproblem：Viterbialgorithm采取动态规划算法，复杂度O(N2T)Learningproblem：forward-backwardalgorithmEM算法一个特例，带隐变量最大似然预计10/10/15第15页StructHMMtypedefstruct{

/*numberofstates;Q={1,2,...,N}*/

intN; /*numberofobservationsymbols;V={1,2,...,M}*/ intM; /*A[1..N][1..N].a[i][j]isthetransitionprobofgoingfromstatei*attimettostatejattimet+1*/

double**A; /*B[1..N][1..M].b[j][k]istheprobabilityofobservingsymbolkinstatej*/

double**B; /*pi[1..N]pi[i]istheinitialstatedistribution.*/ double*pi;

}HMM;10/10/16第16页算法：向前算法（1）10/10/17第17页算法：向前算法（2）定义前向变量为HMM在时间t输出序列O1…Ot，而且位于状态Si概率：10/10/18第18页算法：向前算法（3）迭代公式为：结果为：10/10/19第19页Forwardalgorithm10/10/20第20页算法：向后算法（1）10/10/21第21页算法：Viterbi算法（1）TheViterbialgorithmisadynamicprogrammingalgorithmthatcomputesthemostlikelystatetransitionpathgivenanobservedsequenceofsymbols.Itisactuallyverysimilartotheforwardalgorithm。10/10/22第22页Viterbialgorithm10/10/23第23页Viterbiinc/*1.Initialization*/

for(i=1;i<=phmm->N;i++){

delta[1][i]=phmm->pi[i]*(phmm->B[i][O[1]]);

psi[1][i]=0;

}

/*2.Recursion*/

for(t=2;t<=T;t++){

for(j=1;j<=phmm->N;j++){

maxval=0.0;

maxvalind=1;

for(i=1;i<=phmm->N;i++){

val=delta[t-1][i]*(phmm->A[i][j]);

if(val>maxval){

maxval=val;

maxvalind=i;

}

}

delta[t][j]=maxval*(phmm->B[j][O[t]]);

psi[t][j]=maxvalind;

}

}10/10/24第24页生物学中数学模型10/10/25第25页马氏链10/10/26第26页马氏链10/10/27第27页马氏链10/10/28第28页隐马可夫模型10/10/29第29页隐马可夫模型10/10/30第30页隐马可夫模型profile10/10/31第31页RelatedsoftwareHMMER/SAM(SequenceAlignmentandModelingSystem)/

HMMproAwindowsversionforHMMTheDivisionofBiomedicalInformaticsatCincinnatiChildren'sHospitalMedicalCentermetaMEME:AmotifbasedHiddenMarkovModel10/10/32第32页HMMERProfilehiddenMarkovmodels(profileHMMs)canbeusedtodosensitivedatabasesearchingusingstatisticaldescriptionsofasequencefamily'sconsensus.HMMERisafreelydistributableimplementationofprofileHMMsoftwareforproteinsequenceanalysis.ThecurrentversionisHMMER2.3.2(3Oct),containingminorbugfixesandupdatesfortheMayreleaseofHMMER2.3.

10/10/33第33页HMMER10/10/34第34页HowtocreateaHMM多序列比对相关序列选取模型构建模型训练参数调整应用确立模型10/10/35第35页Example:1.Sequenceselection选取相关序列10/10/36第36页2.AlignmentSaveresultasmsfformat多序列比对10/10/37第37页模型建立3.Hmmbuild4.Hmmt5.Hmmcalibrate模型建立用相关序列对模型进行训练参数调整10/10/38第38页模型文件（1）HMMER2.0[2.3.2]NAMEglobins50LENG162ALPHAminoRFnoCSnoMAPyesCOM./hmmbuildglobins.hmmglobins50.msfNSEQ50DATEThuSep1800:02:14CKSUM4694XT-8455-4-1000-1000-8455-4-8455-4NULT-4-8455NULE595-155885338-294453-1158197249902-1085-142-21-31345531201384-1998-64410/10/39第39页模型文件（2）模型部分：HMMACDEFGHIKLMNPQRSTVWYm->mm->im->di->mi->id->md->db->mm->e-222*-28071-1412-1712-339-321-1729113-1457261-1493-15911181-1737-32-1359-178877-13532620-2119-16974--149-50023343-381399106-626210-466-7202753944596359117-369-294-249--1909-8804-451-894-1115-701-1378-110*2-1118-1371-1805-1237-1464-2231-88925282067-899-510-1267-2325-644-266-1422-1057-63-1884-14865--149-50023343-381399106-626210-466-7202753944596359117-369-294-249--18-6914-7956-894-1115-3550-129**…10/10/40第40页6.未知序列搜索查询Hmmsearch:searchasequenceagainsttheprofileHMM未知查询序列Artemia.faProfileHMM：Globin.hmmCommand:hmmsearchglobin.hmmArtemia.fa查询程序查询未知序列文件所用模型查询命令10/10/41第41页查询结果结果分为2个部分1：说明部分（数听说明、选项、模型说明…）2:结果序列部分10/10/42第42页Result1第一部分：相关信息说明软件信息：版本、权限等HMM文件名称，查询阈值等HMM文件一些描述信息10/10/43第43页Result2.1HIT序列分值，E值，domain数目HITdomains分值、位置、E值等信息10/10/44第44页Result2.2高分匹配序列比对10/10/45第45页Result2.3全部序列HIT分值、E值图形分布10/10/46第46页Result2.4结果统计数据10/10/47第47页ApplicationofHMM:pfam10/10/48第48页ApplicationofHMMTMHMM:Predictionoftransmembranehelicesinproteinshttp://www.cbs.dtu.dk/services/TMHMM/10/10/49第49页PFAMPfamisalargecollectionofproteinmultiplesequencealignments

andprofilehiddenMarkovmodels.PfamisavailableontheWorld

WideWebintheUKathttp://www.sanger.ac.uk/Software/Pfam/,

inSwedenathttp://www.cgb.ki.se/Pfam/,inFranceathttp://pfam.jouy.inra.fr/

andintheUSat/.10/10/50第50页PfamIntroductionPfamisadatabaseofproteindomainfamilies.Pfamcontains

curatedmultiplesequencealignmentsforeachfamily,aswell

asprofilehiddenMarkovmodels(profileHMMs)forfindingthese

domainsinnewsequences.Pfamcontainsfunctionalannotation,

literaturereferencesanddatabaselinksforeachfamily.10/10/51第51页PfamIntroductionVersion14.0,June,7459families

22336uniquePfam-AdomainarchitecturesTwobigfamiliesPfam-A:Ahigh-qualitymanualpartofPfam.Pfam-B:Low-qualityautomaticallygeneratedalignmentsofsequenceclustersinSWISSPROTandTrEMBLthatarenotmodelledinthecuratedpartofPfam.10/10/52第52页PfamIntroductionThere

aretwomultiplealignmentsforeachPfamfamily,theseedalignment

thatcontainsarelativelysmallnumberofrepresentativemembers

ofthefamilyandthefullalignmentthatcontainsallmembers

inthedatabasethatcanbedetected.Allalignmentsusesequences

takenfrompfamseq,whichisanon-redundantproteinsetcomposed

ofSWISS-PROTandSP-TrEMBL.TheprofileHMMisbuiltfromthe

seedalignmentusingtheHMMERpackage,whichisthenusedtosearchthepfamseqsequencedatabase10/10/53第53页PfamGoalsOneofthemaingoalsofPfamwastoaidthea