生物计算技术》第4章多重序列比对分析.ppt

,Biocomputing technology Multiple sequence alignment,第4章 多重序列比对分析,目的要求： 1 掌握多重序列比对的基本概念及意义。 2 掌握多重序列比对的星形比对、树形比对及隐马尔可夫模型。 3 了解多重序列比对的动态规划算法、CLUSTAL W 算法。,教学内容： 4.1 多重序列比对的意义 4.2 多重序列比对算法原理,Multiple sequence alignment,4.1 多重序列比对的意义,目的： 发现多个序列的共性 发现与结构和功能相关的保守序列片段 定义： 设：有k个序列s1, s2, . ,sk，每个序列由同一个字母表中的字符组成，k大于2，通过插入“空位”操作，使得各序列达到一样的长度，从而形成这些序列的多重比对。,Biocomputing technology Multiple sequence alignment,8条免疫球蛋白序列片段的多重比对：,半光氨酸,色氨酸,疏水残基,保守区域,SP得分,Biocomputing technology Multiple sequence alignment,Biocomputing technology Multiple sequence alignment,通过序列的多重比对，可以得到一个序列家族的序列 特征。当给定一个新序列时，根据序列特征，可以判断这 个序列是否属于该家族。 对于多序列比对，现有的大多数算法都基于渐进比对 的思想，在序列两两比对的基础上逐步优化多序列比对的 结果。 进行多序列比对后，可以对比对结果进行进一步处理， 例如构建序列的特征模式，将序列聚类、构建分子进化树等。,4.2 多重序列比对算法原理,4.2.1 SP模型 4.2.2 多重比对的动态规划算法 4.2.3 优化算法 4.2.4 星型比对 4.2.5 树形比对 4.2.6 CLUSTALW算法 4.2.7隐马尔可夫模型,Biocomputing technology Multiple sequence alignment,4.2.1 SP 模型 (Sum-of-Pairs) 逐对加和函数,作用：评价多重序列比对的结果,SP计算的两种方法：,Biocomputing technology Multiple sequence alignment,方法1：先计算多重比对结果的每一列字符的得分， 然后求整体多重比对得分,Biocomputing technology Multiple sequence alignment,假设: 得分函数(代价函数) 具有加和性，即多重比对的得分 是各列得分总和。 思路: 如何给比对的每一列打分，然后将各列的和加起来， 成为一个总得分。 每一列的处理方式: 寻找一个具有k 个变量的打分函数，每一个变量或者是 一个来自特定字母表中的字符，或者是一个空位。 k 是参与多重比对的序列的个数。,Biocomputing technology Multiple sequence alignment,显式函数应满足如下条件： 函数形式简单，具有统一的形式，不随序列的个数 而发生形式的变化。 2. 根据得分函数的意义，函数值应独立于各参数的顺序， 即与待比较的序列先后次序无关。 3. 对相同的或相似字符的比对，奖励的得分值高，而对 于不相关的字符比对或空白，则进行惩罚（得分为负值）。 满足上述条件的一个函数就是常用的逐对加和函数，SP函数 。,方法1：先计算多重比对结果的每一列字符的得分， 然后求整体多重比对得分,其中，c1,c2,ck是一列中的k个字符， p是关于一对字符相似性的打分函数。 SP_score(c1,c2,ck)是多重序列比对中某一列 的得分.,Biocomputing technology Multiple sequence alignment,例：图4.1多重比对的倒数第3列的SP得分：,打分函数： P（a,a）=0 P（a,b）= -1 （ab） P（a,-）=P（-,b）= -1 P（-,-）=0,逐对计算p(1,2)，p (1,3)，.，p(1,8)，p (2,3)，p(2,4)，.p(2,8) .，p(7,8) 的所有得分： （-7-6-5-4-3-2-1）+2 = -26 然后将一个多重比对所有列的得分全部加起来，其和即为该多重比对的得分。 将所有多重比对的得分计算出来进行比较，得分最高的，应该是最好的。,Biocomputing technology Multiple sequence alignment,多重比对在两条特定序列上的投影,Biocomputing technology Multiple sequence alignment,方法2：先计算多重序列结果的序列两两比对得分， 然后计算整体多重比对得分。, 是一个多重比对 ij是由推演出来的序列si和sj的两两比对,方法1和方法2等价的条件：P(-,-)=0,Biocomputing technology Multiple sequence alignment,4.2.2 多重比对的动态规划算法,多重序列比对的最终目标是通过处理得到一个得分最高（或代价最小）的序列对比排列，从而分析各序列之间的相似性和差异。,Biocomputing technology Multiple sequence alignment,4.2.2 多重比对的动态规划算法,s1:VSN-S s2:-SNA- s3:-AS,Biocomputing technology Multiple sequence alignment,前趋节点的个数等于2k - 1,问题： 计算量巨大 时间复杂度为O(2ki=1,.,k si) O(2kNk),Biocomputing technology Multiple sequence alignment,NP-完全问题的定义,Biocomputing technology Multiple sequence alignment,P 类问题为多项式界的问题； NP 类问题是这样一类问题，如果有一个复杂度为多 项式的算法解决其中的某个问题，则所有这些问题都在P类中； 而NP-完全问题是这样一类问题，如果其中的某个问题存在 多项式界的算法，则NP 类中的每一个问题都存在一个多项 式界算法。 NP 完全问题通常被认为是一些人们难以在有限的时间、 空间内对问题求出最佳解的问题，几乎所有专家都认为不可 能在多项式时间内准确求解NP-完全问题。,NP-完全问题的近似求解方法,Biocomputing technology Multiple sequence alignment,舍去寻找最优解的要求，寻找对一般问题比较接近 最优解的近似解； 2. 利用非常规的求解技术求解，例如，利用神经网络、 遗传算法等方法进行问题求解。,生物信息学中NP-完全问题的近似求解方法,Biocomputing technology Multiple sequence alignment,1. 只求解规模比较小的问题； 2. 利用动态规划、分支约束等技术减小搜索空间，提高 求解问题的效率； 3. 针对具体问题的特点，根据实际输入情况，设计实用 求解算法，这样的算法虽然在最坏的情况下其时间复 杂度是非多项式的，但是算法执行的平均效率和复杂 度与多项式的算法相当； 4. 采用近似算法或者启发式方法，如局部搜索、模拟退火、 遗传算法等。 对基于SP 模型寻找最优多重序列比对这样一个问题， 可以用近似的方法求解，其算法的时间复杂度可用多项式表示。,4.2.3 优化计算方法,标准动态规划算法存在的问题： 搜索空间大 剪枝技术：将搜索空间限定在一个较小的区域范围内。 若问题是搜索一条得分最高（或代价最小）的路径，则在搜索时如果当前路径的得分低于某个下限（或累积代价已经超过某个上限），则对当前路径进行剪枝，即不再搜索当前路径的后续空间。,Biocomputing technology Multiple sequence alignment,Biocomputing technology Multiple sequence alignment,在序列两两比对中, Fickett 和Ukkonen 设计了一种 称为定界约束过程 (bounding procedure)的方法来缩小搜 索空间,减少计算量, 其中距离矩阵的上界和下界可以预先 确定或动态变化。 为了在多维空间上使用动态规划算法，Carrillo 和 Lipman 将这种思想引入到多重序列比对，即先进行初步 的序列双重比对，以限制进一步做多重序列全面比对所需 要的多维空间的大小和计算量，从而克服了多序列的维数、 空间和运算量之间的矛盾,Carrillo-Lipman的优化计算方法,Biocomputing technology Multiple sequence alignment,设k 条序列的长度分别为n1n2nk，按照SP 得分模 型计算这些序列的最优比对。依然采用动态规划方法，但 并不计算超晶格空间中所有的节点，而是仅处理与最优路 径“相关”的节点。但是，哪些节点是相关的呢？这需要观 察节点在两条序列上的投影。,确定相关节点的方法： 假设是关于k 条序列s1s2sk 的最优多重比对。 从某个节点向任何两条序列所在的平面投影，如果该投影 是这两条序列两两最优比对的一部分（前面一部分），则 该节点是与最优比对相关的节点。,问题的提出:,一种计算两条序列经过特定断点的最优比对的算法,Biocomputing technology Multiple sequence alignment,设有两条序列s、t, |s|=m, |t|=n; 位点i 将序列s一分为二, 位点j将序列t一分为二, 则:序列 s、t 对于经过特定断点（i、j）的最优比对可分为两个部分: 对应于两条序列前缀0:s:i 与0:t:j 的最优比对; 对应于两条序列后缀 i:s:m与j:t:n 的最优比对;,例:,Biocomputing technology Multiple sequence alignment,比对两条序列: s=ATTCGG, t=GATTC 打分函数: p(a,a)=1, p(a,b)=-1, p(a,-)=p(-,b)=-1,Biocomputing technology Multiple sequence alignment,如果超晶格空间中的一个节点想任意两条序列所在 的平面投影,投影在这些” 断点”中,则超晶格空间中的这 个节点就是与最优路径相关的节点,否则不是相关节点.,小结: 在进行多重序列比对时, 首先要进行序列的两两比对, 其目的就是要找到任意两条序列通过特定断点的最优比对, 找到这些断点,然后,将多重比对中的超晶格空间的节点向 任意两条序列所在的平面投影,看看投影是否在这些断点上, 如果节点向各个平面的投影均在相应的断点上,则这个节点 是与多重序列比对的最优路径相关的节点,否则,就不是相 关节点,要进行”剪枝”处理.,Biocomputing technology Multiple sequence alignment,(1) 设 是关于s1s2sk 的最优比对，如果 SP_score( ) L，则 score( i,j ) Li,j (4-6 ) 其中， score( i,j ) 是 在si 和sj 所在平面投影的得分, 这里,L实际上是最优多重序列比对的一个下限, Lij是序列si 和序列sj 比对得分的一个下限,虽然最优多重比对的投影不一定是两两最优比对， 但是我们可以为投影的得分设立一个下限。,判断超晶格空间中一个节点是否是相关节点的方法：,Biocomputing technology Multiple sequence alignment,(2) 现在的问题: 需要判断超晶格中的一个节点 i = (i1,i2,ik ) 是 否在L 的限制下与最优比对相关。,(3) 简单地说, 如果一个节点的所有投影满足(4-6 )式的 条件，则该节点是相关的。若条件不满足， 即score( i,j )小，则该节点不可能是相关的， 因此，i 肯定不会处于最优路径上。,(4) 公式(4-6)的条件只是一个必要条件，但不是充分条件。 满足条件只是说明i 可能处于最优路径，但不一定处于 最优路径。公式(4-6)条件的作用是限制搜索空间，提高 算法的实施效率。,Biocomputing technology Multiple sequence alignment,(5) 将判断条件公式(4-6)进一步具体化， 则对于所有的1xyk，如果i 满足 Cx,y ix,iy Lx,y （4-7 ） 则i 是相关的。 这里， Cx,y是序列sx 和sy 的总得分矩阵; Cx,y ix,iy 表示在点 ix,iy 处的值,即 经过ix,iy断点的sx和sy的最优比对得分 ; ix,iy是断点; Lx,y是sx 和sy 的比对得分的下限,注意: 尽管不是所有的相关节点均参与最优比对, 但只有与最优路径相关的节点才参与最优比对.,Biocomputing technology Multiple sequence alignment,(6) 判断非相关节点的方法: 假设是关于s1s2sk 的最优比对，其所对应的 路径通过节点 i= ( i1,i2,ix,iy,ik ),则: Cx,y ix,iy Score(i,j) Lx,y 反之, 如果 Cx,y ix,iy Lx,y ,则多重序列最优比对路径不会经过节点 i= ( i1,i2,ix,iy,ik ), 因而,该超晶格节点是非相关节点.,Biocomputing technology Multiple sequence alignment,为了得到一个合理的下限 L，我们可以任选一个包含 所有序列的多重比对，计算其得分，以此作为 L。若选取 的 L 接近于最优值，算法速度将大大提高。 注意：L 的值不能超过最优值，否则算法出错。 在实现上述优化计算方法时，必须非常仔细。不可能 对超晶格中的每一个节点都进行上述判断，否则，计算时 间不会有多大的减少。我们需要一种完全消除无关单元的 办法, 以便不需再处理它们。 下面介绍一种可能的策略:,Biocomputing technology Multiple sequence alignment,在计算机中实现“剪枝”技术的措施-1: 从超晶格的零点0 = (0, 0, , 0) 出发, 此节点总 是相关的, 并影响依赖于它的节点. 每个这样的节点 又有它自己的受影响的节点, 如此展开, 直至达到在 最终角落的节点( n1n2nk ). (2) 以数组a 保存各节点的计算结果. 如果在计算a j 时用到i, 称节点i 影响另一个节点j, 又称,节点j 依赖于节点i。 每个节点依赖于至多2k-1个其它节点,也至多影响 2k-1个其它节点. (3) 节点i 和节点j 关系的另一特征是：b=j- i b是一个非零的二进向量.,Biocomputing technology Multiple sequence alignment,在计算机中实现“剪枝”技术的措施-2:,(4) 为了便于处理，设置一个缓冲区，该缓冲区内仅存放 相关节点的后续节点。 (5) 首先将0 放入其中。 (6) 当节点i 进入缓冲区时，其对应的值a i 被初始化， 然后a i 的值在随后的阶段中被更新。 当节点i 离开缓冲区时，其值即为该点真正的值，并用 于其他节点(依赖于此节点)的计算。 (7) 节点i 的后续节点是否要计算，还取决于i 是否为相关 节点，若不是，则不再计算其后续的其他节点。,具体实现过程:,Biocomputing technology Multiple sequence alignment,设节点j 是一个依赖于节点i 的相关节点， 如果j 不在缓冲区内，则将其放入缓冲区，并计算 a j a i +SP_Score( Colum(s, i, b) ) (3) 如果j 早已在缓冲区中，则按下式更新 a j max( a j , a i +SP_Score( Colum(s,i, b) ) ) 注意: Carrilo-Lipman 算法要求待比较的多个序列具有较大 的相似性，并且序列数不能太多。,4.2.4 星形比对,Biocomputing technology Multiple sequence alignment,* 启发式方法作为首选。,* 启发式方法不一定保证最终能得到最优解，但在大多数 情况下，其计算结果接近于最优结果。,* 启发式这类方法能够大大减少所需的计算时间，加快计 算速度。,* 渐进法是多重序列比对中常用到的启发式方法。其基本 思想是将序列多重比对转化为序列两两比对，逐渐将两 两比对组合起来，最终形成完整的多序列比对。,* 组合两两序列比对的方法有： 星形结构或者树形结构。,1. 星形比对的基本思想：,星形比对是一种启发式方法，首先由Gusfield 提出。 在给定的若干序列中，选择一个核心序列，通过该序列与其它序列的两两比对形成所有序列的多重比对 ，从而使得 在核心序列和任何一个其它序列方向的投影是最优的两两比对。,Biocomputing technology Multiple sequence alignment,2. 星形比对算法概述-1,Biocomputing technology Multiple sequence alignment,* 设s1,s2,sk是k 条待比对的序列。 * 假设已知核心序列是sc，1c k，则可以利用标准的 动态规划算法求出所有sc 和si 的最优两两比对。 这个工作的时间复杂度为O(kn2), 假设所有序列的长度约为n。 * 将这些序列的两两比对聚集起来，并采用“只要是空位，则 永远是空位”的原则。 * 聚集过程从某一个两两比对开始，然后逐步加上其他的两 两比对。在这个过程中，逐步增加sc中的空位字符，以适应 其他的比对，但决不删除sc中已存在的空位字符。,2. 星形比对算法概述-2,Biocomputing technology Multiple sequence alignment,* 随着后续比对的不断加入, 一方面我们有一个由sc指导的、 已经建立好的部分序列的多重比对，另一方面我们有sc与 一个新序列之间的比对。在两种比对中我们插入尽可能少 而必要的空位, 以保持与扩展的sc序列相匹配。然后,将新 扩展的序列加入序列群中, 且它和其它扩展的序列具有相 同的长度。 * 这个过程一直进行到所有的两两比对都加入以后结束，这 一步所需的计算量为O(k2n)。 * 算法总的时间复杂度为O(kn2+k2n)。,sc s1 s2 sk,(sc, s1) (sc, s2) (sc, sk),两两比对 多重比对,Biocomputing technology Multiple sequence alignment,方法1：尝试将每一个序列分别作为核心序列， 进行星形多重序列比对，取比对结果最 好的一个。即SP得分最高的为最好。 方法2：是计算所有的两两比对，取下式值最大 的一个：,3. 如何选择核心序列？,Biocomputing technology Multiple sequence alignment,4.举例,有5个序列： s1 = ATTGCCATT s2 = ATGGCCATT s3 = ATCCAATTTT s4 = ATCTTCTT s5 =ACTGACC,sc=s1,ATTGCCATT ATTGCCATT- ATTGCCATT ATTGCCATT ATGGCCATT ATC-CAATTTT ATCTTC-TT ACTGACC- (s1,s2) (s1,s3) (s1,s4) (s1,s5),ATTGCCATT- ATGGCCATT- ATC-CAATTTT ATCTTC-TT- ACTGACC-,Biocomputing technology Multiple sequence alignment,星形比对是一种多重序列比对近似的方法，然而是一种比较好的近似方法。如果用代价来评判多重序列的比对结果，则可以证明，用该方法所得到多重序列比对的代价不会大于最优多重序列比对代价的两倍。,Biocomputing technology Multiple sequence alignment,4.2.5 树形比对,k个待比对的序列 具有k个叶节点的树 每个叶节点对应一条序列 将序列赋予树的内部节点，可以计算树中每个分支的权值。 权值代表对应分支连接的两个序列之间的相似性。 所有权值的和就是这棵树的得分 树形比对的问题：寻找一种树的内部节点序列赋予方式， 使得树的得分最大。 星型比对可以看作是树形比对的一种特例。,Biocomputing technology Multiple sequence alignment,将CT、CG、CT分别赋予节点 x、y、z，则树的得分为8。,CT,CG,CT,多重序列比对 两两序列比对 合并两个比对（比对的比对） Aligment of AlignmentAA算法,打分函数： P（a,a）=1 P（a,b）= 0 （ab） P（a,-）=P（-,b）= -1,1,1,1,1,2,2,G - T C A T C - G C T G,比对结果：,Biocomputing technology Multiple sequence alignment,AA算法概述-1,Biocomputing technology Multiple sequence alignment,假设有两个多重比对1和2 1代表序列s1,s2,si的多重比对 2代表序列t1,t2,tj的多重比对 并且,( s1,s2,si ) ( t1,t2,tj )= 代表s1和t1的两两比对, 则 计算与相一致的1 和2比对的算法如下:,AA算法概述-2,Biocomputing technology Multiple sequence alignment,标定1的各列, 如果s1在比对中对应位置的编辑操作不是插入或删除,则这些列分别标记为s1对应位置上的字符: a1,a2,al1 s1=l1 标定2的各列, 如果t1在比对中对应位置的编辑操作不是 插入或删除,则这些列分别标记为t1对应位置上的字符: b1,b2,bl2 t1=l2 对a1,a2,al1 和b1,b2,bl2 进行比对, 在所得到的比对中,对于1、2和中原来有插入或删除操 作的位置, 恢复其原有的实际字符或空位字符”-”.,Biocomputing technology Multiple sequence alignment,a1 a2a3a4,b1 b2b3b4b5,对于n个序列的树形比对的基本算法过程如下： （1）初始化，对于每个序列，生成一个叶节点 （2）利用AA算法合并两个节点，形成一个新节点， 合并的结果放在新节点中，原来的两个节点作为 新节点的子节点 （3）反复执行（2），直到形成n个叶节点的树根为止， 根节点中的序列即为最终的多重比对结果。,s1 s2 s3 s4,1,2,Biocomputing technology Multiple sequence alignment,4.2.6 CLUSTAL 算法,Biocomputing technology Multiple sequence alignment,CLUSTAL 算法是一种渐进的比对方法，它是由 D.G.Higgins和 P.M.Sharp 于1988年首次提出的。,目的： 对来自不同物种的功能相同或相似的蛋白进行比对 和聚类分析，可对这些物种的亲缘关系进行判断，并且 对这些蛋白在生物进化过程中的保守性作出估计。,Clustal算法包括三个阶段：,1. 先将多重序列进行两两比对。基于这些比对，计算得到 一个距离矩阵，该矩阵反映每对序列之间的关系。 2. 根据距离矩阵计算产生系统发育树，以此来确定被比较 序列间相似的程度 3. 有了这棵系统发育树的指导，从关系最紧密的两条序列 开始，逐步引入邻近的序列或序列组，并不断重新构建 比对，直到所有序列都被加入为止。,Biocomputing technology Multiple sequence alignment,举例：,Biocomputing technology Multiple sequence alignment,已知三级结构的七个球蛋白序列，分别为：,Hbb_Human : 人的-球蛋白 Hbb_Horse : 马的-球蛋白 Hba_Human : 人的-球蛋白 Hba_Horse : 马的-球蛋白 Myg_phyca : 抹香鲸的血红蛋白 Glb5_petma : 七鳃鳗蓝血红蛋白 Lgb2_Luplu : 羽扇豆的豆血红蛋白,第一阶段：两两比对产生距离矩阵,Biocomputing technology Multiple sequence alignment,用完全动态规划法计算的7个球蛋白序列之间的 77的距离矩阵,第二阶段：产生指导树,Biocomputing technology Multiple sequence alignment,根据第一阶段结果中的距离矩阵，用邻近归并法来 计算产生一棵有分枝长度和序列权重的有根树。,第三阶段：渐进的比对,Biocomputing technology Multiple sequence alignment,这一阶段基本的步骤是通过一系列两两比对来构建更 大的比对序列组。按照指导树中的分支顺序，从有根树的 末梢到树根逐渐进行。根据图4.22的有根树，通过下面的 顺序对序列进行比对：(1)对(2)(3)对(4)(8)对(9) (5)对(10)(6)对(11)(7)对(12)。 * 每一阶段使用了有残基权重矩阵和空位开放及空位扩展罚 分的完全动态规划算法。 * 每一阶段都由对两个已经存在的排列或序列进行比对组成。 * 在先前比对中出现的空位仍然是固定的,图4.23,Biocomputing technology Multiple sequence alignment,序列权重的作用: 计算多重序列比对得分,Biocomputing technology Multiple sequence alignment,peeksavtal geekaavlal padktnvkaa aadktnvkaa egewqlvlhv aaektkirsa,多重序列比对的统计特征分析,对于所得到的多重序列比对，我们往往需要进行归纳分析，总结这些序列的特征，或者给出这些序列共性的表示。,表示方式1：保守序列表示方式 表示出序列每个位置上最可能出现的字符 或所有可能出现的字符 表示方式2；特征统计图谱（profile）表示方式 （或特征统计矩阵表示方法）,Biocomputing technology Multiple sequence alignment,令P=(P1,P2,Pj,PL )；P表示在多重比对的每一列上各种 字符出现的概率分布 Pj=(pj0,pj1,pj|A|)；A代表字母表， Pjk代表字母表A中第k个字符在第j列 出现的概率。 pj0：第0个字符是特殊的空位符号“-”。 称P为多重序列比对的特征统计图谱或特征统计矩阵。 P反映了多重序列比对各个列的特征。,表示方法2：特征统计图谱,Biocomputing technology Multiple sequence alignment,1 2 3 4 5 (位置) A 0.8 0.2 0.2 0.6 0.0 T 0.0 0.4 0.6 0.4 1.0 C 0.2 0.2 0.2 0.0 0.0 G 0.0 0.2 0.0 0.0 0.0 （碱基）,Biocomputing technology Multiple sequence alignment,A T T A T A A C T T C T T A T A C T T T A G A A T,利用保守序列或者特征统计图可以判断一个序列是否满足一定的特征 给定一个序列s=a1a2am，定义字符a在第j位的代价为 其中，|A|代表字母表A的长度， Ak代表A的第k个字符， A0代表空缺字符“-”。整个序列s的代价为,一条序列与特征统计图相对照，如果代价值小，说明该序列具有相应的特征，否则该序列不具备相应的特征。,Biocomputing technology Multiple sequence alignment,4.2.7 隐马尔可夫模型 HMM,Biocomputing technology Multiple sequence alignment,HMMHidden Markov Model,1. 马尔可夫链（复习）, 具有马尔可夫性的离散状态随机过程称为马尔可夫链。 马尔可夫性,即无后效性,若已知现在的状态,将来与过去无关。 离散状态的马尔可夫链的定义: 考虑只取有限个或可数个状态的随机过程: Xn,n=0,1,2, 假设对一切状态:i1,i2,in-1,i,j,一切n0,有: P( Xn+1=jXn=i,X0=i0 ) = P( Xn+1=jXn=i ) 成立, 则称随机过程: Xn,n=0,1,2, 为离散状态的马尔可夫链。,1. 马尔可夫链（复习）,Biocomputing technology Multiple sequence alignment, P( Xn+1=jXn=i ) 称为马尔可夫链的一步转移概率, 记为Pij(n,n+1)。 若转移概率Pij(n,n+1)不依赖n,则可简记为Pij,则称此 马尔可夫链是时齐的,否则是非时齐的. 用P记转移概率所对应的矩阵,即转移矩阵(transition matrix) 一个马尔可夫链的概率分布完全由它的转移矩阵与初始 分布决定.,2. 隐马尔可夫模型 HMM,Biocomputing technology Multiple sequence alignment,隐马尔可夫模型是一种概率模型,它是由马尔可夫链发 展扩充而来的一种随机模型。这种方法已经成功地应用于多 个领域,如语音识别、光学字符识别等。 HMM在生物信息学领域中也有着重要应用。如：序列 分析，基因识别等。 HMM可以被理解为一个双重随机过程： (1) 系统状态变化的随机过程， (2) 由状态决定输出的随机过程。 一个隐马尔可夫链 Xt,Yt 包含两部分： 一个潜在的、不可观察的有限状态马尔可夫链 Xt 一个外显的、可观察的随机过程 Yt ， Yt 的分布依赖于Xt。,2. 隐马尔可夫模型 HMM,Biocomputing technology Multiple sequence alignment,HMM是用概率统计的方法来描述时变信号的过程. HMM是一个动态的统计模型，可以用有限状态机FSM来描述。 ( FSMfinite state machine ) 有限状态机可以从一种状态转移到另一种状态，每个状态 转换有不同的概率。某状态是否转移到下一状态取决于该状态 的状态转移概率，而在某一状态下能看到哪一个观测值，取决 于该状态的观测概率。,隐马尔可夫模型 HMM,Biocomputing technology Multiple sequence alignment,记HMM模型为: 其中: A为状态转移概率矩阵， B为观察概率矩阵， 为初始状态分布。 HMM 模型的建模工作主要为确定这三个参数。,HMM的三个经典问题:,Biocomputing technology Multiple sequence alignment,问题1（评测问题，evaluation）： 已知模型 和输出序列O，求由 生成O的概率。 问题2（译解问题，decoding）： 已知模型 和输出序列O，求最有可能生成O的状态 转移序列。 问题3（学习问题，learning）： 已知模型 和输出序列O，求最有可能生成O时模型 的参数。,Profile概形、谱,Biocomputing technology Multiple sequence alignment,* 概形是对一组序列进行全局多重比对时被发现的,是将比对 中具有较高保守区域变成小的多重比对,然后得到多重比对 记分矩阵. * 概形由更像小的多重排列的列构成,可以包括: 匹配、失配、 插入、缺失. * 概形一旦生成,就可用于寻找一个可能与之匹配的目标序列, 它利用表中记分来评价每个位置的可能性. 例: 25个氨基酸长的概形表格,有25列,每列将有20个记分值. 一列中每个匹配氨基酸的记分都在概形中对应的位置上. 缺点：偏向性,Profile HMM (1) 模型结构,Biocomputing technology Multiple sequence alignment,* 对于生物序列而言,HMM的字符当然是20个字母的氨基酸 或4个字母的核苷酸。但依据不同的问题,其它的一些字符 也可以使用,如64个密码子的三联体字母,3个字母(, coil )的二级结构等. * 编码蛋白质的原始DNA序列,在生物的进化过程中会受到 自然环境和各种因素的影响,使翻译出的蛋白质序列经历 突变、遗失或引入外源序列等变化,最后按不同的进化路 径分化,形成多种功能相近的蛋白质。 所以,可以把这些蛋白质看作由一个基本蛋白质序列经过插 入、删除或替换了某些氨基酸残基而形成的。这个过程可以 用HMM来表示。,图4.9,Biocomputing technology Multiple sequence alignment,图中: 方形代表匹配状态(主状态)，即输出的氨基酸和基本序列中对应 的氨基酸匹配; 圆形表示删除或缺失状态，即从原始蛋白质序列中去掉一个特定 的氨基酸。 菱形表示氨基酸的插入，即在原始蛋白质序列插入一个氨基酸。 各状态间的箭头表示状态间的转换途径。,注意: 不同的参数会使模型以不同的概率产生新的氨基酸。 一个训练好的模型可以代表有共同特征的蛋白质序列。,图4.10,Biocomputing technology Multiple sequence alignment,Profile HMM 与标准的Profile的比较,Biocomputing technology Multiple sequence alignment, Profile HMM有正规的概率作基础，对于序列的删除和 插入状态的记分也有较为可靠的理论依据。而标准的 Profile纯粹是一种启发式的方法。 HMM用统计方法估计序列某一位点核苷酸或氨基酸残基 出现的真正概率，而标准的Profile却是用自身的观察频率 给核苷酸或氨基酸残基指派分数。 由于,Profile HMM方法从10至20个核苷酸序列构成的 比对中提取的信息，相当于用标准的Profile从40至50个 核苷酸序列构成的比对中提取的信息。,(2) 用HMM给序列打分,Biocomputing technology Multiple sequence alignment,* 因训练好的HMM代表一个蛋白质家族，我们可以用它对序 列进行打分，根据分值来衡量这条序列是否属于由此HMM 代表的蛋白质家族。 得分高于阈值，证明这条序列是这个家族的成员； 得分低于阈值，说明它不是家族的成员。,HMM用于分析蛋白质序列的原理:,Biocomputing technology Multiple sequence alignment,分析蛋白质产生不同序列的概率。 对于与模型相符合的序列，能以较大的概率产生该序列； 若不与该模型符合的序例，则按此模型产生该序列的概率会 较小。 任何一条序列都可以由HMM中的一条路径所呈现。 对于给定的模型，计算产生某条序列的概率就是计算这条 路径上的所有输出概率与转移概率的乘积。,图4.11,Biocomputing technology Multiple sequence alignment,Biocomputing technology Multiple sequence alignment,* 如果已知确切的状态路径，这个计算是很简单的。 * 在一个实际模型中，可以通过不同的状态路径产生同一 条序列。因此，产生一条序列的正确的概率是所有可能 的路径产生该序列的概率的总和。 * 这种计算复杂并且难以实施，除非是一条很短的序列。 * 以下介绍两种很好的替代方法： Viterbi算法可以计算出由模型产生某序列的最有 可能的路径。 前向算法(forward algorithm) 用于计算所有路径的概率和。, Viterbi算法,Biocomputing technology Multiple sequence alignment,作为一种优化算法，动态规划算法在生物信息处理中 得到广泛的应用。前期课程，我们已经介绍了动态规划算 法在序列比较中的应用；这里，介绍应用动态规划算法求 解HMM模型的最优路径的方法。这个方法就是著名的 Viterbi算法。,图4.12,D1,D2,D3,I0,I1,I2,I3,M1,M2,M3,0.06,Biocomputing technology Multiple sequence alignment,最优路径：开始 I0 M1 M2 M3 结束 在这条路径上产生“ACCY”的概率为：5.9728*10-5,图4.13,0.0015,0.0046,0.4850,0.2231,Biocomputing technology Multiple sequence alignment, 前向算法 forward algorithm,Biocomputing technology Multiple sequence alignment,给定一个HMM模型M和一个字符序列X(x1,x2,xl),假定 产生序列X的对应路径未知,要求计算模型M产生X的概率 P(X/M). * Viterbi算法是在可以产生序列X的各种路径中,选择一条 最优*,使得P(X/*)最大. * 现在的问题: 既然有多条路径可以产生序列X,那么,模型M 产生序列X总的可能性有多大? * 解决的方法: 类似于Viterbi算法,用求和替代求最大值.,图4.14,1.8*10-4,5.52*10-4,3.0912*10-4,3.3849*10-7,6.8965*10-5,Biocomputing technology Multiple sequence alignment,图4.12的HHM产生“ACCY”序列的概率为： 3.3849*10-7+6.8965*10-5=6.9303*10-5, 局部和全局打分,Biocomputing technology Multiple sequence alignment,有时用局部打分会更合理，即序列的分值由子序列的 最高分值所替代。,例：图4.11中ACCY序列，将概率转换为对数形式，序列 的全局打分就是以下四个分值之和：-13.25。 Score()=ln（0.4）+ ln（0.3）= -2.12 Score()=ln（0.46）+ ln（0.6）= -1.29 Score()=ln（0.97）+ ln（0.5）= -0.72 Score()=ln（0.015）+ ln（0.73）+ ln（0.01）= -9.12 (-2.12)+(-1.29)+(-0.72)+(-9.12)= -13.25,分值太低了，以至于不能判定ACCY是否为被模拟的蛋白质 家族的成员。,Biocomputing technology Multiple sequence alignment,假设ACCY是图4.15中所示蛋白质家族的成员，在这 种情况下，全局分就不是一个很好的衡量标准。如果采用 局部打分，最后分值就会高得多。我们发现最高分的子序 列为CC，分值为 2.01。这个分值就足以判定ACCY为此 家族的成员。这时，局部打分比全局打分更精确。,(3) 模型的训练,Biocomputing technology Multiple sequence alignment, 根据训练用蛋白质序列的平均长度确定模型中各个状态 的数量而建立HMM。 选用训练算法对HMM进行训练，即调整该模型的参数 （转移概率和输出概率），使得该模型适用于训练所用 的序列，并能够以最大的可能性产生参与训练的观察序 列。,(3) 模型的训练-最大似然法 ML,Biocomputing technology Multiple sequence alignment,给定在训练数据集中的一组序列： s(1)，s(2)，s()， 一个模型究竟在多大程度上适合该数据集，可由根据该 模型观测到每一个序列的联合概率来表征：,其中， 是模型产生第j个序列的概率。 上式称为该模型的似然，最大似然（maximum likelihood，ML） 其原理可用于测量模型的性能，即选取模型参数使上式值最大。,(3) 模型的训练-最大后验概率MAP,Biocomputing technology Multiple sequence alignment,基本思想：在给定序列下，使模型的后验概率最大。,假设：对所有可能的模型参数存在优先概率分布，这种 概率分布体现模