生物信息学常用算法简介.ppt

生物信息学常用算法简介 北京大学生物信息中心北京基因组研究所李松岗lsg 010 62756803 常用算法1 动态规划 动态规划算法是一种优化算法 它本质上是一种有效的穷举法 它的基本想法是最优路径上的每一段都应该是局部的最优路径 动态规划算法的典型应用 序列比对 序列比对应用举例 序列组装进化分析保守区发现蛋白质结构与功能预测cDNA的基因组定位基因结构与功能分析 序列比对模型 类型 全局比对与局部比对需考虑的因素 替换 插入 删除例 AGCTA CGTACATACCAGCTAGCGTA TAGC打分系统 替换矩阵 记为 a b 其中a b为我们考虑的字符集中的元素 比对算法的目标 就是找到在给定打分系统下 得分最高的比对方式 动态规划算法 全局比对 两序列 A a1a2a3 amB b1b2b3 bn用Ai Bj分别表示上述序列的前i个和前j个碱基 矩阵元素S i j 表示Ai Bj所有可能比对中的最高得分 则有递推公式 S i j max S i 1 j 1 ai bj S i j 1 b S i 1 j a 局部动态规划 递推公式改为 S i j max 0 S i 1 j 1 ai bj S i j 1 b S i 1 j a 局部动态规划图示 动态规划算法的改进 用动态规划方法进行序列比对 需要nm到nm2的计算时间和nm的存储空间 当序列很长时 常常无法计算 因此人们陆续提出了许多改进算法 能节省空间和时间 有兴趣的同学可参考相关文献 其他DNA打分矩阵及其对比对结果的影响 例如 若得分大于罚分 则可得到长的 有较多插入删除的结果 反之 则得到短的 局部的比对结果 蛋白质序列比对的打分矩阵 PAM矩阵 PersentAcceptedMutation 基于进化模型的打分矩阵 当进化过程中一条序列1 的氨基酸发生了突变 定义该序列在进化的历史上走过了1个PAM单位 此时定义的转移矩阵称为1 PAM的突变矩阵 Dayhoff等 1978 从71个蛋白家族中的1300条近相关 closelyrelated 序列出发 其中任何两对序列之间氨基酸残基差异不大于15 通过构造进化树对序列进行联配 得到氨基酸对之间的联合概率分布 在此基础上得到了1 PAM的突变矩阵 ARNDCQEGHILKMFPSTWYVA989055612911125256921029141217R4990752216438123020355432N349888182856131110112138131D4221990507285600500375010C31109946001111011031112Q411751985618214131461455112E85630028989027111530475113G1149724399523013102102211H14731931989511322122191I21202210298782222671152242L542038213359919348224345519K533135022152822988351469122M3110241021012298595023114F10103100451009992301110283P62220531212301994365012S23517899786127418986232242T115116485276297273398791212W01001100101013000995640Y12213110131212221211099242V1521083412511431252314149884 1 PAMmutationmatrix allvaluesmultipliedby10000 表中各列满足若fi i 1 20 表示20种氨基酸在自然界中的分布 该矩阵还满足 由于fi是自然界中氨基酸经过长期进化后形成的一种稳定分布 因此满足关系也就是说 可以通过对1 PAM突变矩阵外推得到n PAM的突变矩阵 用来表示相距n PAM进化单位的蛋白质之间氨基酸残基的突变概率 即 ARNDCQEGHILKMFPSTWYVA13506777337338847487818846476495977146714618341006899342468803R4601583568493323751589424616304320995361253429512507366351337N4254851092747299534563496601239228552275225366558507197326271D50149688115932606661007549591221213604268188453597536164284273C2121141249126601079411913814513410015315793193166149169187Q359530442468215705558300496243261538303208359390379201253264E58064572210952938631334483635315305768370237525614572215312376G8275838017524665856083387530264264567332224512799569291292346H194271310259173309256170946143153267175231182225221194387149I4803303042394463753142083531460109435410277263203795103814891185L71756647337467365349233961217792390576179815015615757077939421443K53510977136623118407744546673583591240425276510600605264360394M195154138114185183144103170403435165608326130166198181217330F23919320214334022416512439951064919158320411712132469371312419P4843673663852254344103183532512713962591912614495467144222299S774580741671619625636657578394368616439315656997844285393475T7075876876165446216054785805424626355353726328621101274397622W579057401037048521088611059104301416258339833672Y19521523517329322017513053927632720031110541602132118431955253V7084374133526894864463264401415106046410047124545456983825351501 250 PAMmutationmatrix allvaluesmultipliedby10000 对250 PAM突变矩阵 有 即经过250个PAM单位进化后的蛋白质分子 与它的祖先相比较 大约只有20 左右的氨基酸残基保持不变 当我们通过动态规划对两个序列进行联配时 用到PAM突变矩阵的另一种形式 PAM打分矩阵 其中PAM 1打分矩阵定义为 PAM 250打分矩阵定义为 C17S 1912T 22 1312P 33 19 2013A 18 14 18 1911G 25 19 25 25 1811N 24 16 18 24 22 1913D 32 19 20 23 21 21 1413E 35 19 21 22 19 24 20 1312Q 29 18 19 20 19 23 17 19 1314H 22 20 20 23 22 24 15 19 19 1416R 24 20 21 23 22 22 20 24 21 15 1813K 33 20 18 22 21 24 17 20 16 14 19 1312M 22 22 19 31 19 28 25 33 24 18 21 24 2116I 26 26 20 28 25 35 26 34 27 26 25 27 24 1312L 25 26 24 24 22 31 28 35 27 21 24 24 24 13 1410V 19 24 17 25 17 30 27 34 23 24 27 25 24 18 11 1712F 22 28 25 30 25 32 27 33 33 26 20 31 30 17 19 16 2214Y 21 22 25 27 26 30 22 26 27 24 14 23 25 23 24 22 23 1215W 22 25 29 32 29 27 28 38 30 24 23 22 31 23 25 23 29 16 1519CSTPAGNDEQHRKMILVFYW PAM 1scoringmatrix 经过四舍五入 PAM 250scoringmatrix 经过四舍五入 C12S02T 213P 3106A 21112G 310 115N 410 1002D 500 10124E 500 100134Q 5 1 100 11224H 3 1 10 1 221136R 40 10 2 30 1 1126K 500 1 1 21001035M 5 2 1 2 1 3 2 3 2 1 2006I 2 10 2 1 3 2 2 2 2 2 2 225L 6 3 2 3 2 4 3 4 3 2 2 3 3426V 2 10 10 1 2 2 2 2 2 2 22424F 4 3 3 5 4 5 4 6 5 5 2 4 5012 19Y0 3 3 5 3 5 2 4 4 40 4 4 2 1 1 2710W 8 2 5 6 6 7 4 7 7 5 32 3 4 5 2 60017CSTPAGNDEQHRKMILVFYW BLOSUM打分矩阵 BLOSUM矩阵建立在BLOCKS数据库的基础上 这个数据库由很多序列块 block 组成 一个block包含了对SWISSPROT库中蛋白序列进行多序列联配 multiplealignment 得到的一组氨基酸序列片段 这些序列具有很高的相似性 把它们截断为长度相同 中间没有插入或删除空位的一组 就称为一个block 通过对这个数据库中联配氨基酸残基出现频率的统计 生成了BLOSUM矩阵 下表为BLOCKSv 9数据库中某一个块的数据结构片段 其中包括三个类 cluster 类中的每一条序列都能在同一类中至少找到另一条序列 两条序列之间80 以上的氨基酸残基是相同的 FA9 BOVIN 6 LEEFVRGNLERECKEEKCSFEEAREVFENTEKTTEFWKQY29FA9 CANFA 45 LEEFVRGNLERECIEEKCSFEEAREVFENTEKTTEFWKQY29FA9 HUMAN 52 LEEFVQGNLERECMEEKCSFEEAREVFENTERTTEFWKQY30MGP BOVIN 56 ERIRELNKPQYELNREACDDFKLCERYAMVYGYNAAYDRY70MGP HUMAN 56 ERIRERSKPVHELNREACDDYRLCERYAMVYGYNAAYNRY73MGP MOUSE 56 KRVQERNKPAYEINREACDDYKLCERYAMVYGYNAAYNRY65MGP RAT 56 ERVRELNKPAQEINREACDDYKLCERYALIYGYNAAYNRY71OSTC BOVIN 57 LGAPAPYPDPLEPKREVCELNPDCDELADHIGFQEAYRRF33OSTC FELCA 6 LGAPAPYPDPLEPKREICELNPDCDELADHIGFQDAYRRF34OSTC HUMAN 57 LGAPVPYPDPLEPRREVCELNPDCDELADHIGFQEAYRRF34OSTC MACFA 6 LGAPAPYPDPLEPKREVCELNPDCDELADHIGFQEAYRRF33OSTC RABIT 6 QGAPAPYPDPLEPKREVCELNPDCDELADQVGLQDAYQRF45OSTC RAT 55 LGAPAPYPDPLEPHREVCELNPNCDELADHIGFQDAYKRI46 一个类被作为一条序列处理 因此如果一个类中包含n条序列 那么其中每一个氨基酸残基出现的频率被乘以权重因子1 n 对BLOCKS数据库进行统计 得到氨基酸对之间的联合概率分布P i j BLOSUM矩阵定义为其中C为常数 在不同的BLOSUM矩阵 BLOSUM62 BLOSUM80等 中 C取不同的值 1 2 1 3 按照62 的标准对BLOCKS数据库进行聚类 即类中的每一条序列都能在同一类中至少找到另一条序列 两条序列之间62 以上的氨基酸残基是相同的 得到BLOSUM62矩阵 按照70 的标准对BLOCKS数据库进行聚类 得到BLOSUM70矩阵 两种打分矩阵的比较 PAM矩阵基于进化的突变模型 且是从进化距离近的数据外推到远的 BLOSUM矩阵基于蛋白家族中的保守区段 不考虑进化距离远近 PAM矩阵计算了相关序列中的所有位点 而BLOSUM矩阵只计算了一些保守区段中的位点 比对得分的统计检验 这种统计检验主要用于局部比对 因为在局部比对中我们需要判断哪些结果是真正有意义的 在局部比对中 通常能比上的只占参加比对序列总数的一小部分 因此从整体上说 可把比对过程视为反复进行大量小片段间的比对 由于随机序列之间的比对结果服从正态分布 故可用极值分布来对某一具体结果进行检验 极值分布 分布函数 故有 常用比对软件介绍 FASTA计算步骤 1 查找查询序列和数据库序列之间长度大于k 氨基酸1 3 核酸1 6 的精确匹配片段 2 把顺序相同 互相间距离小于给定值的匹配片段连成没有洞 gap 的区间 3 用打分矩阵对其中匹配数和匹配密度最高的区间重新打分 选出其中得分最高的区间 4 通过 gap 把得分最高的区间连接起来 加上它们的匹配分数 并减去加入 gap 带来的罚分 记录下分值最高的连接 5 用动态规划重新处理上述有 gap 的联配 BLAST计算步骤 1 屏蔽序列中低复杂性区域 2 对 字 氨基酸 3 tup 核酸 11 tup 在数据库中每一条目的位置建立索引表 3 对查询序列建立字集 对其中每一个字建立它的邻域 即虽然不同 但比对分值大于给定阈值的字 蛋白的每个字大约有50个近邻 4 在数据库中搜索与查询序列字集邻域中相同的字 5 以找到的字为种子 向两端延伸 直到累计得分开始下降 这样就得到了匹配片段 6 对找到的匹配片段进行统计检验 7 对检验显著的片段重做动态规划的局部比对 8 在BLAST2中 不是对每个匹配片段做动态规划比对 而是把数据库每条序列中的匹配片段连起来做有gap的局部比对 这样有可能得到更长的匹配片段 动态规划算法的其它应用 例1 核酸序列的字典模型已知观测序列S s1s2 sn si A C G T 给定模型其中wi称为单词 观测序列是由这些单词生成的 pwi是单词出现的概率 满足求P S M 解 因为即求序列S在模型M下的概率 需要穷举S在模型M下的所有可能划分方式 定义变量Z 1 i 它表示观测序列S中从1到i的片段s1s2 si在模型M下的概率 则得到递归关系其中 i l 表示序列S中从i l 1到i 长度为l的子字串 如果 i l 为模型中的单词 则p i l 等于该单词的概率pw 否则为0 例2 模体的权重矩阵描述模型已知观测序列S s1s2 sn si A C G T 给定模型其中q0为噪声出现的概率 qi i 1 2 6为六个模体出现的概率 满足p bk 0 bk A C G T 是bk作为噪声时出现的概率 满足 p bk i j bk A C G T 是bk在模体i中位置j出现的概率 满足求P S M 解 同样定义变量Z 1 i 在权重矩阵模型下 递归关系关系为 常用算法2 隐马模型 基本概念 随机过程 一族无穷多个 相互有关联的随机变量 记为 由于参数t经常代表时间 故称为随机过程 T常为自然数 整数或区间 当参数取值为整数时 也称为随机序列 马尔可夫过程 取值为整数的随机过程 若t i时刻的取值只与时刻i 1取值有关 则称为马尔可夫过程 隐马氏模型 HMM 模型的数学描述 存在一个隐序列H 它是不可观测的 且由以下参数生成 其中 为初始状态出现概率 T 为转移概率 即t P hi hi 1 属于 为字符集 即隐序列由哪些字符组成 观测的结果称为明序列O 它由隐序列按照生成概率e a生成 其中e a P a a 为明序列字符集 隐马模型的典型应用 基因识别 基因结构及重要信号 基因识别模型 问题 从基因组序列出发 识别其中的基因 基因组注释 输入 基因组序列输出 各种基因元件 包括启动子 外显子 内含子等 在基因组上的位置 如果需要 也可输出翻译的蛋白质序列 预测结果的可靠程度等 例 基因组编码区的隐马模型 设基因组由两种功能区域组成 即编码区和非编码区 分别由字母c n代表 转移矩阵为同种字母延伸或变为另一种字母的概率 初始状态概率为第一个字母出现c或n的概率 明序列由A C G T四个字母组成 生成概率分别为编码区和非编码区四个字母出现的概率 隐马氏模型的三种典型问题 可能性问题 给定模型参数 当观察到一个明序列时 这一明序列确实由给定模型生成的概率有多大 解码问题 给定模型参数 当观察到一个明序列时 这一明序列所对应的最可能的隐序列是什么 学习问题 观察到足够多明序列时 如何估计转移概率和生成概率 基本算法 可能性问题 已知参数为w的隐马模型和长为T的序列O o1o2 ot oT 求在给定模型下观测到序列O的概率 解 对于任意一条隐序列H 有 其中 为隐序列H第i 1位置的状态 为i位置的状态 当i 1时 t 应由 代替 因此 所求概率为 其中H为一切可能的隐序列 由于H的数量随T的增加以指数增长 这种算法实际并不可行 前向算法 定义在模型w下 隐序列第i位置状态为 明序列为o1 oi时的概率为 则有递推公式 初值为 则有 后向算法 定义在模型w 隐序列第i位置状态为 条件下 明序列为oi 1 oT时的概率为 则有递推公式 初值为 则有 使用前向或后向算法 计算P O w 的复杂度约为O N2 级 联合应用前向和后向算法 可以容易地计算隐序列位置i处于状态 的概率 通过计算上述概率的极大值 可以得到位置i处隐序列最可能的状态 但在解码问题中 我们更关心的是最可能的隐序列 它是一个整体 不能分解开一位一位算 解决这个问题 需要使用Viterbi算法 Viterbi算法 定义其中Hi为一切长为i 以 结束的隐序列 显然 i 对应的隐序列就是最可能生成前i个字母明序列的隐序列 递推公式为 显然就是解码问题的解 学习算法 HMM的学习算法有许多种 包括EM 期望最大化 算法 不同形式的梯度下降法 模拟退火算法等 下面我们主要介绍EM算法 常用算法3 EM算法 EM算法的用途 EM算法 Theexpectationmaximizationalgorithm 是根据不完整数据作最大似然估计的一般方法 这种不完整可以是缺失了某些数据 也可以是存在某些无法观测的隐藏变量 似然函数的概念 我们从某个未知分布得到了一组观察值X x1 x2 xn 此未知分布是由一组参数 1 2 m 决定 定义似然函数为 由于对数函数是单调递增的 因此不会改变极值点的位置 而且它可以把连乘变成连加 从而大大简化计算 因此实际工作中常使用的是对数似然函数 即 最大似然估计 目标 根据已知观测数据估计总体参数 似然函数可视为以未知参数为自变量的函数 它的统计意义又是在未知参数下观测到该组数据的概率 因此我们很自然地取使似然函数达到最大值的未知参数为估计值 这就是最大似然估计 即 为简化计算 实际工作中常使用的是对数似然估计 即 具体计算公式为 如果参数不止一个 相应令偏导数为0即可 例题1正态分布的最大似然估计 设x1 x2 xn是取自正态总体N 2 的简单随机子样 与 2是未知参数 求 和 的极大似然估计 解 由于故有似然函数 取对数 有 似然方程为 由 1 解得 代入 2 得 即 和 2的极大似然估计分别为和 例题2二项分布的最大似然估计 取n粒种子作发芽试验 其中有m粒发芽 求发芽率p的最大似然估计 解 每粒种子发芽与否可视为两点分布 发芽 则X 1 其概率为p不发芽 则X 0 其概率为1 p由似然函数的构造 有 L p P X x1 p P X x2 p P X xn p 由于共有m粒发芽 n m 粒不发芽 L p pm 1 p n m令上式等于0 由于 有 即 发芽率p的极大似然估计为 EM算法的一般概念 X 观测得到的数据 它是不完全的 Z X Y 是完全的 其中Y是缺失的或隐藏的 它的联合分布为 对于这个联合分布函数 我们可以定义它的似然函数 称为完整数据的似然函数 而X的似然函数称为非完整数据的似然函数 注意完整数据似然函数的自变量中 X和 为常数 而Y是随机变量 且服从于某种由X和 所决定的分布 参数X是观测数据 是完全确定的 参数 则是我们需要估计的 在计算过程中它会不断调整 设 i 1 为当前我们使用的 估计值 则我们可以定义完整数据对数似然的期望 函数Q i 1 就是EM算法中 E步骤所要计算的值 在上述函数中 i 1 是上一步计算的结果 是确定的值 而 是接着的M步骤中要调整的参数 调整它以便使似然函数的值增加 即 在EM算法中 上述E步骤和M步骤不断重复 每一次重复似然值都会增加 直到收敛到似然函数的局部极大值 隐马模型的参数估计 利用EM算法进行隐马模型的参数估计 又称为Baum Welch算法 它要解决的主要问题 是求 隐马模型复习 隐马模型 存在隐序列H 由下列参数生成 观测序列O称为明序列 由下列参数生成 直接估计隐马模型参数的算法 定义后验转移概率 注意利用上述表达式可以方便地计算一些重要的后验期望值 隐序列中状态出现的期望数 隐序列中状态转移为的期望次数 按EM算法 我们更关心如何根据观测数据和现有参数估计新的参数 直观上 可以简单地利用上述期望值计算 其中 利用Q函数的估计公式 改用隐马模型的有关参数 有关公式可写为 非完整数据的似然函数 完整数据的似然函数 Q函数 其中为长度为T的所有隐序列的集合 在已知模型参数且针对特定隐序列的条件下 计算观察到明序列的概率是简单的 代入Q函数表达式 得 由于我们要优化的参数分成了相互独立的三部分 故我们可以分别对它们进行优化 对上述表达式中的第一部分 有 由于要在约束下对上式进行优化 引入Lagrange乘子 得 化简 得 由约束条件 得 即 对Q函数表达式第二部分 类似地有 同样引入Lagrange乘子 在约束下 可得 对Q函数表达式第三部分 类似地有 同样引入Lagrange乘子 在约束下 可得 其中比较可知 用Q函数方法求得的叠代公式与直观方法得到的是完全一样的 隐马模型例题1 例1 模体的字典模型已知观测序列S s1s2 sn si A C G T 给定模型其中wi称为单词 观测序列是由这些单词生成的 pwi是单词出现的概率 满足 1 求P S M 2 求路径 使得 解 1 因为即求序列S在模型M下的概率 需要穷举S在模型M下的所有可能划分方式 采用动态规划 定义变量Z 1 i 它