药物生物信息学.ppt

2020/5/29,1,药物生物信息学,序列比对（1）,2020/5/29,2,练习：,Alignment:GACGGATTAGGATCGGAATA,设有两个序列：GACGGATTAG，GATCGGAATA,GA-CGGATTAG,GATCGGAATA-,2020/5/29,3,序列比对的任务与目的,序列比对的根本任务：发现序列之间的相似性辨别序列之间的差异目的：相似序列相似的结构，相似的功能判别序列之间的同源性推测序列之间的进化关系,2020/5/29,4,本章内容,序列的相似性二重序列比对多重序列比对DNA片断组装,2020/5/29,5,第一节序列的相似性,字母表和序列编辑距离矩阵作图法序列的两两比对打分矩阵,2020/5/29,6,同源与相似,序列同源（homology）-具有共同的祖先直向同源（Orthologous）（不同种属）共生同源（paralogous）（同一种属）序列相似（similarity）同源序列一般是相似的相似序列不一定是同源的进化趋同（同功能）,2020/5/29,7,直向同源（a1inspeciesI,a1inspeciesII)共生同源（a1anda2inspeciesI）,进化趋同,水平转移,基因复制,2020/5/29,8,序列的相似性描述定性的描述序列的比对是一种关于序列相似性的定性描述，它反映在什么部位两条序列相似，在什么部位两条序列存在差别。最优比对揭示两条序列的最大相似程度，指出序列之间的根本差异。定量的描述相似度距离,GA-CGGATTAG,GATCGGAATA-,2020/5/29,9,序列比对的基本操作是比对（Alignment）两个序列的比对是指这两个序列中各个字符的一种一一对应关系，或字符的对比排列。,设有两个序列：GACGGATTAG，GATCGGAATAG,Alignment2:GACGGATTAGGATCGGAATAG,Alignment1:GACGGATTAGGATCGGAATAG,2020/5/29,10,1、字母表和序列,将生物分子序列抽象为字符串，其中的字符取自特定的字母表。字母表4字符DNA字母表：A,C,G,T扩展的遗传学字母表或IUPAC编码单字母氨基酸编码上述字母表形成的子集,2020/5/29,11,扩展的遗传学字母表或IUPAC编码,2020/5/29,12,20种标准氨基酸的英文简写,2020/5/29,13,代表字母表A*代表由字母表A中字符所形成的一系列有限长度序列或字符串或序列的集合a、b、c代表单独的字符s、t、u、v代表A*中的序列|s|代表序列s的长度,关于部分特定的符号的规定,2020/5/29,14,为了说明序列s的子序列和s中的单个字符，在s中各字符之间用数字标明分割边界例如，设s=ACCACGTA，则s可表示为0A1C2C3A4C5G6T7A8i:s:j指明第i位或第j位之间的子序列,0ij|s|i:s:i为空序列i-1:s:i表示s中的第i个字符，简记为si子序列0:s:i称为前缀，即prefix(s,i)子序列i:s:|s|称为后缀，即suffix(s,|s|-i+1),2020/5/29,15,子序列与子串s的子序列：选取s中的某些字符（或删除s中的某些字符）而形成s的子序列例如：TTT是ATATAT的子序列。s的子串：是由s中相继的字符所组成。例如：TAC是AGTACA的子串，但不是TTGAC的子串（是子序列）,子串是子序列（连续子序列）子序列不一定是子串,2020/5/29,16,字符串操作,字符串连接操作：两个序列s和t的连接：s+t例如：ACC+CTA=ACCCTA字符串k操作删除字符串两端的字符其定义如下：prefix(s,l)=sk|s|-lsuffix(s,l)=k|s|-lsi:s:j=kisk|s|-j,2020/5/29,17,序列比对的四种基本应用：（1）两条长度相近的序列相似找出序列的差别两个实验室测定同一序列，实验结果比较（2）判断一条序列的前缀与另一条序列的后缀相似大规模DNA测序中序列片断的组装（3）判断一条序列是否是另一条序列的子序列搜索特定模式（4）判断两条序列中是否有非常相似的子序列分析保守序列,2020/5/29,18,2、编辑距离（EditDistance),GCATGACGAATCAGTATGACAAACAGC,GCATGACGAATCAGTATGAC-AAACAGC,说明两条序列的相似程度定量计算,2020/5/29,19,两条序列的相似程度的定量计算相似度，它是两个序列的函数，其值越大，表示两个序列越相似两个序列之间的距离。距离越大，则两个序列的相似度就越小,2020/5/29,20,海明距离：两条长度相等的序列，Hamming距离等于对应位置字符不同的个数。,不足：（1）序列长度不同（2）未必反映两条序列的真正对应关系（3）DNA复制过程中的碱基插入、删除,2020/5/29,21,字符编辑操作（EditOperation）,字符编辑操作可将一个序列转化为一个新序列Match（a，a）字符匹配Delete（a，-）1序列删除1个字符a，或2插入空位Replace（a，b）2序列字符b替换1序列字符aInsert（-，b）1序列插入空位，或2删除1个字符b,2020/5/29,22,扩展的编辑操作,ACCGACAATATGCATAATAGGTATAACAGTCA,ACCGACAATATGCATAACTGACAATATGGATA,第二条序列头尾颠倒,CTAGTCGAGGCAATCTGAACAGCTTCGTTAGT,？,CTAGTCGAGGCAATCTCTTGTCGAAGCAATCA,替换为互补碱基,2020/5/29,23,反向互补序列,RNA发夹式二级结构,2020/5/29,24,3、通过点矩阵进行序列比对“矩阵作图法”或“对角线作图”,2020/5/29,25,2020/5/29,26,2020/5/29,27,序列1,序列2,实例,2020/5/29,28,序列1,序列1,自我比对,2020/5/29,29,滑动窗口技术两条序列中有很多匹配的字符对，因而在点矩阵中会形成很多点标记。,2020/5/29,30,滑动窗口技术使用滑动窗口代替一次一个位点的比对是解决这个问题的有效方法。假设窗口大小为10，相似度阈值为8，则每次比对取10个连续的字符，如相同的字符超过8个，则标记基于滑动窗口的点矩阵方法可以明显地降低点阵图的噪声，并且明确无误的指示出了两条序列间具有显著相似性的区域。,2020/5/29,31,（a）对人类（Homosapiens）与黑猩猩（Pongopygmaeus）的球蛋白基因序列进行比对的完整点阵图。（b）利用滑动窗口对以上的两种球蛋白基因序列进行比对的点阵图，其中窗口大小为10个核苷酸，相似度阈值为8。,(a)(b),2020/5/29,32,作业（一）1、直接比较下面两个序列的相似度；并观察是否有更好的比对。2、尝试用矩阵标记图做出下面序列的最优比对,1:ATTCGAGCCT2:CGTTCAGCTA,2020/5/29,33,4、序列的两两比对,序列的两两比对（PairwiseSequenceAlignment）按字符位置重组两个序列，使得两个序列达到一样的长度匹配和替换，插入和删除,2020/5/29,34,s:AGCACACAAGCACACAt:ACACACTAACACACTAMatch(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)序列AGCACACA和ACACACTA的两种比对结果,Alignment-1Alignment-2,2020/5/29,35,不同编辑操作的代价不同为编辑操作定义函数w，它表示“代价（cost）”或“权重（weight）”。对字母表中的任意字符a、b，定义w(a,a)=0w(a,b)=1abw(a,-)=w(-,b)=1,2020/5/29,36,使用得分（score）函数来评价编辑操作p(a,a)=1p(a,b)=0abp(a,-)=p(-,b)=-1,2020/5/29,37,概念：两条序列s和t的比对的得分（或代价）等于将s转化为t所用的所有编辑操作的得分（或代价）总和；s和t的最优比对是所有可能的比对中得分最高（或代价最小）的一个比对；s和t的真实距离应该是在得分函数p值（或代价函数w值）最优时的距离。,2020/5/29,38,例：,序列比对的目的是寻找一个得分最大（或代价最小）的比对,s:AGCACACAt:ACACACTAscore(s，t）=5,s:AGCACACAt:ACACACTAcost=2,2020/5/29,39,5、打分矩阵（WeightMatrices）,不同类型的字符替换，其代价或得分是不一样的，特别是对于蛋白质序列。某些氨基酸可以很容易地相互取代，而不改变理化性质。缬氨酸丙氨酸（得分高，代价小）赖氨酸丙氨酸（得分低，代价高）比较保守的替换比起较随机的替换更可能维持蛋白质的功能，且更不容易被淘汰不同打分方法对应不同的打分矩阵,2020/5/29,40,(1)核酸打分矩阵,核酸打分矩阵设DNA序列所用的字母表为=A，C，G，Ta.等价矩阵b.BLAST矩阵（最流行的核酸序列比对程序）,等价矩阵表,BLAST矩阵,2020/5/29,41,腺嘌呤A,尿嘧啶U,胞嘧啶C,鸟嘌呤G,胸腺嘧啶T,c.转移矩阵（转换transition，颠换transversion）（嘌呤：腺嘌呤A，鸟嘌呤G；嘧啶：胞嘧啶C，胸腺嘧啶T）,转移矩阵,2020/5/29,42,（2）蛋白质得分矩阵,（i）等价矩阵（ii）氨基酸突变代价矩阵GCM（iii）疏水矩阵（iv）PAM矩阵（PointAcceptedMutation）（v）BLOSUM矩阵（BlocksAminoAcidSubstitutionMatrices）,其中Rij代表打分矩阵元素i、j分别代表字母表第i和第j个字符。,2020/5/29,43,氨基酸突变代价矩阵（遗传密码GCM）,GCM通过计算一个氨基酸残基转变到另一个氨基酸残基所需的密码子变化数目而得到，矩阵元素的值对应于代价。GCM常用于进化距离的计算，其优点是计算结果可以直接用于绘制进化树。但在相似程度很低的蛋白质序列比对中很少使用。,2020/5/29,44,2020/5/29,45,疏水矩阵,该矩阵是根据氨基酸残基替换前后疏水性的变化而得到得分矩阵。若一次氨基酸替换疏水特性不发生太大的变化，则这种替换得分高，否则替换得分低。,2020/5/29,46,2020/5/29,47,PAM矩阵（PointAcceptedMutation）,基于进化的点突变接受模型通过统计相似序列比对中的各种氨基酸替换发生率而得到该矩阵。1978年美国Dayhoff研究小组研究了71个相关蛋白质家族的1572个突变，发现蛋白质家族中氨基酸的替换并不是随机的。这意味着，在进化历程上，相关的蛋白质在某些位置上可以出现不同的氨基酸，这些点突变已经被进化所接受。,2020/5/29,48,矩阵集合-PAM-N一个PAM就是一个进化的变异单位,即1%的氨基酸改变如，PAM120矩阵用于比对相距120个PAM单位的序列。一个PAM-N矩阵元素（i，j）的值：反应两个相距N个PAM单位的序列中第i种氨基酸替换第j种氨基酸的频率。,2020/5/29,49,针对不同的进化距离采用PAM矩阵,序列相似度=40%50%60%,|打分矩阵=PAM120PAM80PAM60,PAM25014%-27%,2020/5/29,50,归一化打分实例,2020/5/29,51,BLOSUM矩阵,BLOcksSUbstitutionisanotheraminoacidsubstitutionmatrix,firstcalculatedbyHenikoffandHenikoff(1992,PNAS)Foritscalculationonlyblocksofaminoacidsequenceswithsmallchangebetweenthemareconsidered.Theseblocksarecalledconservedblocks,2020/5/29,52,BLOSUM相似性度量,BLOSUM与PAM矩阵阶数含义相反，低阶BLOSUM更多用来比较亲缘较远序列。序列相似度=40%62%80%|打分矩阵=BLOSUM40BLOSUM62BLOSUM60,2020/5/29,53,BLOSUM62,2020/5/29,54,第二节二重序列比对基本算法,直接方法生成两个序列所有可能的比对，分别计算代价函数，然后挑选一个代价最小的比对作为最终结果。本质问题：优化动态规划寻优策略动态规划算法（DynamicProgramming）,1、序列两两比对基本算法,2020/5/29,55,最短路径问题,起点,终点,C1,C2,W1,W2,路径1：C1+w1？路径2：C2+w2？取最小值！,算法求解:从起点到终点逐层计算,2020/5/29,56,2020/5/29,57,求解过程,起点,终点,ATTCCGAAGAAGTCGAAGGT,从两个序列前端开始逐步推进直到两个序列的末端。,2020/5/29,58,序列S:i-1ii+1,序列t:j-1jj+1,序列S:i-1ii+1,序列t:j-1jj+1,Case1:匹配（si，tj),中间过程：比对0:S:i与0:T:j,2020/5/29,59,序列S:i-1ii+1,序列t:j-1jj+1,序列S:i-1ii+1,序列t:j-1jj+1,Case2:删除（si，-）,比对0:S:i与0:T:j,2020/5/29,60,序列S:i-1ii+1,序列t:j-1jj+1,序列S:i-1ii+1,序列t:j-1jj+1,Case3:插入（-，tj）,比对0:S:i与0:T:j,2020/5/29,61,设序列s、t的长度分别为m和n。考虑两个前缀0:s:i0: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）。,2020/5/29,62,令：,为序列0:s:i和与序列0:t:j比对的得分按下述方法求解,2020/5/29,63,其初值为：,fori=1,2,.,m,forj=1,2,.,n,2020/5/29,64,距离矩阵按照上述方法，对于给定的得分函数p，两个序列所有前缀的得分定义了一个(m+1)(n+1)的距离矩阵D=(di,j)其中di,j=S(0:s:i,0:t:j),di,j的计算公式如下：,2020/5/29,65,di,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),tjalignedtogap,sialignedtogap,2020/5/29,66,动态规划算法计算过程：计算过程从d0,0开始可以是按行计算，每行从左到右，也可以是按列计算，每列从上到下。当然，任何计算过程，只要满足在计算di,j时di-1,j、di-1,j-1、和di,j-1都已经被计算这个条件即可。在计算di,j后，需要保存di,j是从di-1,j、di-1,j-1、或di,j-1中的哪一个推进的，或保存计算的路径，以便于后续处理。上述计算过程到dm,n结束。,2020/5/29,67,最优路径求解：与计算过程相反从dm,n开始，反向前推。假设在反推时到达di,j，根据保存的计算路径判断di,j究竟是根据di-1,j、di-1,j-1、和di,j-1中的那一个计算而得到的。找到这个点以后，再从此点出发，一直到d0,0为止。走过的这条路径就是最优路径（即代价最小路径），其对应于两个序列的最优比对。,2020/5/29,68,例：s=AGCACACAt=ACACACTA,得分矩阵D（99）,p(a,a)=1p(a,b)=0abp(a,-)=p(-,b)=-1,2020/5/29,69,初始化,初始化函数：,2020/5/29,70,计算d（2，2）,di,j计算公式：,2020/5/29,71,最终的得分矩阵及序列比对,AGCACACA|ACACACTA,2020/5/29,72,作业（二）：1、标出矩阵中黄色封闭框内得分的计算途径（如黄色箭头的标注方法）；2、找出图中所标示出的最优途径的错误之处，以红色线（如图）表示正确的途径。,2020/5/29,73,Traceback,HEAGAWGHE-E-P-AW-HEAE,TracearrowsbackfromthelowerrighttotopleftDiagonalbothUpuppergapLeftlowergap,2020/5/29,74,序列长度的影响：令cw(s,t)表示两个长度分别为m和n的序列的相似性得分设cw(s,t)=99如果m=n=100-则可以说这两个序列非常相似但如果m=n=1000，则仅有10%相同相对长度的得分sim（s，t）=2*cw(s,t)/(m+n)算法分析：数据结构di,j空间复杂度：O(mn)时间复杂度：O(mn),2020/5/29,75,2、子序列与完整序列的比对,-AGCT-ATGCAGCTGCTT,2020/5/29,76,目标：使S(s,i:t:j)最大,序列S:,序列t:,ij,不计前缀0:t:i的得分,也不计删除后缀的j+1:t:|t|得分,2020/5/29,77,不计前缀0:t:i的得分处理第一行,2020/5/29,78,不计删除后缀的j+1:t:|t|得分处理最后一行,2020/5/29,79,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),不计代价,2020/5/29,80,距离矩阵初始化时，对第一行进行如下处理：d0,j=0for0jn最后一行的计算应该是：,同样，dm,n依然是最优局部比对的得分，而匹配的子列i:t:j按如下方式寻找：（1）j=minkdm,k=dm,n（2）反推比对路径，最终通过斜线（非空位）到达（0，i）。,2020/5/29,81,3、寻找最大的相似子序列,目标：使dw(i:s:j,i:t:j)