序列对比-药物生物信息学

1、2022/9/141药物生物信息学序列比对2022/9/142练习：Alignment:GACGGATTAG GATCGGAATA设有两个序列：GACGGATTAG，GATCGGAATAGA- CGGATTAGGATCGGAATA -2022/9/143序列比对的任务与目的序列比对的根本任务：发现序列之间的相似性辨别序列之间的差异目的：相似序列 相似的结构，相似的功能 判别序列之间的同源性推测序列之间的进化关系 2022/9/144本章内容 序列的相似性二重序列比对多重序列比对DNA片断组装2022/9/145第一节 序列的相似性字母表和序列编辑距离矩阵作图法序列的两两比对打分矩阵2022/9

2、/146同源与相似序列同源（homology）- 具有共同的祖先直向同源（Orthologous ）（不同种属） 共生同源（paralogous ）（同一种属）序列相似（similarity）同源序列一般是相似的 相似序列不一定是同源的 进化趋同（同功能）2022/9/147直向同源（a1 in species I, a1 in species II)共生同源（a1 and a2 in species I）进化趋同水平转移基因复制2022/9/148序列的相似性描述定性的描述序列的比对是一种关于序列相似性的定性描述，它反映在什么部位两条序列相似，在什么部位两条序列存在差别。最优比对揭示两条序列

3、的最大相似程度，指出序列之间的根本差异。定量的描述相似度距离GA- CGGATTAGGATCGGAATA -2022/9/149序列比对的基本操作是比对（Alignment） 两个序列的比对是指这两个序列中各个字符的一种一一对应关系，或字符的对比排列 。设有两个序列：GACGGATTAG，GATCGGAATAGAlignment2: GACGGATTAGGATCGGAATAGAlignment1:GACGGATTAG GATCGGAATAG2022/9/14101、字母表和序列将生物分子序列抽象为字符串，其中的字符取自特定的字母表。字母表4字符DNA字母表：A, C, G, T扩展的遗传学字母

4、表或IUPAC编码单字母氨基酸编码上述字母表形成的子集2022/9/1411符 号含 义说 明GGGuanine AAAdenine TTThymine CCCytosineRG or APurine YT or CPyrimidine MA or CAmino KG or TKeto SG or CStrong interaction (3 H bonds) WA or TWeak interaction (2 H bonds) HA or C or TNot-GBG or T or Cnot-AVG or C or Anot-T(not-U) DG or A or Tnot-C NG or

5、 A or T or CAny 扩展的遗传学字母表或IUPAC编码2022/9/1412氨基酸名称英文缩写简 写氨基酸名称英文缩写简 写甘氨酸GlyG丝氨酸SerS丙氨酸AlaA苏氨酸ThrT缬氨酸ValV天冬酰胺AsnN异亮氨酸IleI谷酰胺GlnQ亮氨酸LeuL酪氨酸TyrY苯丙氨酸PheF组氨酸HisH脯氨酸ProP天冬氨酸AspD甲硫氨酸MetM谷氨酸GluE色氨酸TrpW赖氨酸LysK半胱氨酸CysC精氨酸ArgR20种标准氨基酸的英文简写2022/9/1413 代表字母表 A* 代表由字母表A中字符所形成的一系列有限长度序列或字符串或序列的集合 a、b、c代表单独的字符 s、t、u

6、、v代表A*中的序列 |s|代表序列s的长度关于部分特定的符号的规定2022/9/1414为了说明序列s的子序列和s中的单个字符，在s中各字符之间用数字标明分割边界例如，设s=ACCACGTA，则s可表示为 0A1C2C3A4C5G6T7A8 i:s:j 指明第i位或第j位之间的子序列, 0 i j |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) 2022/9/1415子序列与子串s的子序列：选取s中的某些字符（或删除s中的某些字符）

7、而形成s的子序列例如： TTT 是 ATATAT的子序列。 s的子串：是由s中相继的字符所组成。 例如：TAC是AGTACA的子串， 但不是TTGAC的子串（是子序列）子串是子序列（连续子序列）子序列不一定是子串2022/9/1416字符串操作字符串连接操作：两个序列s和t的连接： s + + t例如：ACC+CTA = ACCCTA 字符串k操作 删除字符串两端的字符 其定义如下：prefix(s,l) = sk|s|-lsuffix(s,l) = k|s|-ls i:s:j = kisk|s|-j 2022/9/1417序列比对的四种基本应用：（1）两条长度相近的序列相似 找出序列的差别

8、两个实验室测定同一序列，实验结果比较（2）判断一条序列的前缀与另一条序列的后缀相似 大规模DNA测序中序列片断的组装 （3）判断一条序列是否是另一条序列的子序列 搜索特定模式（4）判断两条序列中是否有非常相似的子序列 分析保守序列2022/9/14182、编辑距离（Edit Distance)GCATGACGAATCAG TATGACAAACAGC GCATGACGAATCAG TATGAC-AAACAGC 说明两条序列的相似程度 定量计算 2022/9/1419 两条序列的相似程度的定量计算相似度，它是两个序列的函数，其值越大，表示两个序列越相似 两个序列之间的距离。距离越大，则两个序列的相

9、似度就越小 2022/9/1420海明距离：两条长度相等的序列，Hamming距离等于对应位置字符不同的个数。不足：（1）序列长度不同（2）未必反映两条序列的真正对应关系（3）DNA复制过程中的碱基插入、删除2022/9/1421字符编辑操作（Edit Operation）字符编辑操作可将一个序列转化为一个新序列 Match（a，a）字符匹配Delete（a，-） 1序列删除1个字符a，或2插入空位 Replace（a，b）2序列字符b替换1序列字符aInsert（-，b）1序列插入空位，或2删除1个字符b2022/9/1422扩展的编辑操作ACCGACAATATGCATA ATAGGTATA

10、ACAGTCAACCGACAATATGCATA ACTGACAATATGGATA 第二条序列头尾颠倒CTAGTCGAGGCAATCTGAACAGCTTCGTTAGT ？CTAGTCGAGGCAATCTCTTGTCGAAGCAATCA替换为互补碱基2022/9/1423反向互补序列RNA发夹式二级结构2022/9/14243、通过点矩阵进行序列比对“矩阵作图法” 或 “对角线作图” 2022/9/14252022/9/14262022/9/1427 序列1 序列2 实 例2022/9/1428 序列1 序列1 自我比对2022/9/1429滑动窗口技术两条序列中有很多匹配的字符对，因而在点矩阵中

11、会形成很多点标记。2022/9/1430滑动窗口技术 使用滑动窗口代替一次一个位点的比对是解决这个问题的有效方法。 假设窗口大小为10，相似度阈值为8，则每次比对取10个连续的字符，如相同的字符超过8个，则标记 基于滑动窗口的点矩阵方法可以明显地降低点阵图的噪声，并且明确无误的指示出了两条序列间具有显著相似性的区域。 2022/9/1431（a）对人类（Homo sapiens）与黑猩猩（Pongo pygmaeus）的球蛋白基因序列进行比对的完整点阵图。（b）利用滑动窗口对以上的两种球蛋白基因序列进行比对的点阵图，其中窗口大小为10个核苷酸，相似度阈值为8。 (a) (b) 2022/9/1

12、432作业（一）1、直接比较下面两个序列的相似度；并观察是否有更好的比对。2、尝试用矩阵标记图做出下面序列的最优比对1: ATTCGAGCCT2: CGTTCAGCTAATCGACTTGCATTCGAGCCT2022/9/14334、 序列的两两比对序列的两两比对（Pairwise Sequence Alignment）按字符位置重组两个序列，使得两个序列达到一样的长度匹配和替换，插入和删除2022/9/1434 s:AGCACACAAGCACACA t:ACACACTAACACACTA Match(A, A)Match(A, A)Delete(G, - )Replace(G, C)Match

13、(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 -1 Alignment -22022/9/1435不同编辑操作的代价不同 为编辑操作定义函数w，它表示“代价（cost）”或“权重（weight）”。对字母表中的任意字符a、b，定义 w (a, a) = 0 w

14、 (a, b) = 1 a b w (a, -) = w ( -, b) = 1 2022/9/1436使用得分（score）函数来评价编辑操作 p (a, a) = 1 p (a, b) = 0 a b p (a, -) = p ( -, b) = -1 2022/9/1437概念：两条序列s 和 t 的比对的得分（或代价）等于将s 转化为t 所用的所有编辑操作的得分（或代价）总和；s 和t 的最优比对是所有可能的比对中得分最高（或代价最小）的一个比对；s 和t 的真实距离应该是在得分函数p值（或代价函数w值）最优时的距离。 2022/9/1438例： 序列比对的目的是寻找一个得分最大（或代

15、价最小）的比对s:AGCACACAt:ACACACTA score (s，t）= 5s:AGCACACAt:ACACACTA cost=22022/9/14395、打分矩阵（Weight Matrices）不同类型的字符替换，其代价或得分是不一样的，特别是对于蛋白质序列。 某些氨基酸可以很容易地相互取代，而不改变理化性质。缬氨酸丙氨酸（得分高，代价小）赖氨酸丙氨酸（得分低，代价高）比较保守的替换比起较随机的替换更可能维持蛋白质的功能，且更不容易被淘汰不同打分方法对应不同的打分矩阵2022/9/1440(1) 核酸打分矩阵核酸打分矩阵设DNA序列所用的字母表为 = A，C，G，T a. 等价矩阵

16、 b. BLAST矩阵 （最流行的核酸序列比对程序）ATCGA1000T0100C0010G0001ATCGA5-4-4-4T-45-4-4C-4-45-4G-4-4-45 等价矩阵表 BLAST矩阵2022/9/1441腺嘌呤A尿嘧啶U胞嘧啶C鸟嘌呤G胸腺嘧啶Tc. 转移矩阵（转换transition，颠换transversion） （嘌呤：腺嘌呤A，鸟嘌呤G；嘧啶：胞嘧啶C，胸腺嘧啶T）ATCGA1-5-5-1T-51-1-5C-5-11-5G-1-5-51 转移矩阵2022/9/1442（2）蛋白质得分矩阵（i）等价矩阵（ii） 氨基酸突变代价矩阵GCM （iii）疏水矩阵 （iv）PA

17、M矩阵（Point Accepted Mutation）（v） BLOSUM矩阵 （Blocks Amino Acid Substitution Matrices）其中Rij代表打分矩阵元素i、j分别代表字母表第i和第j个字符。2022/9/1443氨基酸突变代价矩阵（遗传密码GCM）GCM通过计算一个氨基酸残基转变到另一个氨基酸残基所需的密码子变化数目而得到，矩阵元素的值对应于代价。GCM常用于进化距离的计算，其优点是计算结果可以直接用于绘制进化树。但在相似程度很低的蛋白质序列比对中很少使用。2022/9/14442022/9/1445疏水矩阵该矩阵是根据氨基酸残基替换前后疏水性的变化而得到

18、得分矩阵。若一次氨基酸替换疏水特性不发生太大的变化，则这种替换得分高，否则替换得分低。2022/9/14462022/9/1447PAM矩阵（Point Accepted Mutation）基于进化的点突变接受模型通过统计相似序列比对中的各种氨基酸替换发生率而得到该矩阵。1978年美国Dayhoff研究小组研究了71个相关蛋白质家族的1572个突变，发现蛋白质家族中氨基酸的替换并不是随机的。这意味着，在进化历程上，相关的蛋白质在某些位置上可以出现不同的氨基酸，这些点突变已经被进化所接受。2022/9/1448矩阵集合- PAM-N一个PAM就是一个进化的变异单位, 即1%的氨基酸改变如，PAM

19、120矩阵用于比对相距120个PAM单位的序列。一个PAM-N矩阵元素（i，j）的值： 反应两个相距N个PAM单位的序列中第i种氨基酸替换第j种氨基酸的频率。2022/9/1449针对不同的进化距离采用PAM 矩阵序列相似度 = 40% 50% 60% | | |打分矩阵 = PAM120 PAM80 PAM 60PAM250 14% - 27% 2022/9/1450 归一化打分实例2022/9/1451BLOSUM矩阵BLOcks SUbstitution is another amino acid substitution matrix, first calculated by Heni

20、koff and Henikoff (1992, PNAS)For its calculation only blocks of amino acid sequences with small change between them are considered. These blocks are called conserved blocks 2022/9/1452BLOSUM相似性度量BLOSUM与PAM矩阵阶数含义相反，低阶BLOSUM更多用来比较亲缘较远序列。 序列相似度 = 40% 62% 80% | | | 打分矩阵 = BLOSUM40 BLOSUM62 BLOSUM 60202

21、2/9/1453BLOSUM 622022/9/1454第二节 二重序列比对基本算法直接方法 生成两个序列所有可能的比对，分别计算代价函数，然后挑选一个代价最小的比对作为最终结果。本质问题：优化动态规划寻优策略动态规划算法（Dynamic Programming）1、序列两两比对基本算法2022/9/1455最短路径问题起点终点C1 C2 W1 W2路径1：C1 + w1 ？路径2：C2 + w2 ？ 取最小值！算法求解: 从起点到终点逐层计算2022/9/14562022/9/1457求解过程起点终点ATTCCGAAGA AGTCGAAGGT 从两个序列前端开始 逐步推进 直到两个序列的末端

22、。2022/9/1458序列S: i-1 i i+1序列t: j-1 j j+1序列S: i-1 i i+1序列t: j-1 j j+1Case1:匹配（si，tj )中间过程：比对0:S:i 与 0:T:j 2022/9/1459序列S: i-1 i i+1序列t: j-1 j j+1序列S: i-1 i i+1序列t: j-1 j j+1Case2:删除（si， -） 比对0:S:i 与 0:T:j2022/9/1460序列S: i-1 i i+1序列t: j-1 j j+1序列S: i-1 i i+1序列t: j-1 j j+1Case3:插入（ -，tj ） 比对0:S:i 与 0:T

23、:j2022/9/1461设序列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 ）。 2022/9/1462令：为序列0:s:i和与序列 0:t:j

24、比对的得分按下述方法求解 2022/9/1463其初值为：for i=1 , 2 ,., mfor j=1 , 2 ,., n2022/9/1464距离矩阵按照上述方法，对于给定的得分函数p，两个序列所有前缀的得分定义了一个(m+1)(n+1)的距离矩阵D = ( d i , j )其中d i , j = S (0:s:i , 0:t:j ) d i , j的计算公式如下：2022/9/1465d i , j 最小值的三种选择决定了各矩阵元素之间的关系，用下图表示：di,jdi,j-1di-1,jdi-1,j-1 距离矩阵元素d i , j 的计算S (0:s:i , 0:t:j )S (0:

25、s:i-1 , 0:t:j )S (0:s:i-1 , 0:t:j-1 )S (0:s:i , 0:t:j-1 )tj aligned to gapsi aligned to gap2022/9/1466动态规划算法计算过程： 计算过程从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中的哪一个推进的，或保存

26、计算的路径，以便于后续处理。上述计算过程到d m , n结束。 2022/9/1467最优路径求解：与计算过程相反 从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为止。 走过的这条路径就是最优路径（即代价最小路径），其对应于两个序列的最优比对。 2022/9/1468tsgapACACACTAgapAGCACACA例：s = AGCACACAt = ACACACTA 得分矩阵D （99）p (a,

27、 a) = 1p (a, b) = 0 a bp (a, -) = p ( -, b) = -12022/9/1469tsgapACACACTAgap0-1-2-3-4-5-6-7-8A-1G-2C-3A-4C-5A-6C-7A-8初始化初始化函数：2022/9/1470计算d（2，2）tsACACACTA0-1-2-3-4-5-6-7-8A-110-1-2-3-4-5-6G-201C-3A-4C-5A-6C-7A-8di,j计算公式：2022/9/1471最终的得分矩阵及序列比对tsACACACTA0-1-2-3-4-5-6-7-8A-110-1-2-3-4-5-6G-2010-1-2-3-

28、4-5C-3-11110-1-2-3A-4-2021210-1C-5-3-1132321A-6-4-2024333C-7-5-3-113543A-8-6-4-202455AGCACACA| | |ACACACTA2022/9/1472tsACACACTA0-1-2-3-4-5-6-7-8A-110-1-2-3-4-5-6G-2010-1-2-3-4-5C-3-11110-1-2-3A-4-2021210-1C-5-3-1132321A-6-4-2024333C-7-5-3-113543A-8-6-4-202455作业（二）：1、标出矩阵中黄色封闭框内得分的计算途径（如黄色箭头的标注方法）；2、

29、找出图中所标示出的最优途径的错误之处，以红色线（如图）表示正确的途径。2022/9/1473TracebackHEAGAWGHEE0-8-16-24-32-40-48-56-64-72-80P-8-2-9-17-25-33-42-49-57-65-73A-16-10-3-4-12-20-28-36-44-52-60W-24-18-11-6-7-15-5-13-21-29-37H-32-14-18-13-8-9-13-7-3-11-19E-40-22-8-16-16-9-12-15-73-5A-48-30-16-3-11-11-12-12-15-52E-56-38-24-11-6-12-14-1

30、5-12-91HEAGAWGHE-E-P-AW-HEAETrace arrows back from the lower right to top left Diagonal both Up upper gap Left lower gap2022/9/1474序列长度的影响：令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) 算法分析：数据结构d i , j空间复杂度：O (mn)时间复杂度：

31、O (mn)2022/9/14752、子序列与完整序列的比对-AGCT-ATGCAGCTGCTT2022/9/1476目标：使S(s, i:t:j ) 最大序列S:序列t: i j不计前缀0:t:i 的得分, 也不计删除后缀的j+1:t:|t|得分2022/9/1477不计前缀0:t:i 的得分处理第一行tsACACACTA000000000AGCACACA2022/9/1478不计删除后缀的j+1:t:|t|得分 处理最后一行2022/9/1479dm,jdm,j-1dm-1,jdm-1,j-1S (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 )不计代价2022/9/1480距离矩阵初始化时，对第一行进行如下处理：d0,j = 0 for 0 j n 最后一