（概率论与数理统计专业论文）生物序列比对.pdf

生物序列比对 中文摘要 生物序列比对 中文摘要 随着人类基因组计划( h g p ) 的实施，生物信息学成为当今最重要的学科之一 生物信息学可以看作是集统计学，数学和计算科学于一体的学科主要用来分析 生物序列，即：d n a 、r n a 和氨基酸( 蛋白质) 序列序列比对是生物信息学中最 基本、最重要的方法通过序列比对，可以发现生物序列的功能、结构和进化信 息，因而序列比对是生物信息学的基础 在本文的第一章，首先介绍了序列比 对方法和得分方案序列比对包括全局序列比对和局部序列比对从比对序 列的数目来划分，又包括两两序列比对和多重序列比对接着，我们给出了一些 序列比对的算法第二章，介绍了b l a s t 算法( 基本的局部比对搜索工具) ，这个 算法是建立在统计理论基础上的，并从两个方面出发，证明了b l a s t 算法第三 章，给出了隐马尔可夫模型，介绍了怎样从多条序列建立隐马尔可夫模型，并用 隐马尔可夫模型来进行多重序列比对第四章，我们讨论了生物网络中的局部图 比对和m o t i f 搜索在得分函数的基础上，我们为拓扑m o t i f s 发展一个搜索算法， 该算法称为图比对，这是一个与序列比对有些相似的方法 在文章最后，对b l a s t 算法、隐马尔可夫模型( 删) 和图比对作了讨论，分 析了算法存在的缺点，并且为将来进一步研究给出了几点建议 关键词：序列比对最大得分片段b l a s t 算法隐马尔可夫模型图比对 作者：黄宁 指导老师：汪四水( 副教授) 生物序列比对英文摘要 b i o l o g i c a lse q u e n c e sa l i g n m e n t a b s t r a c t w i t ht h ed e v e l o p m e n to fh u m a ng e n o m ep r o j e c t ( h g p ) ，b i o i n f o r m a t i e sh a s d e v e l o p e do n eo ft h em o s ti m p o r t a n ts u b j e c t b i o i n f o r m a t i c sc a nb ed e f i n e da sa c o l l e c t i o no fm a t h e m a t i c a l ，s t a t i s t i c a la n dc o m p u t a t i o n a lm e t h o d sf o ra n a l y i n g b i o l o g i c a ls e q u e n c e s ，t h a ti s ，d n a ，r n aa n da m i n oa c i d ( p r o t e i n ) s e q u e n c e s s e q u e n c ea l i g m e n ti st h em o s tb a s i ca n di m p o r t a n tm e t h o d b ys e q u e n c ea l i g n m e n t , w ec a nd i s c o v e rf u n c t i o n a l ，s t r u c t u r a la n de v o l u t i o n a r yi n f o r m a t i o ni nb i o l o g i c a l s e q u e n c e s s o ，s e q u e n c ea l i g m e n ti st h eb a s i co f b i o i n f o r m a t i c s 二 i nt h ec h a p t e r1o ft h i sa r t i c l e ，w ef i r s ti n t r o d u c et h em e t h o d so ft h es e q u e n c e s a l i g n m e n ta n ds c o r i n gs c h e m e s s e q u e n c e sa l i g n m e n tc o n t a i ng l o b a la l i g n m e n ta n d l o c a la l i g n m e n t f r o mt h en u m b e ro ft h es e q u e n c e ，s e q u e n c e sa l i g n m e n tc o n t a i n p a i r - w i s es e q u e n c ea l i g n m e n ta n dm u r i p l ea l i g n m e n t t h e n ，w eg i v es o m em e t h o d s o fs e q u e n c e sa l i g n m e n ta l g o r i t h m i nt h ec h a p t e r2 ，w eg i v et h eb l a s ta l g o r i t h m ( b a s i cl o c a la l i g n m e n ts e a r c ht 0 0 1 ) b l a s ta l g o r i t h mi sb a s e do n s t a t i s t i c s ，a n dw ew i l lp r o v eb l a s ta l g o r i t h mf r o mt w od i f f e r e n tw a y s i nt h ec h a p t e r 3 ，w ed e s c r i b eh m m a n de x p l a i nh o wah m mc a nb eu s e dt op r o d u c eam u l t i p l e a l i g n m e n ta l g o r i t h mf o rn u m b e r so fs e q u e n c e s i nt h ec h a p t e r4 ，w ed i s c u s sl o c a l g r a p ha l i g n m e n ta n dm o t i fs e a r c hi nb i o l o g i c a ln e t w o r k s b a s e d0 1 1t h i ss c o r i n g f u n c t i o n ，w ed e v e l o pas e a r c ha l g o r i t h mf o rt o p o l o g i c a lm o t i f sc a l l e dg r a p ha l i g n m e n t ， ap r o c e d u r e 埘t hs o m ea n a l o g i e st os e q u e n c ea l i g n m e n t i nt h ee n d ，w ed i s c u s st h ep r o b l e mo ft h eb l a s t a l g o r i t h m ，h i d d e nm a r k o v m o d e l sa n dg r a p ha l i g n m e n t a n dt h e nw eg i v es o m es u g g e s t i o n sf o r t h ef u t u r e k e yw o r d s ：s e q u e n c ea l i g n m e n t ，s c o r i n gs c h e m e ，m a x i m a ls e g m e n tp a i r s ， b l a s t a l g o r i t h m ，h i d d e nm a r k o vm o d e l s ，g r a p ha l i g n m e n t w r i t t e n b yh u a n g n i n g s u p e r v i s e db ya s s o c i a t ep r o f w a n gs i s h u i 苏州大学学位论文独创性声明及使用授权声明 学位论文独创性声明 本人郑重声明：所提交的学位论文是本人在导师的指导下，独立 进行研究工作所取得的成果。除文中已经注明引用的内容外，本论文 不含其他个人或集体已经发表或撰写过的研究成果，也不含为获得苏 州大学或其它教育机构的学位证书而使用过的材料。对本文的研究作 出重要贡献的个人和集体，均已在文中以明确方式标明。本人承担本 声g 目的法律责任。 研究生签名： 啦日 学位论文使用授权声明 苏州大学、中国科学技术信息研究所、国家图书馆、清华大学论 文合作部、中国社科院文献信息情报中心有权保留本人所送交学位论 文的复印件和电子文档，可以采用影印、缩印或其他复制手段保存论 文。本人电子文档的内容和纸质论文的内容相一致。除在保密期内的 保密论文外，允许论文被查阅和借阅，可以公布( 包括刊登) 论文的 全部或部分内容。论文的公布( 包括刊登) 授权苏州大学学位办办理。 研究生签名： 啦e l 期：主! 1 2 ：堡曼 导师签名：趁吐日期：犁：置 生物序列比对第一章序列比对 第一章序列比对 1 1 序列比对的背景 生命最重要的物质基础是核酸( d n a 和r n a ) 和蛋白质核酸和蛋白质可以看 成是由字符表中选出的字母组成的，比如：d n a 分子中的碱基分别是腺嘌呤( a ) 、 鸟嘌呤( g ) 、胸腺嘧啶( t ) 、胞嘧啶( c ) 因此，我们一般将d n a 序列字母表示为： a 、c 、g 、t ，同样，蛋白质序列也可以用2 0 个字母表示为： a 、v 、d 、e 、f 、 g 、h 、i 、k 、l 、m 、n 、p 、q 、r 、s 、t 、v 、w 、y ) 随着人类基因组计划的实施，生物信息学成为- i 3 新的学科，成为当今最重 要，最热门的科学发展之一，极大的推动了分子生物学的发展，已经被广泛应用 于基因序列数据的处理等方面，对分子生物学，医学等方面的发展起了巨大的推 动作用 生物信息学最首要的任务就是从大量的生物信息数据中，提取有生物价值的 信息，而序列比对是当前最重要，最基本的方法所谓序列的比对，是指两个或 多个序列按字母比较，尽可能确切地反映它们之间的相似和相异性，于阐明序列 之间的同源关系 进行序列比对的目的之一，就是通过比较，找出序列之间的相似性，发现与 结构相联系的保守序列片段，以及从己知序列预测新序列的结构和功能 1 2 序列比对的方法 序列的比对方法可以按不同的标准进行划分，目前，已知的序列比对方法很 多，包括全局序列比对，局部序列比对根据参与比对序列的数目，可以把序列 比对分为两序列比对和多序列比对 1 2 1 全局序列比对( g l o b a la li g n m e n t ) 全局序列比对：对序列进行全程扫描，考察两条序列整体之间的相似性情况， 在给定得分值的情况下进行序列比较 例1 1 两条d n a 序列，s = c t t a g a ，t = g t a a ，经过序列全局比对，得到下 生物序列比对第一章序y i j l g 对 面三个比对结果 s ：ct ta g as ：ct ta gas ：ct ta g a t ：g ta a t ：gt a at ：一gta a 全局序列比对利用了动态规划的思想，在给定的两条序列全部长度上进行比 对，得到全长序列最优比对例1 1 的两条序列在全局范围内，只有两个字母是 完全匹配的，其他位置都没有好的匹配，说明这两条序列经过全局比对，没有大 的相似性 - 二 1 2 2 局部序列比对( l o c a la l i g n m e n t ) 局部序列比对：当两条序列进行比对时，找出待比对序列中的某一子片段的最优 比对，但是这个最优比对，不一定是全局最优比对的片段 例姚 z ：ee c c 盒暑2 许多蛋白质在全局范围内并不具有相似性，就需要在局部范围里考虑，在大 多数情况下，使用局部比对是较为合理的，这种比对方法可能会揭示一些匹配的 序列段，而本来这些序列段是被一些完全不相关的残基所淹没的因此如果使用 全局比对，很可能会掩埋一些局部的相似性 1 2 3 两两序列比对( p a i r w i s es e q u e n c ea l i g n m e n t ) 两两序列比对：就是把两条未知的序列进行排列，通过字母的匹配( m a t c h ) ，删 除( d e l e t e ) 和插入( i n s e r t ) 操作，使得两条序列达到同样长度，在操作的过程中， 尽可能保持相同的字母对应在同一个位置 例1 2 两条序列s t a c c a g t ，l ：c c c g t a a 的长度虽然一样，若比对不允许 空格，只有一种比对方法 占：t a cc a g t t ：c c c g t a a 但是，在比对中通过加入空格，就可以得到比较好的比对 占：ta c c a gt 一 t ：c cc gta a 2 生物序列比对 第一章序列比对 另外，也可以得到另种比较合理的比对 占：ta cca gt 一一 t ：一一cc cgtaa 1 2 4 多重序列比对( m u l t i p l ea li g n m e n t ) 在序列比对时，两两比对远远不能满足当今生物研究的需要，难以找出多 条序列的共性，就要求我们进行多重序列比对 多重序列比对：就是参加比对的序列数目不止两条，通过字母的匹配( m a t c h ) ， 删除( d e l e t e ) 和插入( i n s e r t ) 操作，通过比对找出多条序列的共性 与序列两两比对不同，多重序列比对的目标是找出多条序列的共性多重序 列比对，是生物信息学研究一个主要的方法，随着d n a 测序方法的快速发展， 将未知的序列同整个数据库的已知序列进行比对，从而找出相似序列，就成了 最常用的方法同时多序列比对还可以辅助检查一个序列家族中的全局相似性 和进化亲缘关系通过序列的多重比对，找出相对保守的子序列，就是家族的 特征序列当遇到一条未知序列时，通过比对，判断这条序列是不是属于这个 家族 例1 3有五条d n a 序列 工1 = a t t g c c a t t x 2 = a t g g c c a t t x 3 = a t c c a a t t t t ，= a t c t t c t t = a c t g a c c 经过多重序列比对我们可以得到一个比较好的比对，如下 工1 = at tgccatt 一一 工2 = atg gc catt 一一 x 3 = atc c aat t tt 一 一= atctt ctt 一一 工5 = ac tg a cc 一一一 字母的排列，构建个动态规划矩阵，由于实际数据利用多维的动态规划矩 阵进行序列比对相当困难，因此就需要利用好的算法来降低复杂度，多序列比对 生物序列比对第一章序列比对 算法基本上都是基于动态规划思想的相对于两两序列的算法，多序列比对算法 发展的还不够成熟比较常用的有隐马尔可夫模型，它是一种基于统计理论的多 序列比对模型 1 3 得分计划与得分值 在例1 2 中，两条序列经过比对，得到两个比对结果一般的在进行序列比 对时，不同的比对方法可能得到多个比对结果，甚至对同一组序列用同一种比对 方法，也可能会得到多个比对结果，那么怎样从这些比对结果中筛选出更加合理 的比对就非常重要，我们可以为每个比对结果进行打分，得分越高的比对，说明 序列的匹配效果越好，就越容易找出特征片段，进行序列比对的目的就转化为找 出一个得分最高的比对，并看作是最优比对 序列比对的得分方案很多，比较简单的得分方案可以按照下面的方法定义： f s ( a ，6 ) = o ，a ；e b s ( 口，n ) = l i s ( - ，n ) = s ( a ，一) = - i 在这个得分方案下，计算例1 2 两个比对结果的得分，找出得分高的比对 s ：tac ca gt 一一 比对得分等于3 ( ：_ 1 + 0 + l + 卜0 + l + 1 + 0 + 0 ) 占：t ac ca gt 一一 比对得分等于4 ( ：0 + 0 + 1 + 1 + 0 + 1 十1 十0 + 0 ) t ：一一c c c g t a a 因而，我们认为第二个比对优于第一个比对，并且把 s ：tac cagt 一一 t ：一一cc c g ta a 看作是序列j ：t a c c a g t ，t ：c c c g t a a 的最优比对 现实中的得分安排需要考虑到各种各样的因素，得分标准可能由生物化学特 性( 如带电性、疏水性) 、物理特性( 如分子质量、形状) 、动力学特性( 如翻 转率) 或次级结构( 如q 一螺旋，1 3 一链，转角，开放卷曲) 决定氨基酸分类 同样基于密码子的差异和相似的三级结构 4 生物序列比对第一章序列比对 下面以单序列得分为例介绍得分计划，两条序列的比对得分可化成单序列得 分例如两条d n a 序列比对，每条序列均含有四个字母类型( a ，c ，g ，t ) ，定义a 与a 的比对为字母a l ，a 与c 的比对为字母a 2 ，依次类推可以共定义4 x 4 = 1 6 个 字母类型，两条序列的比对得分即转化为单序列各字母得分的情况给定字母集 口l ，口2 ，口，) 以及对应字母得分 岛，s 2 , - - ，s ，) 在单序列得分的情况下，对于核苷 酸来说，r = 4 ( 核苷酸只有四种字母：a ，c ，g ，t 对于两条核苷酸序列比对r = 1 6 ) ； 对于嘌呤与嘧啶比较，r = 2 ：对于密码子，r = 6 1 ：对于标准氨基酸，r = 2 0 ：对于一个 氨基酸的化学分类 脂肪族的，芳香族的 ，r = 8 ：对于可能由于氨基酸引起的带 电性，r = 3 不同的得分方案，可以构建不同的得分矩阵，如：等价矩阵，疏水矩阵，p a m 矩阵，b l o s u m 矩阵等得分矩阵是序列进行比较的基础，有了得分矩阵，可以较 快的找出最优比对，同时，在不同的得分矩阵下，最优比对结果可能不同，这就 要求我们根据比对目的，合理构造出得分矩阵 1 4 序列比对的算法 序列比对是生物信息学研究的重点，也是难点，进行序列比对最直接的方法 就是将待比对的两条序列所有可能的比对，一一列举出来，再利用得分函数，计 算出每个比对的得分，从中选出得分最高的比对，也就是我们寻找的最优比对 假定有两条序列j ，t ，长度分别为m ，一，两条序列之间的比对数目是m ，以的 一个指数函数当序列长度比较短时，我们可以较快的列举出所有比对但是，当 序列长度比较长时，即，m ，n 比较大时，要想列举出所有可能的比对就非常困难， 就是在计算机上也是很难实现的，必将花费大量的时间和人力因而，就要求我 们利用有效的序列比对的算法，利用算法对生物学中的海量数据进行处理，较快 的找出序列比对的最优比对生物学发展至今，对序列比对算法的研究也逐步加 大，已经产生了很多序列比对的算法 两序列的比对算法有：基于动态规划思想的最优全局全局比n e e d l e m a n 生物序列比对第一章序列比对 w u n s c h 算法，最优局部比对的s m i t h w a t e r m a n 算法：基于启发式思想的f a s t a 算法和b l a s t 算法；m u m e r 算法；遗传算法等在现有的两两序列比对算法中， b l a s t 算法是当今发展的比较成熟的序列比对算法 多序列比对方法有：c l u s t a l 算法，星比对算法，树比对算法，基于隐马尔 可夫模型的多序列比对方法等其中，隐马尔可夫模型是一种智能化算法，隐马 尔可夫模型应用在生物信息学领域上，已经取得了很大成就，有力的推动生物信 息学的发展 6 生物序列比对第二章b l a s t 算法 第二章b l a s t 算法 。 b l a s t 算法，即基本的局部序列比对搜索工具，是由h l t s c h u l 等人在1 9 9 0 年提出的，主要研究对于一个固定长度的无间隙比对，从数据库中找出与研究序 列相似的子序列，同时，b l a s t 算法可以给出所有比对得分高于门槛参数的序列 片段b l a s t 算法的一个非常重要的特点就是建立在严格的统计学理论基础上， 我们将要从两个方面来证明这个理论 2 1b l a s t 算法的统计理论基础 针对单序列比对的情况，设字母q = d i 口2 ，口，) ，q ，吒，口，之间是相互独 立的特别地，对d n a 序列，q = a ，g ，c ，t ) 假定，每个字母q 的得分记为毛，i = l ，2 ，字母q 的在序列中出现的概率记 为研，i = l ，2 ，r 字母q 的背景概率为研，i = l ，2 ， 我们要找出具有最大得分的 序列片段，这些片段称为h i g h s c o r i n gs e g m e n tp a i r s ( 最大得分片段) ，这些 片段的得分称为最大片段得分片段的长度并不是预先给定的，它取决于所给的 数据 为了证明b l a s t 算法，我们首先给出两个合理的假设： 假设一，至少有一个得分为负，即：存在墨“，s 2 ，墨) ，使得墨 o ； r 假设二，得分的数学期望为负，即：e = b 墨 o i = i 对于上面的两个假设，我们作下面的解释 注1 对于假设i ，若任意岛 1 这时 层= p i s i f ；l rr = a 一口乃凶) f = li = l = e a e 0 因而，假设l 和假设2 都是比较合理的 2 2 利用极值分布证明b l a s t 算法 下面我们来证明，对于给定长度n 的序列，墨，i = l ，2 ，n 为序列中字母的得 分，五看作是独立同分布的随机变量，其最大得分片段得分膨。的分布为极值分 布考虑部分和过程：岛= o ，瓯= ：l 五， m n = s u p ( s _ f 一品) o y 对按如下方式划分，令 k o = 0 ，k v = m i n 七：七j o l + 1 ，叉一芦b 0 ) ，v = l ，2 ， g 2 景( s i 一) ，y = 1 ，2 ，“ 由最大得分片段的定义可以知道帆= m a x o , ，幺，鳊( 。) ，q 。) 其中，r ( n ) 表示从1 到玎的全部偏移次数，q 为从( 。) 到n 的最大残差 令f ( ) ，) = p ( g 力，a = e ( k i ) ，由波动理论知 彳= e x p & o ) ) 七专l n 生物序列比对第二章b l a s t 算法 因为e ( x ) 0 ，从而彳有限，即：彳 o o z + 为 晶) 首次取到正数所对应的值，盯+ 为此时既所对应的下标，因而z + = 毛， 令t ( 4 ) 2 p ( z 1 善) ，z 为 晶 首次驭负数对应的值，则z 一2 黾 引理2 1 假定x 满足前面提到的两个假设，则 当x 是不等距离散化时，c 牛= j ；! 骢 1 一f ( n s e x p ( 8 n 8 ) = 多：岳： 一v一 当x 是等距离散化时，c 宰2 j 翌 1 _ f ( 纠e x p ( 8 * y ) 5 手 v， 其中，0 e e 出】= 1 的正根， c：1-l(oo)(1-eexp(oz-) e z + e x p ( 伊z + ) ；矿 o o 】 4 e x p - 2 ( e 【e x p ( 矿& ) ；品 o 】+ p 最o ) ) ) ：，：1：!一 e x e x p ( 目x ) 定理2 1 设两条序列s 的长度分别为n 利用肘。表示两条序列比对的最大得分片 段的得分则 当x 是不等距离散化时 l i r ap m 一警5 工 = e x p ( - k e x p ( 一目x ) ) ( 2 1 ) h o 当x 是等距离散化时，且间距是万 e x p ( - k e x p ( - 矿工) ) ! 受1 n f p 坂一等工) 一， 黝l i r au p p 一警如h e x p ( 一k e x p ( 一0 工) ) 其中疋= e x p ( 一0 + 0 3 z 【 证明：考i 患m ( k m ) = m a x ( q l ，q 2 ，瓯) ，由于q l ，q ，级是独立的，得到 尸彤( ) 等州书( 警训】m 南引理2 1 9 生物序列比对 第二章b l a s t 算法 肼l i 嘲mp m ( k j 警)肼o。， = l i r ae x p m l o g 1 一f 粤+ 圳 m - - ) m = e x p ( 一c e x p ( 一矿功 由于彳= e ( k 1 ) ，令朋= n a 】表示，l 彳整数部分的整数部分，则 。l i r a 州( o 学州= 唧( 一c e x p ( n 又因为型兰坐骘1 ，且膨。是单调的r 因此 l i mp m ( k ) i n n n + z ) = e x p ( 一c e x p ( 一9 x ) 月 其中 万e x p 2 ( 研e x p ( 矿s k ) ；s k o i ；1j ；l 定理2 2 设两条序列s ，t 的长度分别为聊，力利用m ( s ，t ) 表示两条序列比对的最 大得分片段的得分当x 是不等距离散化时， 1 0 生物序列比对第二章b l a s t 算法 l i r ap m m - 雩x ：e x p ( 一k m ne x p ( 一矿工) h 其中，矿是乃p 鸭= l 的唯一正根 驴 证明方法与定理2 1 的证明方法类似，这里略去 例2 1 q = a ，c ，g ，t ) 为d n a 序列的字母集，各字母出现的概率均为，两条 序列比对得分计划为：相同字母比对得分为s l = l ，不同字母比对得分为s 2 = - 1 ， 则得分为1 懈f o p , = i i ，得分为一1 的概率为p 2 = i 1 ，利用公式可得r = l n 3 ， 由k 的计算公式可得k ：! 3 解：利用公式，可以计算最大片段得分帆超过某给定阈值的概率大于 j ：下1 1 17 + 工的概率为l e x p 一k e x p ( - a 功，对于不同的玎与x 可得相应的结果 五k刀js p 1 0 9 8 6 1 1 3 1 0 0 01 37 5 8 7 7 1 0 0 7 6 8 1 0 9 8 6 1 1 3 1 0 0 01 4 7 6 8 7 7 10 0 6 9 1 1 0 9 8 6 11 3 1 0 0 01 5 7 7 8 7 7 10 0 6 2 1 1 0 9 8 6 11 31 0 0 01 67 8 8 7 7 10 0 5 5 9 1 0 9 8 6 1i 31 0 0 01 77 9 8 7 7 10 0 5 0 2 l - 0 9 8 6 1 i 3 1 0 0 01 88 0 8 7 7 1 0 0 4 5 1 1 0 9 8 6 1 i 3 1 0 0 0 1 98 1 8 7 7 10 0 4 0 5 此时，在0 9 5 水平下，局部片段比对得分超过7 9 8 7 7 1 时应视为高得分 2 3 利用启发式证明k a r li n - a i t s c h u l 定理 设一列独立同分布的随机变量， y ，i z ) ，并且满足层( z ) o 令 x j = s u p i i + l - 。y j ，m = s u px j 。 1 s f s 一 为了求出m 。的分布，下面利用泊松分布族的启发式来证明，把 墨2 纠看 生物序列比对 第二章b l a s t 算法 作是来自独立同分布的强度为名6 的泊松过程n m l 引理2 2 设 m ：f r ) 是 0 ，t 上的泊松分布，则 p ( m r b ) e x p ( - 以t ) 引理2 3 设矿是满足e ( e x p 护z ) = l 的唯一正根，胁= p ( 五6 ) ，则 = f 机 ( 2 2 ) ( 2 3 ) 对于布朗运动，五= 一+ 仃蜀，我们来计算，这个布朗运动在正半轴上， 所花时间的分布情况 引理2 4 布朗运动置= t + a b t ，为其漂移速度，盯2 为方差，且为标准布朗运 动，ac r ，r ( 彳) 表示( 五：f 0 ) 在a 上所花时间的总和，则 姗扩1 岛r ( 。，回2 南， 。 声l o ”、77 i f i 在这里，我们只需用到在正半轴逗留时间的结论 令p o ，五：o ，r ( o ，o d ) ：r ，其中 “。 盯= 三，盯2 = - 三 e e x p ( - g f ) = 2 ( ( 1 + 2 b ) 2 + 1 ) 一1 后( x ) = 2 x2 ( x 2 ) 一2 0 ( x 2 ) iil 矽和是标准正态分布的密度和概率其中面( 功= l - o ( x ) 因此后( g ) = 日( 。= 吾盯= 丢 引理2 5 设五为平稳过程，6 是常数，叮2 ( 6 ) 为方差，s ( t ，6 ) 为在水平6 下觚总 的逗留时间，贝i js ( t ，6 ) 可以看成是参数为u 的泊松分布( g ) ，其中v = t 2 b ，即 s ( t ，6 ) 缈f s ( o ) ，d = f 九 ，既乃e g ，其中见= 以五6 ) 因而，对于偏移为- a ( b ) ，方差为盯2 ( 6 ) 的布朗运动。假定分布族为 生物序列比对第二章b l a s t 算法 脚 口c 6 = 鬻r ，因此我 门推导出 e c b = e 离胪嚣阶器 亿4 ， 定理2 3 长度为n 的生物序列，序列中字母的得分为岛，i = 1 ，2 ，忍，x 为随机变 量，取值为墨，i = 1 ，2 ，刀，令品= o ，最= ：l 墨，帆= s u p ( 昌一最) ，b 为常数， o 女列” 则当n 足够大时 p ( m 。6 ) e x p - k ne x p ( - e b ) ( 2 5 ) 其中k ：珲，：e x ，矿2 ：d x 仃。 证明：由引理知尸( 鸩 这就得到了最大得分片段得分的近似分布，但是推导过程与定理l 的推导过 程完全不同对于两两序列的情况，设两条序列s , t 长度分别记为m ，n s ：a l a 2 a m ，f ：岛如也，同样我们引入m 册表示最大得分片段的对应的得分，假定 比对不容许空格，因而，所得到的最优比对只能是序列上连续的一块区域可以 假定序列比对的最优比对为：a i a f + l 口f + 七，岛6 f + 1 也+ 七 定理2 4 设两条d h a 序列s ，t 长度分别记为m ，n ，变量x 取值为序列比对的得分 生物序列比对 第二章b l a s t 算法 值，则当聊，刀足够大时，帆。满足： p ( m 。刀b ) e x p 一k 朋ne x p ( - a b ) ( 2 5 ) 其中k ：珲，：e x ，盯2 ：d x 例2 2 设两条序列序列j = a l a 2 一口。，f = h a 屯，进行两两序列比对后，字母比对 得分为：等于l 的概率为p = 1 4 ，得分等于一1 的概率为r = _ i ，得分等于0 的概 率为g = 1 4 利用上面的结论，算得= e y = 一i i ，叮2 = d y = e y 2 一( 五y ) 2 = 嚣 k = 瓦b = 等= 嚣，喜见e x p c 锥，= ，代入数值算得口= h 暑= t n 2 ， p ( m 。b ) e x p 一k m ne x p ( - t b ) = e x p ( 一青m 珂e x p ( 一6 i n 2 ) 由此可以算得最大得分片段得分的近似分布 1 4 生物序列比对第三章基于隐马尔可夫模型( h m m ) 的多重序列比对 第三章基于隐马尔可夫模型( h m m ) 的多重序列比对 隐马尔可夫模型( h i d d e nm a r k o vm o d e l s ，h m m ) 是种基于概率的统计模型， 应用范围非常广泛，如：语言识别、生物学等本节主要考虑h m m 在生物序列比 对问题上的应用，利用h m m 来解决多重序列比对问题 3 1隐马尔可夫模型 离散有限的马尔可夫链( d i s c r e t e t i m ef i n i t em a r k o vc h a i n s ) 定义如下： ( 1 ) 有限集合s = 侈l ， ，其中s 。，是，为个所有可能的状态，对任 意时间点t = l ，2 ，3 ，马尔可夫链均对应状态空间s 中的一个状态，用s 表示 t 时刻的状态，其中墨爷i ，s 2 9 ooo ，曲 ( 2 ) 从时间f f + l ，状态从s 转移到邑对应概率毋= 尸( 墨专s ) ，如果毋仅 依赖于歹，与时间f 无关，那么所有既构成一个矩阵记为p ，成为状态转移 概率矩阵 p = p l l p t 2 p l n p 2 t p 2 2 p 2 n p n i p n 2 p n n ( 3 ) 将来状态仅与现在状态有关，而与过去状态没有关系，即 p = p ( 墨+ 1 = s si 葺= s ，五一1 = s k ，) = p ( 五+ 1 = s si x , = 墨) p ( x t n = si ，x t = s 0 = = _ _ _ _ - _ _ - _ - - 。_ _ _ _ ；。_ _ _ _ i _ _ _ - - _ 一 p ( 五+ l = 马) 五就形成一条马尔可夫链，可以看出，一条马尔可夫链仅与初始状态和转 移概率矩阵有关 设一条d n a 序列工= 西砭。，其中而，屯，而s ，假定对任意的字符口，b ，转 换概率定义如下：= p ( 五+ l = 6 i 墨= 口) 生物序列比对 第三章基于隐马尔可夫模型( h m m ) 的多重序列比对 我们把 薯) 看作一个随机过程，随机变量未来的取值，仅仅与当前取值有 关，而与过去取值无关，那么产生序列x = x l x 2 的概率记为p ( z ) 则 尸( x ) = p ( x o p ( x2 i 而) p ( z3 lx 2 ，而) p ( x 。i 一l ，x 1 ) = e ( x i ) p ( 工2 ix 1 ) p ( x 3i 恐) ，( 而i 毛一1 ) = p ( x i ) p x i x2 p 札m 假定一条原核生物的d n a 序列x = x l x 2 ，其中而，而，x , n s ，包含玎个基 因，我们很自然的按照如下形式定义转移概率p 曲为：序列中字母曲连续出现的 次数与字母a 出现次数的比值，即 # a b 毒 2 y a k 特别地定义 序列中a 的出现的次数 儿2 1 丽丽矿 为了处理问题的方便，需要在马尔可夫链上增加一个开始状态和结束状态 下面用字母b 表示开始状态，用字母e 表示结束状态这样，d n a 序列x 就可以 写成x = 砚而矗e ，产生序列x = 而恐。乇的概率就可以写成： 尸( x ) = p ( 而ib ) p ( x2 ix t ) p ( x3 1x 2 ) p ( x 。i 一i ) p ( ei ) =p b x ipj i j 2 p 工j i 工，p _ i ，层 图3 1 1 6 生物序列比对第三章基于隐马尔可夫模型( h m m ) 的多重序列比对 例3 1 有如下四条d n a 序列，利用隐马尔可夫模型估计相关参数 l p e a 2 j ， 4 p a a2 i ， 1 j l p c 2 万 5 p g a2 i ， ， x i ：c t a t t g a t 吃：a a a g a c t t c 而：c g c g a t g a a c c c x 4 ：a a g c a t g p e c = i 1 ，p e g = o ，p e t = o = 西2 ，= 西2 ， = 一5 ，p a o p a cp a gp a t 1 3 = o 2 西，2 西 2 一， 2 o 2 2 22 p c c2 万p c o 2 歹p e r2 歹p c o5 万 p g c = 号，p g o = o ，p 研= o ，p = o 1l 33l m 2 而趾2 而，胤2 而舶2 而，2 ； 因此，上面的计算结果表明，各个字母间不是全连通的，下面我们用图形来 表示出来 图3 2 一个隐马尔可夫模型( b m m ) 描述为：三元素集合a = ( a ，b ，万) ，其中的字母 作如下说明： ( 1 ) 有n 个状态，记为s = s 。，岛，) ，墨表示t 时刻的状态，其中 s t 俗l ，岛，) ( 2 ) 每个状态有m 个观测符号，矿= k ， 如：基因的编码序列称为外 1 7 生物序列比对第三章基于隐马尔可夫模型( h m m ) 的多重序列比对 显子，非编码序列称为内含子，那么基因的状态就可以用集合集合 s = s ，品 来表示，其中墨表示外显子，岛表示内含子 5 端 ，其中 b j ( ) = p ( 咋i 吼= s ) ，歹= l ，2 ，k = 1 ，2 ，m ( 5 ) 初始状态分布：t = 娩) ，_ ，= l ，2 ，n ，乃= p ( g l = 墨) 可以用一个有向无环图( d a g ) ，把隐马尔可夫模型简单的表示为： 隐马尔可夫链( h m m ) 观测值0 1 0 2 0 3 有向的无环( d i r e c t e d - a c i r c l e g r a p h ，d a g ) 图3 3 3 2 隐马尔可夫模型的三个基本问题 ( 1 ) 评价问题：给定观测序列d = q 0 2 呀和一个模型名= ( 彳，b ，万) ，怎样计算由 此产生的观测序列的概率p ( oia ) ? ( 2 ) 解码问题：给的观测序列0 = o 1 0 2 听和一个模型元= ( 彳，b ，万) ，如何确定一 条适合产生序列的状态序列q = 吼能4 r ，即找q 。= m a x p ( q i o , , , 1 ) ，这里的 生物序列比对第三章基于隐马尔可夫模型( h m m ) 的多重序列比对 求和是对所有q 求和，共有n r 利，情况 ( 3 ) 学习问题：对于给定的一个观测序列0 = d 。0 2 0 t ，如何确定模型 名= 一，b ，石) ，使得p ( o i 五) 最大? 三个基本问题相应的解法 ( 1 ) 评价问题的解法：向前或向后算法 给定一个隐马尔可夫模型a = ( 4 ，丑，石) 和给定观测序列0 = o 1 0 2 o r ，假定产生观 测序列的路径未知，怎样计算模型兄= ( a ，b ，刀) 产生的概率p ( oi 兄) 记q = q l q 2 。聊，o = o 1 0 2 o t ， 所以 p ( 0 1 2 ) = p ( o ，q i a ) q = p ( o i q ，2 ) p ( q i a ) q = p ( o 1 0 2 o r g i q 2 q r ，元) p ( g 眦q q r12 ) q p ( q ia ) = p ( q l q 2 q ri 见) = 尸( 口ll2 ) p ( q 2l 1 , ，q 1 ) 一，( 劬lq l q 2 口卜i ，2 ) p ( oiq ，2 ) = p ( 0 1 0 2 0 tq l q 2 q 7 - ，1 , ) = p ( o ll q ，2 ) p ( 0 2i d l ，q ，a ) p ( d 7 i o l 0 2 0 t - i ，q ，a ) = p ( o liq l2 ) p ( 0 2i 吼，z ) p ( 卟iq r , 1 , ) p ( o l a ) = p ( 吼12 ) p ( q 2la ，q i ) p ( q rlq i q 2 q r - i ，；o p ( o iiq i ，名) p ( d tq 2 ，a ) p ( o riq r ，a ) 口 = 芝：以吼1 2 ) p ( o - j q l ，五) p ( 仍i 五，q 1 ) p ( 0 2 q 2 。五) 。p ( q r i 细q 2 目l - i ，2 ) p ( o rj 牙r ，2 ) q ：艺羔p ( g 。12 ) p ( d lig i ，旯) j p ( 9 2i 无9 1 ) p (