生物信息学基础 第三章课件.ppt

第三章 序列比对,主讲教师：丁彦蕊单位：信息工程学院,序列比对：寻找两条或者两条以上的序列中各个字符的一一对应关系。序列比对的根本任务：发现序列之间的相似性寻找共同区域辨别序列之间的差异目的：相似序列 相似的结构，相似的功能 判别序列之间的同源性推测序列之间的进化关系,第一节 序列的相似性,同源（homology）- 具有共同的祖先(P81)直向同源（Orthologous ）:共同祖先的直接后代（没有发生基因复制事件）之间的同源基因称为直向同源。 共生同源（paralogous ）:两个物种A和B的同源基因，分别是共同祖先基因组中由复制事件而产生的不同拷贝的后代，这被称为共生同源基因。相似（similarity）同源序列一般是相似的 相似序列不一定是同源的 进化趋同（同功能）,序列比较的基本操作是比对（Alignment） 两个序列的比对是指这两个序列中各个字符的一种一一对应关系，或字符的对比排列 。,设有两个序列：GACGGATTAG，GATCGGAATAG,Alignment2: GA CGGATTAGGATCGGAATAG,Alignment1:GACGGATTAG GATCGGAATAG,1、字母表和序列,字母表4字符DNA字母表：A, C, G, T扩展的遗传学字母表或IUPAC编码单字母氨基酸编码,扩展的遗传学字母表或IUPAC编码,特定的符号 代表字母表 A* 代表由字母表A中字符所形成的一系列有限长度序列或字符串的集合 a、b、c代表单独的字符 s、t、u、v代表A*中的序列 |s|代表序列s的长度,为了说明序列s的子序列和s中单个字符，在s中各字符之间用数字标明分割边界例如，设s=ACCACGTA，则s可表示为 0A1C2C3A4C5G6T7A8 i:s:j 指明第i位或第j位之间的子序列, 当然，0 i j |s|。 子序列0:s: i 称为前缀，即prefix(s,i)子序列 i:s:|s|称为后缀，即suffix(s, |s|-i+1),i:s: i 为空序列j-1:s:j 表示s 中的第j 个字符，简记为sj子序列与子串（p82）子序列：选取s中的某些字符（或删除s中的某些字符）而形成s的子序列例如： TTT 是 ATATAT的子序列。,s的子串：是由s中相继的字符所组成。 例如：TAC是AGTACA的子串，但不是TTGAC的子串（是子序列）。 子串是子序列子序列不一定是子串,字符串操作,字符串连接操作：两个序列s和t的连接： s + + t例如：ACC+CTA = ACCCTA 字符串k操作 删除字符串两端的字符 其定义如下：prefix(s,l) = sk|s|-lsuffix(s,l) = k|s|-ls i:s:j = ki-1sk|s|-j,序列比较可以分为四种基本情况(P83)（1）两条长度相近的序列相似 找出序列的差别（2）判断一条序列的前缀与另一条序列的后缀相似（3）判断一条序列是否是另一条序列的子序列（4）判断两条序列中是否有非常相似的子序列,2、编辑距离（Edit Distance),GCATGACGAATCAG TATGACAAACAGC,GCATGACGAATCAG TATGAC-AAACAGC,说明两条序列的相似程度 定量计算,两条序列的相似程度的定量计算相似度，它是两个序列的函数，其值越大，表示两个序列越相似 两个序列之间的距离。距离越大，则两个序列的相似度就越小,字符编辑操作（Edit Operation）,字符编辑操作可将一个序列转化为一个新序列 Match（a，a）Delete（a，-） Replace（a，b）Insert（-，b）,直接距离计算的不足,扩展的编辑操作,ACCGACAATATGCATA ATAGGTATAACAGTCA,ACCGACAATATGCATA ACTGACAATATGGATA,第二条序列头尾颠倒,CTAGTCGAGGCAATCTGAACAGCTTCGTTAGT,？,反向互补序列,RNA发夹式二级结构,3、通过点矩阵进行序列比较“矩阵作图法” 或 “对角线作图”, 序列1 ,序列2 ,实 例, 序列1 ,序列1 ,自我比较,滑动窗口技术两条序列中有很多匹配的字符对，因而在点矩阵中会形成很多点标记。,滑动窗口技术 使用滑动窗口代替一次一个位点的比较是解决这个问题的有效方法。 假设窗口大小为10，相似度阈值为8，则每次比较取10个连续的字符，如相同的字符超过8个，则标记 基于滑动窗口的点矩阵方法可以明显地降低点阵图的噪声，并且明确无误的指示出了两条序列间具有显著相似性的区域。,（a）对人类（Homo sapiens）与黑猩猩（Pongo pygmaeus）的球蛋白基因序列进行比较的完整点阵图。（b）利用滑动窗口对以上的两种球蛋白基因序列进行比较的点阵图，其中窗口大小为10个核苷酸，相似度阈值为8。,(a) (b),具有连续相似区域的两条DNA序列的简单点阵图,4、 序列的两两比对,序列的两两比对(Pairwise Sequence Alignment）:通过字符匹配和替换，或者插入和删除字符，使两条序列达到一样的长度，并使两条序列中相同的字符尽可能一一对应。,s:AGCACACAAGCACACA t:ACACACTAACACACTA Match(A, A)Match(A, A)Delete(G, - )Replace(G, C)Match(C, C)Insert( -, A)Match(A, A)Match(C, C)Match(C, C)Match(A, A)Match(A, A)Match(C, C)Match(C, C)Replace(A, T)Insert( -, T)Delete(C, -)Match(A, A)Match(A, A)图3.6 序列AGCACACA和ACACACTA的两种比对结果,Alignment -1 Alignment -2,不同编辑操作的代价不同编辑操作定义函数w，它表示“代价（cost）”或“权重（weight）”。对字母表中的任意字符a、b，定义 w (a, a) = 0 w (a, b) = 1 a b w (a, -) = w ( -, b) = 1,也可以使用得分（score）函数来评价编辑操作 p (a, a) = 1 p (a, b) = 0 a b p (a, -) = w ( -, b) = -1,概念,两条序列s 和 t 的比对的得分（或代价）等于将s 转化为t 所用的所有编辑操作的得分（或代价）总和；s 和t 的最优比对是所有可能的比对中得分最高（或代价最小）的一个比对；s 和t 的真实距离应该是在得分函数p值（或代价函数w值）最优时的距离。,例如：s:AGCACACAt:ACACACTA cost=2 s:AGCACACA t:ACACACTA score (s，t）= 5序列比对的目的是寻找一个得分最大（或代价最小）的比对。,5、打分矩阵（Weight Matrices）(P87),（1）核酸打分矩阵设DNA序列所用的字母表为 = A，C，G，T a. 等价矩阵（相同核苷酸得分为1，不同核苷酸替换得分为0） b. BLAST矩阵（相同核苷酸得分为+5，不同核苷酸得分为-4）c. 转移矩阵（transition，transversion） （嘌呤：腺嘌呤A，鸟嘌呤G；嘧啶：胞嘧啶C，胸腺嘧啶T）,表3.1 等价矩阵表,表3.3 转移矩阵,表3.2 BLAST矩阵,（2）蛋白质打分矩阵,（i）等价矩阵,其中Rij代表打分矩阵元素i、j分别代表字母表第i和第j个字符。,（ii） 氨基酸突变代价矩阵GCM GCM矩阵通过计算一个氨基酸残基转变到另外一个氨基酸残基所需的密码子变化数目而得到，矩阵元素的值对应于代价。如果变化一个碱基就可以使一个氨基酸的密码子改变为另一个氨基酸的密码子，则这两个氨基酸的替换代价为1，如果需要两个碱基的改变，则替换代价为2，以此类推。GCM矩阵常用于进化距离的计算，其优点是计算结果可以直接用于绘制进化树，但在蛋白质序列比对尤其是相似程度较低的序列比对中很少使用。,（iii）疏水矩阵根据氨基酸残基替换前后疏水性的变化得到的矩阵。如果氨基酸A被氨基酸B替换后，疏水性变化不大则替换得分高，否则替换得分低。,（iv）PAM矩阵（Point Accepted Mutation）统计自然界中各种氨基酸残基的相互替换率。如果两种特定的氨基酸之间替换发生得比较频繁，则这一对氨基酸在得分矩阵中的互换得分就高。PAM矩阵基于进化原理，建立在进化的点接受突变模型基础上，通过统计相似序列中的各种氨基酸替换发生率而得到的矩阵。,PAM矩阵（Point Accepted Mutation） 基于进化的点突变模型 一个PAM就是一个进化的变异单位, 即1%的氨基酸改变,这类矩阵里列出同源蛋白质在进化过程中氨基酸变化的可能性。 这类矩阵式基于进化原理的 证据： 编码相同蛋白质的基因随着进化发生分歧，相似度降低。 科学 用得多,矩阵集合- PAM-N如，PAM120矩阵用于比较相距120个PAM单位的序列。一个PAM-N矩阵元素（i，j）的值： 反应两个相距N个PAM单位的序列中第i种氨基酸替换第j种氨基酸的频率。,针对不同的进化距离采用PAM 矩阵,序列相似度 = 40% 50% 60%,| | |打分矩阵 = PAM120 PAM80 PAM 60,PAM250 14% - 27%,（v） BLOSUM矩阵 （Blocks Amino Acid Substitution Matrices）通过统计相似蛋白质序列的替换率得到的。PAM矩阵是从蛋白质序列的全局比对结果推导出来的，而BLOSUM矩阵是从蛋白质序列块比对而推导出来的。,BLOSUM 62,第二节 两两比对算法,1、序列两两比对基本算法,直接方法 生成两个序列所有可能的比对，分别计算代价函数，然后挑选一个代价最小的比对作为最终结果。本质问题：优化动态规划寻优策略动态规划算法（Dynamic Programming）（P93）,最短路经问题,起点,终点,C1,C2,W1,W2,路径1：C1 + w1 ？路径2：C2 + w2 ？ 取最小值！,算法求解: 从起点到终点逐层计算,利用动态规划方法求解序列的两两比对,起点,终点,ATTCCGAAGA AGTCGAAGGT,ATTCCGAAG AGTCGAAGG,AT,+,（1）,ATTCCGAAGA AGTCGAAGG,-T,+,（2）,ATTCCGAAG AGTCGAAGGT,A-,+,（3）,求解过程,起点,终点,ATTCCGAAGA AGTCGAAGGT,从两个序列前端开始 逐步推进 直到两个序列的末端。,序列S: i-1 i i+1,序列t: j-1 j j+1,序列S: i-1 i i+1,序列t: j-1 j j+1,Case1:匹配（si，tj ),中间过程：比对0:S:i 与 0:T:j,序列S: i-1 i i+1,序列t: j-1 j j+1,序列S: i-1 i i+1,序列t: j-1 j j+1,Case2:删除（si， -）,序列S: i-1 i i+1,序列t: j-1 j j+1,序列S: i-1 i i+1,序列t: j-1 j j+1,Case3:插入（ -，tj ）,设序列s、t的长度分别为m和n。考虑两个前缀 0:s:i 0:t:j 假如已知序列0:s:i 和0:t:j 所有较短子列的最优比对，即已知：（1）0:s:(i-1) 和 0:t:(j-1) 的最优比对（2） 0:s:(i-1) 和 0:t:j 的最优比对（3） 0:s:i 和 0:t:(j-1) 的最优比对则0:s:i和 0:t:j 的最优比对一定是上述三种情况之一的扩展（1）替换（si，tj）或匹配（si，tj ) ，这取决于si 是否等于tj ；（2）删除（si， -）；（3）插入（ -，tj ）。,令：,为序列0:s:i和与序列 0:t:j 比对的得分按下述方法求解,其初值为：,for i=1 , 2 ,., m,for j=1 , 2 ,., n,距离矩阵按照上述方法，对于给定的得分函数p，两个序列所有前缀的得分定义了一个(m+1)(n+1)的距离矩阵D = ( d i , j )其中d i , j = S (0:s:i , 0:t:j ),d i , j的计算公式如下：,d i , j 最小值的三种选择决定了各矩阵元素之间的关系，用下图表示：,di,j,di,j-1,di-1,j,di-1,j-1,S (0:s:i , 0:t:j ),S (0:s:i-1 , 0:t:j ),S (0:s:i-1 , 0:t:j-1 ),S (0:s:i , 0:t:j-1 ),动态规划算法计算过程： 计算过程从d 0 , 0开始 可以是按行计算，每行从左到右，也可以是按列计算，每列从上到下。 当然，任何计算过程，只要满足在计算d i , j时 d i-1 , j、d i-1 , j-1、和d i, j-1都已经被计算这个条件即可。 在计算d i , j后，需要保存d i , j是从d i-1 , j、d i-1 , j-1、或d i, j-1中的哪一个推进的，或保存计算的路径，以便于后续处理。上述计算过程到d m , n结束。,最优路径求解：与计算过程相反 从d m , n开始，反向前推。 假设在反推时到达d i ,j，根据保存的计算路径判断d i , j究竟是根据d i-1 , j、d i-1 , j-1、和d i, j-1中的那一个计算而得到的。找到这个点以后，再从此点出发，一直到d 0 , 0为止。 走过的这条路径就是最优路径（即代价最小路径），其对应于两个序列的最优比对。,计算过程：（1）初始化,计算过程：（2）反复计算按列计算,计算过程：（2）反复计算按行计算其他方式,计算过程：（3）求最佳路径,例：s = AGCACACAt = ACACACTA,得分矩阵D （99）,初始化,计算d（2，2）,最终的得分矩阵及序列比对,AGCACACA| | |ACACACTA,举例,例 1: Si 和 Tj对齐,S: C A T T C A C T: C - T T C A G,i - 1,i,j,j -1,S: C A T T C A - C T: C - T T C A G -,例 2: Si 中加入空位,i - 1,i,j,S: C A T T C A C - T: C - T T C A - G,例 3: Tj 是入空位,i,j,j -1,计算过程,C(n,m),C(0,0),C(i,j), C T C G C A G C,A,C,T,T,C,A,C,+10 匹配, -2 不匹配, -5 空位,0 -5 -10 -15 -20 -25 -30 -35 -40,10,5, C T C G C A G C,A,C,T,T,C,A,C,回溯得到最佳的比对,C - T C G C A G CC A T - T C A - CC - T C G C A G CC A T T - C A - C,第一种比对方式,第二种比对方式,2、子序列与完整序列的比对,-AGCT-ATGCAGCTGCTT,目标：使S(s, i:t:j ) 最大,序列S:,序列t:,i j,不计前缀0:t:i 的得分, 也不计删除后缀的j+1:t:|t|得分,局部序列比对给定两条序列0:s:m和0:t:n，从t中寻找一个子序列i:t:j使得S(s,i:t:j)最大.,不计前缀0:t:i 的得分处理第一行,不计删除后缀的j+1:t:|t|得分 处理最后一行,dm,j,dm,j-1,dm-1,j,dm-1,j-1,S (0:s:i , 0:t:j ),S (0:s:i-1 , 0:t:j ),S (0:s:i-1 , 0:t:j-1 ),S (0:s:i , 0:t:j-1 ),不计代价,距离矩阵初始化时，对第一行进行如下处理：d0,j = 0 for 0 j n 最后一行的计算应该是：,同样，d m, n依然是最优局部比对的得分，而匹配的子列i:t:j 按如下方式寻找： （1） j = min k d m ,k = d m ,n （2）由（m,j）反推比对路径，最终通过斜线（非空位）到达（0，i）。,（3-10）,（3-11）,3、寻找最大的相似子序列,目标：使dw (i:s :j, i:t:j ) 最大,序列S:,序列t:,i j,i j,数据结构：（m+1）（n+1）的矩阵 D但是，对数组元素含义解释与基本算法有所不同每个元素的值代表序列0:s:i 某个后缀和序列0:t:j 某个后缀的最佳比对。,这种局部比对不计前缀的得分，所以新的边界条件是：d0,j = 0 for 0 j n （3-12）di,0 = 0for 1 i m另外，由于0:s:i 和0:t:j 总有一个得分为“0”的空后缀比对，因此矩阵D中的所有元素大于或等于“0”，于是，新的递归计算公式为：,(3-13),寻找最佳比对的子序列在矩阵中找最大值该值就是最优的局部比对得分该值对应的点为序列局部比对的末点然后反向推演前面的最优路径，直到局部比对的起点。,TATA|TATA,4、准全局比较,所谓准全局比较就是在评价序列比对时不计终端“空缺”（end space，或空位）的得分或代价,序列1 长度为8序列2 长度为18（a）6个匹配，1个失配，1个空位 （b）8个匹配,情况1 ：不记s后面的空位与 t 后缀比对的得分在矩阵di,j中取最后一行的最大值.,序列S:,序列t:,i j,i j,空位后缀,情况2 ：不记s前面的空位与 t 前缀比对的得分将矩阵di,j中的第一行各元素值置为“0”,序列S:,序列t:,i j,i j,空位前缀,情况3：情况4：半全局比较算法与基本算法在计算di,j时的区别归纳为下列四个方面：（1）第一行初始值为“0”，表示不计第一个序列的前端空位；（2）寻找最后一行的最大值，表示不计第一个序列的末端空位；（3）第一列初始值为“0”，表示不计第二个序列的前端空位；（4）寻找最后一列的最大值，表示不计第二个序列的末端空位。,对于最后一行和最后一列的另一种处理办法是：最后一行的横向移动不被空位罚分最后一列的纵向移动也不被罚分这样，就可以允许在两条序列终端自由存在空位。 当矩阵D所有元素计算完以后，其右下角得值即为两条序列最终准全局比对的得分。,ACACTGATCG|ACACTG,5、关于连续空位的问题,K 阶空位 K个连续的空位字符 “-”ATG-A-T-C-A-GATG-ATCAGATGCAGTGCAATGATGTTTTTATCAG生物学意义 “插入” 或“删除” 突变突变次数连续空位可能对应于一次突变非连续空位对应于 多次突变,对于连续空位的代价是一个线性的函数。设p（k）代表空位得分函数，其中k是连续空位的个数，则： p（k）= -bk 这里b（0）是单个“空位”得分的绝对值。处理方法：任何一个比对可以被唯一地分为若干个相继的块。有三类块：（1）两个字符的比对（2）与序列s空位进行比对的t的最大连续字符序列（3） 与序列t空位进行比对的s的最大连续字符序列,上述算法的时间复杂度为O(n3)。比起标准算法，其多花的时间主要用于处理连续的空位。那么，是否可以改进连续空位的得分函数，而使得算法的时间复杂度降低为O(n2)呢？如果认为k个连续空位比k个孤立空位出现的可能性更大，则p(k) kp(1) （3-22）或更一般地，p(k1 + k2 + + kn ) p(k1) + p(k2) + +p(kn) (3-23)可以用下式重新计算连续“空位”的得分：p(0)=0（3-24） p(k) = h g(k-1), k1 （3-25） h0，g0，hg。,6、比较相似序列,相似序列快速比较算法 例如，有两个序列： s=GCGCATGGATTGAGCGA t=TGCGCCATGGATGAGCA 最优比对所对应的路径偏离主对角线，经过一段以后重新返回主对角线。,经验法则（针对蛋白质序列）： 如果两个序列的长度都大于100，在适当地加入空位之后，它们配对的相同率达到25%以上，则两个序列相关； 如果配对的相同率小于15%，则不管两个序列的长度如何，它们都不可能相关； 如果两个序列的相同率在15%25%之间，它们可能是相关的。,第三节 序列多重比对,多序列比对的定义设有k个序列s1, s2, . ,sk，每个序列由同一个字母表中的字符组成，k大于2。通过插入操作，使得各序列达到一样的长度。,多序列比对的目的发现多个序列的共性发现与结构和功能相关的保守序列片段,多序列比对的应用研究蛋白质之间的关系研究一个家族中的相关蛋白质研究相关蛋白质的保守区域研究蛋白质的结构和功能进行蛋白质结构预测推测各个序列的进化历史,1、SP（Sum-of-Pairs）模型,评价多重序列比对的结果,按照每个对比的列进行打分，然后加和(P106)处理每一列： k个变量的打分函数 用一个k维数组来表示该显式函数（类似于打分矩阵）期望：函数在形式上应该简单具有统一的形式不随序列的个数而发生形式变化,其中，c1,c2,ck是一列中的k个字符，p是关于一对字符相似性的打分函数。,逐对加和SP（sum-of-pairs）函数(P106),逐对计算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,另一种计算方式：先处理每一个序列对在处理序列对时，逐个计算字符对，最后加和则SP得分模型的计算公式如下：, 是一个多重比对ij是由推演出来的序列s i 和s j的两两比对,2、多重比对的动态规划算法,多重序列比对的最终目标是得到一个得分最高（或代价最小）的序列对比排列，从而分析各序列之间的相似性和差异。,前趋节点的个数等于2k - 1,假设以k维数组A存放超晶格，则计算过程如下： a 0, 0, ,0 = 0 a i = max a i - b + SP-score(Column(s, i, b),(3-37),(3-38),if bj = 1,if bj = 0,多序列比对的过程实际是一个递推过程，在计算每个晶格节点得分的时候，将其各前趋节点的值分别加上从前趋节点到当前点的SP得分，然后取最大值作为当前节点的值。,计算量问题对于k条序列的比对，动态规划算法需要处理k维空间里的每个节点，计算量与晶格中的节点数成正比，而节点数等于各序列长度的乘积，况且计算的每个节点依赖于前趋节点的个数，因此用动态规划进行多重序列比对的时间复杂度为O(2ki=1,.,k si)即O(2kNk),利用SP模型寻找最优多重序列比对是一个NP-完全类问题。NP-完全问题通常被认为是一些人们难以在有限的时间、空间内对问题求出最佳解得问题，几乎所有专家都认为不可能在多项式时间内准确求解的问题。,解决办法只求解规模比较小的问题利用动态规划、分支约束等技术减小搜索空间，提高求解问题的效率。针对具体问题的特点，根据实际情况，设计实用求解算法。采用近似算法或者启发式方法，如局部搜索、模拟退火、遗传算法等。,3、 优化计算方法（p110),标准动态规划算法存在的问题： 搜索空间大对于两两序列比对，最优的路径常常处于对角线附近，而对于多序列比对，最优的路径不可能在某个平面上，而是在某一个区域范围内。,利用人工智能空间搜索策略的剪枝技术，根据问题本身的特殊性将搜索空间限定在一个较小的区域范围内。若问题是搜索一条得分最高（或代价最小）的路径，则在搜索时如果当前路径的得分低于某个下限（或累积代价已经超过某个上限），则对当前路径进行剪枝，即不再搜索当前路径的后续空间。,4、星形比对(P112),星形比对的基本思想是：在给定的若干序列中，选择一个核心序列，通过该序列与其它序列的两两比对形成所有序列的多重比对，从而使得在核心序列和任何一个其它序列方向的投影是最优的两两比对。利用标准的动态规划方法求出所有si和sc的最优两两比对时间为O（2knk）将这些两两比对聚集起来并采用“只要是空位 ， 则永远是空位”的原则。,scs1s2sk,(sc, s1) (sc, s2) (sc, sk),两两比对 多重比对,如何选择核心序列？尝试将每一个序列分别作为核心序列，进行星形多重序列比对，取比对结果最好的一个。另一种方法是计算所有的两两比对，取下式值最大的一个： sim( si, sc ),例如，有5个序列： s1 = ATTGCCATT s2 = ATGGCCATT s3 = ATCCAATTTT s4 = ATCTTCTT s5 = ACTGACC,sc=s1,ATTGCCATT ATTGCCATT- ATTGCCATT ATTGCCATTATGGCCATT ATC-CAATTTT ATCTTC-TT ACTGACC-,ATTGCCATT-ATGGCCATT- ATC-CAATTTT ATCTTC-TT- ACTGACC-,星形比对是一种近似的方法，可以证明，用该方法所得到多重序列比对的代价不会大于最优多重序列比对代价的两倍,5、树形比对,k个待比对的序列 具有k个叶节点的树每个叶节点对应一个序列将序列赋予树的内部节点，可以计算树中每个分支的权值。权值代表对应分支连接的两个序列之间的相似性。所有权值的和就是这棵树寻找一种树的内部节点序列赋予方式，使得树的得分最大。,将CT、CG、CT分别赋予节点x、y、z，则树的得分为8。这里假设如果a=b，则p(a,b)=1， 否则p(a,b)=0，p(a,-)=-1。,CT,CG,CT,多重序列比对 两两序列比对 合并两个比对（比对的比对）,Alignment of alignments, AA算法 假设：有两个多重序列比对1、2，1代表序列s1、s2、si的多重比对，2代表序列t1、t2、tj的多重比对，（s1，s2，si）（t1，t2，tj）=代表s1和t1的两两比对，则计算与相一致的1和2比对的算法如下 ：（1）标定1的各列，如果s1在比对中对应位置的编辑操作不是插入或删除，则这些列分别标记为s1对应位置上的字符a1、a2、als1（ls1为序列s1的长度）；（2）标定2的各列，如果t1在比对中对应的位置编辑操作不是插入或删除，则这些列分别标记为t1对应位置上的字符b1、b2、blt1（lt1为序列t1的长度）；（3）对a1、a2、als1和b1、b2、blt1进行比对；（4）在所得到的比对中，对于1、2和中原来有插入或删除操作的位置，恢复其原有的实际字符或空位字符“-”。,例：1: s1 -H-LVV 2: t1 L-HCLV ：s1 -H-LVV s2 G-VLVC t2 VLHCL- t1 LHCLV- s3 GN-LVVAA算法的输出为-H-LVV-G-VLVC-GN-LVVL-HC-LV-V-HC-L分别对第1、2列和4、5列进行压缩，则最后结果为,HLVVGVLVGGNLVVLHCLV-VHCL-,对于n个序列的树形比对的基本算法过程如下：（1）初始化，对于每个序列，生成一个叶节点（2）利用AA算法合并两个节点，形成一个新节 点，合并的结果放在新节点中，原来的两 个节点作为新节点的子节点（3）反复执行（2），直到形成n个叶节点的树 根为止，根节点中的序列即为最终的多重 比对结果。,s1 s2 s3 s4,1,2,6、其它多重序列比对算法,一般渐进式比对方法所采用的过程：（1）先将多个序列进行两两比对，基于这些比较，计算得到一个距离矩阵，该矩阵反映每对序列的关系；（2） 利用距离矩阵，建立一棵“相关树”；（3）从最接近的一对序列出发，逐步归并形成比对的聚类，直到所有序列处理完。,例：,(LYCES, SPIOL 84), (YEAST, (XENLA,(RAT, MOUSE 96), HUMAN 83), CHICK71) 66), DROVI 58),相关树,多序列比对,目前使用最广泛的多重序列比对程序是ClustalW ClustalW是一种渐进的比对方法，先将多个序列进行两两比对，基于这些比较，计算得到一个距离矩阵，该矩阵反映了每对序列的关系,EBI的CLUSTALW网址是： http:/www.ebi.ac.uk/clustalw/,7、统计特征分析,对于所得到的多重序列比对，我们往往需要进行归纳分析，总结这些序列的特征，或者给出这些序列共性的表示,HLVVGVLVGGNLVVLHCLV-VHCL-,（1）保守序列表示序列每个位置上最可能出现的字符（或者所有可能出现的字符）ATNTSC (N - A,T,C,G ; S - G,C),（2）特征统计图（Profile） 令P=(P1,P2,PL)，P表示在的每一列上各种字符出现的概率分布Pj=(pj0,pj1,pj|A|)A代表字母表，Pjk代表字母表A中第k个字符在第 j 列出现的概率。 第0个字符是特殊的空位符号“-”。,ATTATAACTTCTTATACTTTAGAAT,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 （碱基）,利用保守序列或者特征统计图可以判断一个序列是否满足一定的特征给定一个序列s=a1a2am，定义字符a在第j位的代价为 其中，|A|代表字母表A的长度，Ak代表A的第k个字符，特别地A0代表空缺字符“-”。整个序列s的代价为,一条序列与特征统计图相对照，如果代价值小，说明该序列具有相应的特征，否则该序列不具备相应的特征。,第四节 DNA片段组装,大规模基因组测序得到待测序列的一系列序列片段这些序列片段覆盖待测序列序列片段之间也存在着相互覆盖或者重叠。,目标序列序列碎片,1、片段组装问题定义：给定一组取自特定字母表的字符串集合F，寻找一个最短的字符串s，使得F中的每一个字符串都是s的一个连续子串。这里，集合F的字符串相当于待组装的序列片段，而s则是序列片段组装的结果。,Input Answer ACCGT -ACCGT- CGTGC -CGTGC TTAC TTAC- TACCGT -TACCGT- TTACCGTGC,（1）碱基标识错误,4个主要问题,（2）不知道片段的方向,（3）存在重复区域,. .,（4）缺少覆盖,2、序列片段组装模型,序列片段组装过程：三个步骤（1）首先进行序列片段的两两比较，确定可能的片段之间的覆盖（或者重叠）；（2）确定所有片段统一的覆盖模式，即确定各个序列片段的相对位置；（3）最后确定片段组装结果，即确定目标序列。,（1）最短公共超串模型,三种片段组装模型,给定一个字符串集合F，求出一个最短的字符串S，使得对于所有属于F 的字符串f，S是 f 的超串（或者 f 是 S 的子串）。设F=ACT，CTA，AGT 则S=ACTAGT 是 F 的超串由于S必须是各片段严格的超串，因此不允许片段的实验误差，各片段的方向必须是已知的。,（2）重建模型,考虑到片段的误差和未知方向的问题近似子串假设f、g是代表两条序列的字符串，f 作为 g 近似子串的代价为：,S(g)代表 g 所有子串的集合，d为一般编辑距离。,设 f=GCGATAG, g =CAGTCGCTGATCGTACG, 则最佳的子序列比对如下-GC-GATAG-CAGTCGCTGATCGTACG,设 是一个介于0和1之间的数，称串f 是在误差下S 的近似子串，如果ds(f, S) f重建模型：给定一个字符串集合F，求一个最短的字符串S，使得对于所有属于F的字符串f，下式成立：min (ds(f, S), ds(f, S) ) f其中 f 是 f 的反向互补串。,（3）多重连续区模型,称一个多重序列比对是 t-contig，如果其最弱连接的交叠长度至少为 t。如果能够根据序列片段集合F构造一个t-contig，称F允许一个t-contig。多重连续区模型：给定一个片段集合F和一个整数 t（0），将F分割为最小数目的子集Ci，1ik，每个Ci允许一个t-contig。,目标序列序列碎片,不连续区域,设 F = GTAC，TAATG，TGTAA,