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

Biocomputingtechnology Multiplesequencealignment 第4章多重序列比对分析 目的要求 1掌握多重序列比对的基本概念及意义 2掌握多重序列比对的星形比对 树形比对及隐马尔可夫模型 3了解多重序列比对的动态规划算法 CLUSTALW算法 教学内容 4 1多重序列比对的意义4 2多重序列比对算法原理 Multiplesequencealignment 4 1多重序列比对的意义 目的 发现多个序列的共性发现与结构和功能相关的保守序列片段定义 设 有k个序列s1 s2 sk 每个序列由同一个字母表中的字符组成 k大于2 通过插入 空位 操作 使得各序列达到一样的长度 从而形成这些序列的多重比对 Biocomputingtechnology Multiplesequencealignment 8条免疫球蛋白序列片段的多重比对 半光氨酸 色氨酸 疏水残基 保守区域 SP得分 Biocomputingtechnology Multiplesequencealignment Biocomputingtechnology Multiplesequencealignment 通过序列的多重比对 可以得到一个序列家族的序列特征 当给定一个新序列时 根据序列特征 可以判断这个序列是否属于该家族 对于多序列比对 现有的大多数算法都基于渐进比对的思想 在序列两两比对的基础上逐步优化多序列比对的结果 进行多序列比对后 可以对比对结果进行进一步处理 例如构建序列的特征模式 将序列聚类 构建分子进化树等 4 2多重序列比对算法原理 4 2 1SP模型4 2 2多重比对的动态规划算法4 2 3优化算法4 2 4星型比对4 2 5树形比对4 2 6CLUSTALW算法4 2 7隐马尔可夫模型 Biocomputingtechnology Multiplesequencealignment 4 2 1SP模型 Sum of Pairs 逐对加和函数 作用 评价多重序列比对的结果 SP计算的两种方法 Biocomputingtechnology Multiplesequencealignment 方法1 先计算多重比对结果的每一列字符的得分 然后求整体多重比对得分 Biocomputingtechnology Multiplesequencealignment 假设 得分函数 代价函数 具有加和性 即多重比对的得分是各列得分总和 思路 如何给比对的每一列打分 然后将各列的和加起来 成为一个总得分 每一列的处理方式 寻找一个具有k个变量的打分函数 每一个变量或者是一个来自特定字母表中的字符 或者是一个空位 k是参与多重比对的序列的个数 Biocomputingtechnology Multiplesequencealignment 显式函数应满足如下条件 函数形式简单 具有统一的形式 不随序列的个数而发生形式的变化 2 根据得分函数的意义 函数值应独立于各参数的顺序 即与待比较的序列先后次序无关 3 对相同的或相似字符的比对 奖励的得分值高 而对于不相关的字符比对或空白 则进行惩罚 得分为负值 满足上述条件的一个函数就是常用的逐对加和函数 SP函数 方法1 先计算多重比对结果的每一列字符的得分 然后求整体多重比对得分 其中 c1 c2 ck是一列中的k个字符 p是关于一对字符相似性的打分函数 SP score c1 c2 ck 是多重序列比对中某一列的得分 Biocomputingtechnology Multiplesequencealignment 例 图4 1多重比对的倒数第3列的SP得分 打分函数 P a a 0P a b 1 a b P a P b 1P 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然后将一个多重比对所有列的得分全部加起来 其和即为该多重比对的得分 将所有多重比对的得分计算出来进行比较 得分最高的 应该是最好的 Biocomputingtechnology Multiplesequencealignment 多重比对在两条特定序列上的投影 Biocomputingtechnology Multiplesequencealignment 方法2 先计算多重序列结果的序列两两比对得分 然后计算整体多重比对得分 是一个多重比对 ij是由 推演出来的序列si和sj的两两比对 方法1和方法2等价的条件 P 0 Biocomputingtechnology Multiplesequencealignment 4 2 2多重比对的动态规划算法 多重序列比对的最终目标是通过处理得到一个得分最高 或代价最小 的序列对比排列 从而分析各序列之间的相似性和差异 Biocomputingtechnology Multiplesequencealignment 4 2 2多重比对的动态规划算法 s1 VSN Ss2 SNA s3 AS Biocomputingtechnology Multiplesequencealignment 前趋节点的个数等于2k 1 问题 计算量巨大时间复杂度为O 2k i 1 k si O 2kNk Biocomputingtechnology Multiplesequencealignment NP 完全问题的定义 Biocomputingtechnology Multiplesequencealignment P类问题为多项式界的问题 NP类问题是这样一类问题 如果有一个复杂度为多项式的算法解决其中的某个问题 则所有这些问题都在P类中 而NP 完全问题是这样一类问题 如果其中的某个问题存在多项式界的算法 则NP类中的每一个问题都存在一个多项式界算法 NP完全问题通常被认为是一些人们难以在有限的时间 空间内对问题求出最佳解的问题 几乎所有专家都认为不可能在多项式时间内准确求解NP 完全问题 NP 完全问题的近似求解方法 Biocomputingtechnology Multiplesequencealignment 舍去寻找最优解的要求 寻找对一般问题比较接近最优解的近似解 2 利用非常规的求解技术求解 例如 利用神经网络 遗传算法等方法进行问题求解 生物信息学中NP 完全问题的近似求解方法 Biocomputingtechnology Multiplesequencealignment 1 只求解规模比较小的问题 2 利用动态规划 分支约束等技术减小搜索空间 提高求解问题的效率 3 针对具体问题的特点 根据实际输入情况 设计实用求解算法 这样的算法虽然在最坏的情况下其时间复杂度是非多项式的 但是算法执行的平均效率和复杂度与多项式的算法相当 4 采用近似算法或者启发式方法 如局部搜索 模拟退火 遗传算法等 对基于SP模型寻找最优多重序列比对这样一个问题 可以用近似的方法求解 其算法的时间复杂度可用多项式表示 4 2 3优化计算方法 标准动态规划算法存在的问题 搜索空间大剪枝技术 将搜索空间限定在一个较小的区域范围内 若问题是搜索一条得分最高 或代价最小 的路径 则在搜索时如果当前路径的得分低于某个下限 或累积代价已经超过某个上限 则对当前路径进行剪枝 即不再搜索当前路径的后续空间 Biocomputingtechnology Multiplesequencealignment Biocomputingtechnology Multiplesequencealignment 在序列两两比对中 Fickett和Ukkonen设计了一种称为定界约束过程 boundingprocedure 的方法来缩小搜索空间 减少计算量 其中距离矩阵的上界和下界可以预先确定或动态变化 为了在多维空间上使用动态规划算法 Carrillo和Lipman将这种思想引入到多重序列比对 即先进行初步的序列双重比对 以限制进一步做多重序列全面比对所需要的多维空间的大小和计算量 从而克服了多序列的维数 空间和运算量之间的矛盾 Carrillo Lipman的优化计算方法 Biocomputingtechnology Multiplesequencealignment 设k条序列的长度分别为n1n2 nk 按照SP得分模型计算这些序列的最优比对 依然采用动态规划方法 但并不计算超晶格空间中所有的节点 而是仅处理与最优路径 相关 的节点 但是 哪些节点是相关的呢 这需要观察节点在两条序列上的投影 确定相关节点的方法 假设 是关于k条序列s1s2 sk的最优多重比对 从某个节点向任何两条序列所在的平面投影 如果该投影是这两条序列两两最优比对的一部分 前面一部分 则该节点是与最优比对相关的节点 问题的提出 一种计算两条序列经过特定断点的最优比对的算法 Biocomputingtechnology Multiplesequencealignment 设有两条序列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的最优比对 例 Biocomputingtechnology Multiplesequencealignment 比对两条序列 s ATTCGG t GATTC打分函数 p a a 1 p a b 1 p a p b 1 Biocomputingtechnology Multiplesequencealignment 如果超晶格空间中的一个节点想任意两条序列所在的平面投影 投影在这些 断点 中 则超晶格空间中的这个节点就是与最优路径相关的节点 否则不是相关节点 小结 在进行多重序列比对时 首先要进行序列的两两比对 其目的就是要找到任意两条序列通过特定断点的最优比对 找到这些断点 然后 将多重比对中的超晶格空间的节点向任意两条序列所在的平面投影 看看投影是否在这些断点上 如果节点向各个平面的投影均在相应的断点上 则这个节点是与多重序列比对的最优路径相关的节点 否则 就不是相关节点 要进行 剪枝 处理 Biocomputingtechnology Multiplesequencealignment 1 设 是关于s1s2 sk的最优比对 如果SP score L 则score i j Li j 4 6 其中 score i j 是 在si和sj所在平面投影的得分 这里 L实际上是最优多重序列比对的一个下限 Lij是序列si和序列sj比对得分的一个下限 虽然最优多重比对的投影不一定是两两最优比对 但是我们可以为投影的得分设立一个下限 判断超晶格空间中一个节点是否是相关节点的方法 Biocomputingtechnology Multiplesequencealignment 2 现在的问题 需要判断超晶格中的一个节点i i1 i2 ik 是否在L的限制下与最优比对相关 3 简单地说 如果一个节点的所有投影满足 4 6 式的条件 则该节点是相关的 若条件不满足 即score i j 小 则该节点不可能是相关的 因此 i肯定不会处于最优路径上 4 公式 4 6 的条件只是一个必要条件 但不是充分条件 满足条件只是说明i可能处于最优路径 但不一定处于最优路径 公式 4 6 条件的作用是限制搜索空间 提高算法的实施效率 Biocomputingtechnology Multiplesequencealignment 5 将判断条件公式 4 6 进一步具体化 则对于所有的1 x y k 如果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的比对得分的下限 注意 尽管不是所有的相关节点均参与最优比对 但只有与最优路径相关的节点才参与最优比对 Biocomputingtechnology Multiplesequencealignment 6 判断非相关节点的方法 假设 是关于s1s2 sk的最优比对 其所对应的路径通过节点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 因而 该超晶格节点是非相关节点 Biocomputingtechnology Multiplesequencealignment 为了得到一个合理的下限L 我们可以任选一个包含所有序列的多重比对 计算其得分 以此作为L 若选取的L接近于最优值 算法速度将大大提高 注意 L的值不能超过最优值 否则算法出错 在实现上述优化计算方法时 必须非常仔细 不可能对超晶格中的每一个节点都进行上述判断 否则 计算时间不会有多大的减少 我们需要一种完全消除无关单元的办法 以便不需再处理它们 下面介绍一种可能的策略 Biocomputingtechnology Multiplesequencealignment 在计算机中实现 剪枝 技术的措施 1 从超晶格的零点0 0 0 0 出发 此节点总是相关的 并影响依赖于它的节点 每个这样的节点又有它自己的受影响的节点 如此展开 直至达到在最终角落的节点 n1n2 nk 2 以数组a 保存各节点的计算结果 如果在计算a j 时用到i 称节点i影响另一个节点j 又称 节点j依赖于节点i 每个节点依赖于至多2k 1个其它节点 也至多影响2k 1个其它节点 3 节点i和节点j关系的另一特征是 b j ib是一个非零的二进向量 Biocomputingtechnology Multiplesequencealignment 在计算机中实现 剪枝 技术的措施 2 4 为了便于处理 设置一个缓冲区 该缓冲区内仅存放相关节点的后续节点 5 首先将0放入其中 6 当节点i进入缓冲区时 其对应的值a i 被初始化 然后a i 的值在随后的阶段中被更新 当节点i离开缓冲区时 其值即为该点真正的值 并用于其他节点 依赖于此节点 的计算 7 节点i的后续节点是否要计算 还取决于i是否为相关节点 若不是 则不再计算其后续的其他节点 具体实现过程 Biocomputingtechnology Multiplesequencealignment 设节点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星形比对 Biocomputingtechnology Multiplesequencealignment 启发式方法作为首选 启发式方法不一定保证最终能得到最优解 但在大多数情况下 其计算结果接近于最优结果 启发式这类方法能够大大减少所需的计算时间 加快计算速度 渐进法是多重序列比对中常用到的启发式方法 其基本思想是将序列多重比对转化为序列两两比对 逐渐将两两比对组合起来 最终形成完整的多序列比对 组合两两序列比对的方法有 星形结构或者树形结构 1 星形比对的基本思想 星形比对是一种启发式方法 首先由Gusfield提出 在给定的若干序列中 选择一个核心序列 通过该序列与其它序列的两两比对形成所有序列的多重比对 从而使得 在核心序列和任何一个其它序列方向的投影是最优的两两比对 Biocomputingtechnology Multiplesequencealignment 2 星形比对算法概述 1 Biocomputingtechnology Multiplesequencealignment 设s1 s2 sk是k条待比对的序列 假设已知核心序列是sc 1 c k 则可以利用标准的动态规划算法求出所有sc和si的最优两两比对 这个工作的时间复杂度为O kn2 假设所有序列的长度约为n 将这些序列的两两比对聚集起来 并采用 只要是空位 则永远是空位 的原则 聚集过程从某一个两两比对开始 然后逐步加上其他的两两比对 在这个过程中 逐步增加sc中的空位字符 以适应其他的比对 但决不删除sc中已存在的空位字符 2 星形比对算法概述 2 Biocomputingtechnology Multiplesequencealignment 随着后续比对的不断加入 一方面我们有一个由sc指导的 已经建立好的部分序列的多重比对 另一方面我们有sc与一个新序列之间的比对 在两种比对中我们插入尽可能少而必要的空位 以保持与扩展的sc序列相匹配 然后 将新扩展的序列加入序列群中 且它和其它扩展的序列具有相同的长度 这个过程一直进行到所有的两两比对都加入以后结束 这一步所需的计算量为O k2n 算法总的时间复杂度为O kn2 k2n scs1s2 sk sc s1 sc s2 sc sk 两两比对 多重比对 Biocomputingtechnology Multiplesequencealignment 方法1 尝试将每一个序列分别作为核心序列 进行星形多重序列比对 取比对结果最好的一个 即SP得分最高的为最好 方法2 是计算所有的两两比对 取下式值最大的一个 3 如何选择核心序列 Biocomputingtechnology Multiplesequencealignment 4 举例 有5个序列 s1 ATTGCCATTs2 ATGGCCATTs3 ATCCAATTTTs4 ATCTTCTTs5 ACTGACC sc s1 ATTGCCATTATTGCCATT ATTGCCATTATTGCCATTATGGCCATTATC CAATTTTATCTTC TTACTGACC s1 s2 s1 s3 s1 s4 s1 s5 ATTGCCATT ATGGCCATT ATC CAATTTTATCTTC TT ACTGACC Biocomputingtechnology Multiplesequencealignment 星形比对是一种多重序列比对近似的方法 然而是一种比较好的近似方法 如果用代价来评判多重序列的比对结果 则可以证明 用该方法所得到多重序列比对的代价不会大于最优多重序列比对代价的两倍 Biocomputingtechnology Multiplesequencealignment 4 2 5树形比对 k个待比对的序列 具有k个叶节点的树每个叶节点对应一条序列将序列赋予树的内部节点 可以计算树中每个分支的权值 权值代表对应分支连接的两个序列之间的相似性 所有权值的和就是这棵树的得分树形比对的问题 寻找一种树的内部节点序列赋予方式 使得树的得分最大 星型比对可以看作是树形比对的一种特例 Biocomputingtechnology Multiplesequencealignment 将CT CG CT分别赋予节点x y z 则树的得分为8 CT CG CT 多重序列比对 两两序列比对 合并两个比对 比对的比对 AligmentofAlignment AA算法 打分函数 P a a 1P a b 0 a b P a P b 1 1 1 1 1 2 2 G TCATC GCTG 比对结果 Biocomputingtechnology Multiplesequencealignment AA算法概述 1 Biocomputingtechnology Multiplesequencealignment 假设有两个多重比对 1和 2 1代表序列s1 s2 si的多重比对 2代表序列t1 t2 tj的多重比对并且 s1 s2 si t1 t2 tj 代表s1和t1的两两比对 则计算与 相一致的 1和 2比对的算法如下 AA算法概述 2 Biocomputingtechnology Multiplesequencealignment 标定 1的各列 如果s1在比对中对应位置的编辑操作不是插入或删除 则这些列分别标记为s1对应位置上的字符 a1 a2 al1 s1 l1 标定 2的各列 如果t1在比对中对应位置的编辑操作不是插入或删除 则这些列分别标记为t1对应位置上的字符 b1 b2 bl2 t1 l2对a1 a2 al1和b1 b2 bl2进行比对 在所得到的比对中 对于 1 2和 中原来有插入或删除操作的位置 恢复其原有的实际字符或空位字符 Biocomputingtechnology Multiplesequencealignment a1a2a3a4 b1b2b3b4b5 对于n个序列的树形比对的基本算法过程如下 1 初始化 对于每个序列 生成一个叶节点 2 利用AA算法合并两个节点 形成一个新节点 合并的结果放在新节点中 原来的两个节点作为新节点的子节点 3 反复执行 2 直到形成n个叶节点的树根为止 根节点中的序列即为最终的多重比对结果 s1s2s3s4 1 2 Biocomputingtechnology Multiplesequencealignment 4 2 6CLUSTAL算法 Biocomputingtechnology Multiplesequencealignment CLUSTAL算法是一种渐进的比对方法 它是由D G Higgins和P M Sharp于1988年首次提出的 目的 对来自不同物种的功能相同或相似的蛋白进行比对和聚类分析 可对这些物种的亲缘关系进行判断 并且对这些蛋白在生物进化过程中的保守性作出估计 Clustal算法包括三个阶段 1 先将多重序列进行两两比对 基于这些比对 计算得到一个距离矩阵 该矩阵反映每对序列之间的关系 2 根据距离矩阵计算产生系统发育树 以此来确定被比较序列间相似的程度3 有了这棵系统发育树的指导 从关系最紧密的两条序列开始 逐步引入邻近的序列或序列组 并不断重新构建比对 直到所有序列都被加入为止 Biocomputingtechnology Multiplesequencealignment 举例 Biocomputingtechnology Multiplesequencealignment 已知三级结构的七个球蛋白序列 分别为 Hbb Human 人的 球蛋白Hbb Horse 马的 球蛋白Hba Human 人的 球蛋白Hba Horse 马的 球蛋白Myg phyca 抹香鲸的血红蛋白Glb5 petma 七鳃鳗蓝血红蛋白Lgb2 Luplu 羽扇豆的豆血红蛋白 第一阶段 两两比对产生距离矩阵 Biocomputingtechnology Multiplesequencealignment 用完全动态规划法计算的7个球蛋白序列之间的7 7的距离矩阵 第二阶段 产生指导树 Biocomputingtechnology Multiplesequencealignment 根据第一阶段结果中的距离矩阵 用邻近归并法来计算产生一棵有分枝长度和序列权重的有根树 第三阶段 渐进的比对 Biocomputingtechnology Multiplesequencealignment 这一阶段基本的步骤是通过一系列两两比对来构建更大的比对序列组 按照指导树中的分支顺序 从有根树的末梢到树根逐渐进行 根据图4 22的有根树 通过下面的顺序对序列进行比对 1 对 2 3 对 4 8 对 9 5 对 10 6 对 11 7 对 12 每一阶段使用了有残基权重矩阵和空位开放及空位扩展罚分的完全动态规划算法 每一阶段都由对两个已经存在的排列或序列进行比对组成 在先前比对中出现的空位仍然是固定的 图4 23 Biocomputingtechnology Multiplesequencealignment 序列权重的作用 计算多重序列比对得分 Biocomputingtechnology Multiplesequencealignment peeksavtalgeekaavlalpadktnvkaaaadktnvkaaegewqlvlhvaaektkirsa 多重序列比对的统计特征分析 对于所得到的多重序列比对 我们往往需要进行归纳分析 总结这些序列的特征 或者给出这些序列共性的表示 表示方式1 保守序列表示方式表示出序列每个位置上最可能出现的字符或所有可能出现的字符表示方式2 特征统计图谱 profile 表示方式 或特征统计矩阵表示方法 Biocomputingtechnology Multiplesequencealignment 令P P1 P2 Pj PL P表示在多重比对 的每一列上各种字符出现的概率分布Pj pj0 pj1 pj A A代表字母表 Pjk代表字母表A中第k个字符在第j列出现的概率 pj0 第0个字符是特殊的空位符号 称P为多重序列比对 的特征统计图谱或特征统计矩阵 P反映了多重序列比对 各个列的特征 表示方法2 特征统计图谱 Biocomputingtechnology Multiplesequencealignment 12345 位置 A0 80 20 20 60 0T0 00 40 60 41 0C0 20 20 20 00 0G0 00 20 00 00 0 碱基 Biocomputingtechnology Multiplesequencealignment ATTATAACTTCTTATACTTTAGAAT 利用保守序列或者特征统计图可以判断一个序列是否满足一定的特征给定一个序列s a1a2 am 定义字符a在第j位的代价为其中 A 代表字母表A的长度 Ak代表A的第k个字符 A0代表空缺字符 整个序列s的代价为 一条序列与特征统计图相对照 如果代价值小 说明该序列具有相应的特征 否则该序列不具备相应的特征 Biocomputingtechnology Multiplesequencealignment 4 2 7隐马尔可夫模型HMM Biocomputingtechnology Multiplesequencealignment HMM HiddenMarkovModel 1 马尔可夫链 复习 具有马尔可夫性的离散状态随机过程称为马尔可夫链 马尔可夫性 即无后效性 若已知现在的状态 将来与过去无关 离散状态的马尔可夫链的定义 考虑只取有限个或可数个状态的随机过程 Xn n 0 1 2 假设对一切状态 i1 i2 in 1 i j 一切n 0 有 P Xn 1 j Xn i X0 i0 P Xn 1 j Xn i 成立 则称随机过程 Xn n 0 1 2 为离散状态的马尔可夫链 1 马尔可夫链 复习 Biocomputingtechnology Multiplesequencealignment P Xn 1 j Xn i 称为马尔可夫链的一步转移概率 记为Pij n n 1 若转移概率Pij n n 1 不依赖n 则可简记为Pij 则称此马尔可夫链是时齐的 否则是非时齐的 用P记转移概率所对应的矩阵 即转移矩阵 transitionmatrix 一个马尔可夫链的概率分布完全由它的转移矩阵与初始分布决定 2 隐马尔可夫模型HMM Biocomputingtechnology Multiplesequencealignment 隐马尔可夫模型是一种概率模型 它是由马尔可夫链发展扩充而来的一种随机模型 这种方法已经成功地应用于多个领域 如语音识别 光学字符识别等 HMM在生物信息学领域中也有着重要应用 如 序列分析 基因识别等 HMM可以被理解为一个双重随机过程 1 系统状态变化的随机过程 2 由状态决定输出的随机过程 一个隐马尔可夫链 Xt Yt 包含两部分 一个潜在的 不可观察的有限状态马尔可夫链 Xt 一个外显的 可观察的随机过程 Yt Yt的分布依赖于Xt 2 隐马尔可夫模型HMM Biocomputingtechnology Multiplesequencealignment HMM是用概率统计的方法来描述时变信号的过程 HMM是一个动态的统计模型 可以用有限状态机FSM来描述 FSM finitestatemachine 有限状态机可以从一种状态转移到另一种状态 每个状态转换有不同的概率 某状态是否转移到下一状态取决于该状态的状态转移概率 而在某一状态下能看到哪一个观测值 取决于该状态的观测概率 隐马尔可夫模型HMM Biocomputingtechnology Multiplesequencealignment 记HMM模型为 其中 A为状态转移概率矩阵 B为观察概率矩阵 为初始状态分布 HMM模型的建模工作主要为确定这三个参数 HMM的三个经典问题 Biocomputingtechnology Multiplesequencealignment 问题1 评测问题 evaluation 已知模型和输出序列O 求由生成O的概率 问题2 译解问题 decoding 已知模型和输出序列O 求最有可能生成O的状态转移序列 问题3 学习问题 learning 已知模型和输出序列O 求最有可能生成O时模型的参数 Profile 概形 谱 Biocomputingtechnology Multiplesequencealignment 概形是对一组序列进行全局多重比对时被发现的 是将比对中具有较高保守区域变成小的多重比对 然后得到多重比对记分矩阵 概形由更像小的多重排列的列构成 可以包括 匹配 失配 插入 缺失 概形一旦生成 就可用于寻找一个可能与之匹配的目标序列 它利用表中记分来评价每个位置的可能性 例 25个氨基酸长的概形表格 有25列 每列将有20个记分值 一列中每个匹配氨基酸的记分都在概形中对应的位置上 缺点 偏向性 ProfileHMM 1 模型结构 Biocomputingtechnology Multiplesequencealignment 对于生物序列而言 HMM的字符当然是20个字母的氨基酸或4个字母的核苷酸 但依据不同的问题 其它的一些字符也可以使用 如64个密码子的三联体字母 3个字母 coil 的二级结构等 编码蛋白质的原始DNA序列 在生物的进化过程中会受到自然环境和各种因素的影响 使翻译出的蛋白质序列经历突变 遗失或引入外源序列等变化 最后按不同的进化路径分化 形成多种功能相近的蛋白质 所以 可以把这些蛋白质看作由一个基本蛋白质序列经过插入 删除或替换了某些氨基酸残基而形成的 这个过程可以用HMM来表示 图4 9 Biocomputingtechnology Multiplesequencealignment 图中 方形代表匹配状态 主状态 即输出的氨基酸和基本序列中对应的氨基酸匹配 圆形表示删除或缺失状态 即从原始蛋白质序列中去掉一个特定的氨基酸 菱形表示氨基酸的插入 即在原始蛋白质序列插入一个氨基酸 各状态间的箭头表示状态间的转换途径 注意 不同的参数会使模型以不同的概率产生新的氨基酸 一个训练好的模型可以代表有共同特征的蛋白质序列 图4 10 Biocomputingtechnology Multiplesequencealignment ProfileHMM与标准的Profile的比较 Biocomputingtechnology Multiplesequencealignment ProfileHMM有正规的概率作基础 对于序列的删除和插入状态的记分也有较为可靠的理论依据 而标准的Profile纯粹是一种启发式的方法 HMM用统计方法估计序列某一位点核苷酸或氨基酸残基出现的真正概率 而标准的Profile却是用自身的观察频率给核苷酸或氨基酸残基指派分数 由于 ProfileHMM方法从10至20个核苷酸序列构成的比对中提取的信息 相当于用标准的Profile从40至50个核苷酸序列构成的比对中提取的信息 2 用HMM给序列打分 Biocomputingtechnology Multiplesequencealignment 因训练好的HMM代表一个蛋白质家族 我们可以用它对序列进行打分 根据分值来衡量这条序列是否属于由此HMM代表的蛋白质家族 得分高于阈值 证明这条序列是这个家族的成员 得分低于阈值 说明它不是家族的成员 HMM用于分析蛋白质序列的原理 Biocomputingtechnology Multiplesequencealignment 分析蛋白质产生不同序列的概率 对于与模型相符合的序列 能以较大的概率产生该序列 若不与该模型符合的序例 则按此模型产生该序列的概率会较小 任何一条序列都可以由HMM中的一条路径所呈现 对于给定的模型 计算产生某条序列的概率就是计算这条路径上的所有输出概率与转移概率的乘积 图4 11 Biocomputingtechnology Multiplesequencealignment Biocomputingtechnology Multiplesequencealignment 如果已知确切的状态路径 这个计算是很简单的 在一个实际模型中 可以通过不同的状态路径产生同一条序列 因此 产生一条序列的正确的概率是所有可能的路径产生该序列的概率的总和 这种计算复杂并且难以实施 除非是一条很短的序列 以下介绍两种很好的替代方法 Viterbi算法 可以计算出由模型产生某序列的最有可能的路径 前向算法 forwardalgorithm 用于计算所有路径的概率和 Viterbi算法 Biocomputingtechnology Multiplesequencealignment 作为一种优化算法 动态规划算法在生物信息处理中得到广泛的应用 前期课程 我们已经介绍了动态规划算法在序列比较中的应用 这里 介绍应用动态规划算法求解HMM模型的最优路径的方法 这个方法就是著名的Viterbi算法 图4 12 D1 D2 D3 I0 I1 I2 I3 M1 M2 M3 0 06 Biocomputingtechnology Multiplesequencealignment 最优路径 开始I0M1M2M3结束在这条路径上产生 ACCY 的概率为 5 9728 10 5 图4 13 0 0015 0 0046 0 4850 0 2231 Biocomputingtechnology Multiplesequencealignment 前向算法forwardalgorithm Biocomputingtechnology Multiplesequencealignment 给定一个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 Biocomputingtechnology Multiplesequencealignment 图4 12的HHM产生 ACCY 序列的概率为 3 3849 10 7 6 8965 10 5 6 9303 10 5 局部和全局打分 Biocomputingtechnology Multiplesequencealignment 有时用局部打分会更合理 即序列的分值由子序列的最高分值所替代 例 图4 11中ACCY序列 将概率转换为对数形式 序列的全局打分就是以下四个分值之和 13 25 Score ln 0 4 ln 0 3 2 12Score ln 0 46 ln 0 6 1 29Score ln 0 97 ln 0 5 0 72Score ln 0 015 ln 0 73 ln 0 01 9 12 2 12 1 29 0 72 9 12 13 25 分值太低了 以至于不能判定ACCY是否为被模拟的蛋白质家族的成员 Biocomputingtechnology Multiplesequencealignment 假设ACCY是图4 15中所示蛋白质家族的成员 在这种情况下 全局分就不是一个很好的衡量标准 如果采用局部打分 最后分值就会高得多 我们发现最高分的子序列为CC 分值为 2 01 这个分值就足以判定ACCY为此家族的成员 这时 局部打分比全局打分更精确 3 模型的训练 Biocomputingtechnology Multiplesequencealignment 根据训练用蛋白质序列的平均长度确定模型中各个状态的数量而建立HMM 选用训练算法对HMM进行训练 即调整该模型的参数 转移概率和输出概率 使得该模型适用于训练所用的序列 并能够以最大的可能性产生参与训练的观察序列 3 模型的训练 最大似然法ML Biocomputingtechnology Multiplesequencealignment 给定在训练数据集中的一组序列 s 1 s 2 s 一个模型究竟在多大程度上适合该数据集 可由根据该模型观测到每一个序列的联合概率来表征 其中 是模型产生第j个序列的概率 上式称为该模型的似然 最大似然 maximumlikelihood ML 其原理可用于测量模型的性能 即选取模型参数使上式值最大 3 模型的训练 最大后验概率MAP Biocomputingtechnology Multiplesequencealignment 基本思想 在给定序列下 使模型的后验概率最大 假设 对所有可能的模型参数存在优先概率分布 这种概率分布体现模型的特点 为模型优先概率分布 为归一化常数 其中 3 模型的训练 Baum Welch算法 Biocomputingtechnology Multiplesequencealignment 方法 Baum Welch算法 目的 给定观测值序列 通过计算确定一个模型 使得P O 最大 Baum Welch算法的基本思路 初始模型 0 待训练模型 基于 0以及观察值序列 训练新模型 如果logP X logP X 0 说明训练已经达到预期效果 算法结束 否则 令 0 继续第 步工作 4 序列加权 Biocomputingtechnology Multiplesequencealignment 如果训练集中有一小组序列具有高度相似性 那么训练出来的HMM就会过于专门化 如果HMM过于专门化 那么在对其它序列进行打分时 有可能将那些本该属于此家族的序列打分过低 而据此判断其不属于本家族蛋白质 出现偏差 为防止这种情况发生 在估计模型参数时要考虑序列加权问题 序列加权方法有2种 基于树结构的加权法 位置特定的加权法 基于树结构的加权方法 Biocomputingtechnology Multiplesequencealignment 原理 一个家族中的序列被放置在树的各个分枝上 代表它们从同一个祖先分化而来 形象化地将树用导线制成 树根电压为 树叶电压为0 流经每根树枝的电流用基尔霍夫定律 Kirchoff sLaw 来计算 这些电流就被当作序列的权重 图4 16 Biocomputingtechnology Multiplesequencealignment 图4 174 18 Biocomputingtechnology Multiplesequencealignment 位置特定的加权方法 Biocomputingtechnology Multiplesequencealignment 原理 基于多序列比对中每一列观察到的变化 序列中每一个位置的权重可用下式计算 其中 m 该列中不同的氨基酸残基数目 k 感兴趣的氨基酸残基在该列中出现的次数 整条序列的权重为所有位置的权重的平均值 将各条序列的权重规范化 使其和为1 4 18 Biocomputingtechnology Multiplesequencealignment 5 过分拟合 overfitting 与一般化 generalization Biocomputingtechnology Multiplesequencealignment 能精确估计氨基酸分布的方法对于建立好的HMM是非常必要的 但也有潜在的缺陷 如图4 18所示 假设一个训练集由这四条序列组成 图中第一列只存在一种氨基酸C 若用前述的方法 则其他19种氨基酸在这个位点出现的概率都为0 如果我们建模时将这些概率都设为0 那么训练出来的模型就不可能识别以其他氨基酸开头的家族成员 这就称为过分拟合 解决办法 一般化的方法可以解决这个问题 方法 利用假记数 pseudocount 假记数 pseudocount Biocomputingtechnology Multiplesequencealignment 即使在训练集比对的某列未出现给定的氨基酸 我们也给它分配一个虚假计数 在计算概率时 对假计数的处理与对实际观察的计数的处理完全相同 假记数 pseudocount Biocomputingtechnology Multiplesequencealignment Biocomputingtechnology Multiplesequencealignment 思考题 1 1 多重序列比对的定义 2 多重序列比对分析的目的是什么 3 建立SP模型的目的是什么 4 在具体计算SP得分时 有两种方案 它们是怎样进行的 这两种方案等价的条件是什么 5 在多重序列比对中 允许某两条序列比对中出现空位与空位的比对 为什么 在计算比对得分时 怎样处理空位与空位的比对 6 如果利用标准的动态规划算法对10条序列进行比对 那么 在计算超晶格某个节点的得分值时 需要考虑多少个前趋节点 Biocomputingtechnology Multiplesequencealignment 思考题 2 7 如果给你一个多重比对 如何将这个比对投影到特定的两条序列上 要求画出投影路径 8 简述星形多重序列比对的基本思想 其中 核心序列Sc的选择方法 9 简述Carrillo Lipman的优化计算方法的基本思想 10 已知打分函数 p a a 1 p a b 1 p a p b 2 求序列s ACTG t ATCT经过特定断点的序列比对 11 简述树形多重序列比对的基本思想 其中 AA算法原理是什么 12 简述CLUSTAL算法原理 Biocomputingtechnology Multiplesequencealignment 思考题 3 13 多重序列比对的特征统计矩阵的概念14 HMM模型的基本原理15 HMM的三个经典问题是什么 这三个问题用什么算法解决 算法的原理是什么 图4 194 20 Biocomputingtechnology Multiplesequencealignment Biocomputingtechnology Multiplesequencealignment