1、生物信息学,第三章 序列比对,为什么要序列比对？,寻找进化过程中的同源序列; 基于同源物鉴定的功能预测; 基本假设：序列的保守性 功能的保守性注意： 1. 蛋白质一般在三级结构的层面上执行功能； 2. 蛋白质序列的保守性决定于其编码DNA的保守性；,通常,本章内容提要,第一节：数学基础：概率及概率模型 第二节：双序列比对算法的介绍 Dot matrix 动态规划算法 (Needleman-Wunsch, Smith-Waterman算法) FASTA和BLAST算法 第三节：打分矩阵及其含义 第四节：多序列比对,第一节 序列比对的数学基础,排列组合,从N个物品中取出k个物品的排列数：从N个物品
2、中取出k个物品的组合数：,概率模型,概率模型:一个能够通过不同的概率产生不同结果的模型。概率模型可以模拟或者仿真某一类型的所有事件，并且对每个事件赋予一个概率。,色子模型：一个色子存在6个概率值：p1,p2, p6，其中，掷出i的概率为pi(i=1,2, ,6)。因此： pi0，且,考虑三次连续的掷色子，结果为1，6，3，则总概率为：p1p6p3,概率分布,考虑连续变量x，例如：物体的重量。则当重量确切为1公斤时的概率，为0。 变量的区间：P(x0xx1)当区间无限小 - 0时，上式： P(x -x/2 xx + x/2 ) = f(x)xf(x)称为概率密度函数 因此：,且,二项分布,1.
3、事件只有两种可能出现的结果。例如掷硬币，正面记为“1”，反面记为“0”。 2. 则掷硬币N次，有k次是1的概率为：,二项分布的期望值与标准方差,期望值 E(x) = ,方差 Var X=2,泊松分布 (Poisson distribution),1. 稀有事件发生的概率：在一个连续的时间或空间中，稀有离散变量出现的概率 2. N - , E(x)= ,e = 2.71828,泊松分布与二项分布的近似,对于大的N及小的p值的二项分布，能够相当准确地用一个参数为=Np的泊松分布近似。 当实验次数很多而概率很小时：二项分布泊松分布,例1：鸟枪法的覆盖率,假设：需要测序的BAC长度200kbp; 总共
4、测序的序列数量：N; 每次测序：500bp； 每次测序的覆盖率 p：500/200kbp=0.0025 因此：每个点平均覆盖到的次数: =N * p k: 测序能够覆盖到点X的次数。,鸟枪法：覆盖率,点X被覆盖k次的概率：(二项分布泊松分布),当点X一次都不被覆盖时，k=0; 此时的概率为：,覆盖率 vs. 准确性,例2：泊松分布,Prof. Gene发现一种序列上的调控信号，在人的基因组上平均每500kbp一个。那么，随机给一条1mbp的序列，在上面发现5个这样的信号，完全是随机产生的概率是多少？,本例中， E(x)= 2 (1mbp/500kbp),统计显著性：p-value 0.05,超
5、几何分布,与二项式分布的区别：不放回抽样。 例：有N个球，其中红球M个，白球N-M个，每次拿出一个球再放回，总共n次，其中有m个球是红球的概率为 (二项式分布)：,p=M/N,超几何分布 (2),上例改为：有N个球，其中红球M个，白球N-M个，每次拿出一个球不放回，总共n次，其中有m个球是红球的概率为：,并且，0mMN,超几何分布右尾概率,上例再改为：有N个球，其中红球M个，白球N-M个，每次拿出一个球不放回，总共n次，其中至少有m个球是红球的概率为：,并且，0mMN,超几何分布左尾概率,上例再改为：有N个球，其中红球M个，白球N-M个，每次拿出一个球不放回，总共n次，其中最多有m个球是红球的
6、概率为：,并且，0mMN,超几何分布双尾概率,所有出现概率 = 观察表概率的概率之和,Fishers Exact Test,超几何分布的精确概率计算。 前提是固定边际分布，即a+b 、c+d、a+c与b+d的值不变。 RA Fisher, 1935年文章示例：,Fishers Exact Test,计算公式：,=,统计显著性 假设检验中的P 值 (P value),P value：一种在原假设为真的前提下出现观察样本以及更极端情况的概率。 显著性水平A：认为预先设定的显著性水平阈值， P A 为显著。 一般以P 0.05 为显著， P 0.01 为非常显著，其含义是样本间的差异由抽样误差所致的
7、概率小于0.05 或0.01。,假设检验,本例中， 零假设H0：该女同事只是随便乱猜答案； 备择假设H1：该女同事所言虚；P value计算：P (a = 3 | a+b=c+d=a+c=b+d=4 ) = 0.229P (a = 4 | a+b=c+d=a+c=b+d=4 ) = 0.014,例 3：超几何分布,Prof. Gene从人的26873个蛋白质中预测了2264个能结合某类金属离子X。现已知，人的26873个蛋白质中有421个蛋白质具有某种功能结构域D，而在预测的2264个X金属蛋白中，有94个具有结构域D。问：结构域D在2264个X金属蛋白中是显著出现，显著不出现，还是随机出现？
8、,例 3：超几何分布,问题转化：在26873个蛋白质的体系中，取出2264个蛋白质，其中至少有94个蛋白质具有功能结构域D的概率是多少？ N = 26873; n= 2264; M = 421; m = 94;,例 3：超几何分布,例 3：超几何分布,a+b+c+d=26873 c+d=2264 b+d=421 d=94,http:/ Dot Matrix，点阵法 2. 动态规划算法： Global: Needleman-Wunsch Local: Smith-Waterman 3. Word or k-tuple算法：FASTA, BLAST,1. 点阵法,1970年，Gibbs & McI
9、ntyre； 寻找两条序列间所有可能的比对； 发现蛋白质或者DNA序列上正向或者反向的重复； 发现RNA上可能存在的互补区域。 工具： http:/myhits.isb-sib.ch/cgi-bin/dotlet http:/arbl.cvmbs.colostate.edu/molkit/dnadot/,点阵法：自身的比对,点阵法：重复序列,点阵法：反向重复/回文,点阵法：不同序列的比对,1:PKDFCKALV,2:PK-FTKAIV,Seq 1,Seq2,点阵法的序列比对,计算效率,用CPU的计算时间和内存占用量来衡量； 对于需要解决的问题，其单位数量n在某算法下运算的基本操作重复执行次数表
10、示为f(n)； 时间复杂度: T(n)=O(f(n)； 如果需要解决的问题的大小与单位数量n的平方成正比，则O(n2) 对于算法来说： O(1) O(log(n) O(n) O(n2) O(an) O(n!),NP问题,1. 一般的，O(nk), 当k3 时，为多项式时间，较为容易处理。 2. 当O(an)，则难以处理。 3. NP完全问题（NPC）：无法找到能够在多项式时间复杂度内解决方法的问题； 4. 近似算法/优化算法，求近似解。,P/NP问题-千禧年大奖难题之一,1900年，德国数学家David Hilbert提出的23个历史性数学难题。 千禧年大奖难题美国克雷数学研究所（Clay M
11、athematics Institute,CMI）于2000年5月公布七个世界数学难题。,P/NP/NPC问题,P问题： Polynomial Problems 可以在多项式( polynomial )时间内解决的问题; NP: “Non-deterministic Polynomial”,并非 “Non-Polynomial”可以在多项式的时间里验证一个解的问题; NPC: NP-complete,2. 动态规划算法,1. 打分模型、替代矩阵以及空位罚分。 2. 比对算法：递归及动态规划算法； 3. 全局优化比对：Needleman-Wunsch 4. 局部优化比对：Smith-Waterm
12、an 5. 工具资源： http:/www.ludwig.edu.au/course/lectures2005/Likic.pdf http:/ http:/zhanglab.ccmb.med.umich.edu/NW-align/ http:/ 字符相同：identity 2. 字符替代：similarity，相似性，氨基酸/碱基之间的替代和突变 3. 插入和缺失 4. 空位罚分,BLOSUM62替代矩阵,空位罚分,1. 线性罚分：d, 每次罚分的分数；g，空位数,2. 修正的罚分：d, 第一次罚分的分数； g，空位数；e, 修正后的参数,递归和动态规划算法,两条序列的比较，无空位：时间复杂
13、度为O(n2);两条序列比对，允许空位，时间复杂度为：因此，有空位的双序列比对，时间复杂度为：O(22n)，指数增加，NPC问题！,递归和动态规划算法 (2),数学上保证提供最优解。动态规划算法：比较所有可能的字符对，考虑匹配、错配以及空位罚分，并且将比对次数控制在多项式时间内。替代矩阵：BLOSUM62， 空位罚分：11 延伸的空位罚分：1 (BLAST工具),例：全局比对,序列1： V D S C Y 序列2： V E S L C Y 替代矩阵中的分数： 4 2 4 -11 9 7,两序列比对的总分： Score=(AA pair scores) gap penalty = 15,动态规划
14、算法：全局比对,本例：线性罚分,全局比对 (2),要求解Sij的分数，我们必须先知道Si-1, j-1, Si-1, j, 以及Si, j-1的分数，这种方法叫做递归算法；采用这种方法，可以把大的问题分割成小的问题逐一解决，即动态规划算法；需要存储如何得到Sij分数的过程。,全局比对 (3),i,j,Needleman-Wunsch算法； Sij = max of Si-1, j-1 + (xi, yj)Si-1, j - d (从左到右) Si, j-1 - d (从上到下),全局比对 (4),i,j,4,-11,-11,Needleman-Wunsch算法； Sij = max of Si
15、-1, j-1 + (xi, yj)Si-1, j - d (从左到右) Si, j-1 - d (从上到下),全局比对 (5),4,-11,-11,全局比对 (6),-3,-11,-11,Needleman-Wunsch算法； Sij = max of Si-1, j-1 + (xi, yj)Si-1, j - d (从左到右) Si, j-1 - d (从上到下),全局比对 (7),-3,-11,-11,全局比对 (8),4,2,回溯：比对结果,4,2,V D S C Y V E S L C Y,局部优化比对,下例：局部优化打分 两条序列如下：,L D S C H G E S L C K,
16、目标：使用局部优化算法寻找比对的结果,局部优化比对 (2),1. Smith-Waterman算法； 2. 时间复杂度O(n2)； 3. Sij = max of 0Si-1,j-1 + (xi, yj)Si-1,j - d (从左到右) Si,j-1 - d (从上到下) 本例中：gap: 12，线性罚分模型。,局部优化比对 (3),Smith-Waterman算法； Sij = max of Si-1, j-1 + (xi, yj)Si-1, j - d (从左到右) Si, j-1 - d (从上到下)0,局部优化比对 (4),-12,-12,-3,局部优化比对 (5),-12,-12,-4,局部优化比对 (6),-12,-2,局部比对结果：,L D S C H G E S L C K,打分有何不同？,1.局部优化比对：9分 2.全局优化比对：2+4-12+9=3 3. 为何不同？,L D S C H G E S L C K,http:/zhanglab.ccmb.med.umich.edu/NW-align/,
