生物序列比对中的算法课件

生物序列比对中的算法中科院计算所生物信息学实验室华大—曙光联合实验室张法提纲背景知识序列相似性的比较两条序列的比对问题多序列的比对问题一些启发式的算法生物序列比对中的并行算法DNA(1)脱氧核糖核酸DNA的分子组成核甘(nucleotides)磷酸盐(phosphate)糖(sugar)一种碱基腺嘌呤(Adenine)鸟嘌呤(Guanine)胞嘧啶(Cytosine)胸腺嘧啶(Thymine)DNA(2)碱基的配对原则A(腺嘌呤)—T(胸腺嘧啶)C(鸟嘌呤)—G(胞嘧啶)一个嘌呤基与一个嘧啶基通过氢键联结成一个碱基对。DNA分子的方向性5'→3'DNA(3)DNA的双螺旋结构

碱基对之间的互补能力DNA(4)DNA的复制在DNA解旋酶的作用下两条链分离开，分别作为一个模板，在聚合酶的作用下合成一条新链。RNA、转录和翻译RNA(核糖核酸)：单链结构、尿嘧啶U代替胸腺嘧啶T、位于细胞核和细胞质中。转录：DNA链→RNA链信使RNA(mRNA)，启动子。翻译：mRNA上携带遗传信息在核糖体中合成蛋白质的过程。变异进化过程中由于不正确的复制，使DNA内容发生局部的改变。变异的种类主要有以下三种：替代(substitution)插入或删除(insertionordeletion)

indel重排(rearrangement)蛋白质由氨基酸依次链接形成在生物体中总共有20种氨基酸。蛋白有十分复杂的三维结构。其三维机构决定了蛋白质的功能。基因什么是基因？DNA上具有特定功能的一个片断，负责一种特定性状的表达。一般来讲，一个基因只编码一个蛋白质。基因组任何一条染色体上都带有许多基因，一条高等生物的染色体上可能带有成千上万个基因，一个细胞中的全部基因序列及其间隔序列统称为genomes（基因组）。DNA上的基因

基因基因的编码基因编码是一个逻辑的映射，表明存储在DNA和mRNA中的基因信息决定什么样的蛋白质序列。每个碱基三元组称为一个密码子(codon)碱基组成的三元组的排列共有43＝64种，而氨基酸共有20种类型，所以不同的密码子可能表示同一种氨基酸。

带来的问题序列排列问题基因组的重排问题蛋白质结构和功能的预测基因(外显子、内含子)查找问题序列装配(SequenceAssembly)问题动机在生物学的研究中，将未知序列同已知序列进行比较分析已经成为一种强有力的研究手段，生物序列相似性比较中绝大部分的问题在计算机科学领域中主要体现为字符串的匹配和查找。图谱分析生物序列片断的拼接基因区域的预测生物序列同源的比较蛋白质结构和功能预测两条序列比对问题的分类全局比对(GlobalAlignment)局部比对(LocalAlignment)空位罚分(GapPenalty)全局比对(1)－定义定义1：两个任意的字符x和y，(x,y)表示表x和y比较时的分值。

(x,x)=2,(x,y)=(x,-)=(-,y)=-1定义2：S=s1…sn和T=t1…tm，其全局比对A可以用序列S´和T´来表示，其中：(1)|S´

|=|T´|；

(2)将S´和T´中的空字符除去后所得到的序列分别为S和T；比对A的分值Score为：全局比对(2)－原始算法输入：序列S和T，其中|S|=|T|=n

输出：S和T的最优比对fori=0tondofor(S的所有的子序列A，其中|A|=i)dofor(T的所有的子序列B，其中|B|=i)do

……全局比对(3)动态规划(DynamicProgramming)Bellman在50年代提出的理论基础－最优化原理弱点：thecurseofdimensionalitySmith-WatermanAlgorithm全局比对(4)Smith-Waterman算法计算出两个序列的相似分值，存于一个矩阵中。(editmatrix、DP矩阵)根据此矩阵，按照回溯的方法寻找最优的比对序列。全局比对(5)前提条件递归关系计算editmatrix：fori=0tondoforj=0tomdoCalculateV(i,j)usingV(i-1,j-1),V(i,j-1),V(i-1,j)end全局比对(6)计算完editmatrix后，得到序列比对的最优分值，通过回溯（Traceback）的方法可获得序列的最优比对序列。回溯的算法：for(i=|S|，j=|T|;i>0&&j>0;){if(V[i,j]=V[i–1,j–1]+σ(S[i],T[j])){i––;j––;}elseif(V[i,j]=V[i–1,j]+σ(S[i],–)){i––;insert(‘–’,T,j);}elseif(V[i,j]=V[i,j–1]+σ(–,T[i])){j––;insert(‘–’,S,i);}}例：S=“acgctg”和T=“catgt”

(x,x)=2,(x,y)=(x,-)=(-,y)=-1

ji01c2a3t4g5t00-1-2-3-4-51a-1-110-1-22c-2100-1-23g-300-1214c-4-1-1-1115t-5-2-21036g-6-3-3032三种可能的最优比对序列：S:acgctg-T:-c–atgtS:acgctg-T:-ca–tgtS:-acgctgT:catg-t-

局部比对(1)两条序列在一些局部的区域内具有很高的相似度。在生物学中局部比对比全局比对更具有实际的意义。局部比对(2)前提条件：

V(i,0)=0;

V(0,j)=0;递归关系：

找出i*和j*，使得：

局部比对(3)对全局比对策略稍作修改可得到局部最优比对算法。比对的路径不需要到达搜索图的尽头,如果某种比对的分值不会因为增加比对的数量而增加时，这种比对就是最佳的。依赖于记分系统的性质：因为某种路径的记分会在不匹配的序列段减少,当分值降为零时，路径的延展将会终止，一个新的路径就会产生。局部比对(4)S=“abcxdex”，T=“xxxcde”记分函数(x,y)：

(x,x)=2,(x,y)=(x,-)=(-,y)=-1

ji01x2x3x4c5d6e000000001a00000002b00000003c00002104x02221105d01111326e00000257x0222114S=“abcxdex”，T=“xxxcde”

局部最优比对是：

cxde

c-de或

x-de

xcde

空位罚分(1)空位:序列中任意连续的尽可能长的空格。空位的引入是为了补偿那些插入或缺失，但是在序列的比对中引入的空位不能太多，否则会使序列的排列变得面目全非。每引入一个空位，比对的分值都会有所扣除，常见的罚分规则主要有两种：权值恒定模型仿射罚分模型空位罚分(2)权值恒定模型：在每个空位中的空格的分值为零，即：(x,-)=(-,y)=0

比对的分值为：其中，Wg为开放一个空位的罚分空位罚分(3)仿射罚分模型(AffineGapModel)

用一个附加的罚分比例去乘空位的长度Wg+q×Ws

比对的分值为：空位罚分(4)初始条件：

V(0,0)=0;

V(i,0)=E(i,0)=Wg

+iWs;

V(0,j)=F(0,j)=Wg

+jWs;递归条件：

V(i,j)=max{G(i,j),E(i,j),F(i,j)};

G(i,j)=V(i-1,j-1)+σ(S[i],T[j]);

E(i,j)=max{E(i,j-1)+Ws,V(i,j-1)+Wg+Ws}

F(i,j)=max{F(i-1,j)+Ws,V(i-1,j)+Wg+Ws}.利用动态规划计算序列最优比对的算法的复杂度分析：时间复杂度；O(mn)

多序列比对(1)两条序列比对问题的一般化推广动机：携带更多的消息。DNA或蛋白质数据库容量的指数级的增长，即使使用很简单的记分函数也很难找到一种在有效时间内的解决方法，而几乎所有的近似算法和许多的启发式算法，实际上都是在算法的计算速度和获得最佳比对结果的敏感性之间寻找一种权衡策略。多序列比对(2)rigorous算法两条序列多条序列，NPhardtree_based算法和iterative算法首先从序列中找出两条相似度最大的比对，然后按照相似度递减的顺序连续添加一些序列到最优的比对序列中。序列非常接近原始算法基础上的近似算法，但是它们也是非常耗时的。多序列比对(3)记分方法：用函数δ(x,y)来计算字符x和y之间的距离，两个序列的比对的距离分值我们用V来表示：

k条序列比对的分值：为k条序列中任意两条序列(共有Ck2条)的分值(距离)V之和，用SP(Sum_of_Pairs)来表示:动态规划的方法在多序列比对中的时间复杂度为O(2n)k中心星算法输入：一组含有k条序列的集合Ω

找出序列St，St∈Ω，使得∑i≠tD(

Si,St)的值最小，令A={St}逐次地添加Si∈Ω-{St}到A中，并使Si与St的比对的值最小；假设S1，S2，…Si-1已添加到A中，由于在分别和St进行比对的过程中需要加入一些空格，故此时A为：A={S1’，S2’，…Si-1’，S’t}。添加Si到A中，按照两条序列比对的动态规划算法比较S’t和Si，分别产生新的序列St’’和Si’，再按照St’’中添加空格的位置调节序列S1’，S2’，…Si+1’为S1’’，S2’’，…Si-1’’，并用St’’替换St。启发式算法-FASTA基本思想是：一个能够揭示出真实的序列关系的比对至少包含一个两个序列都拥有的字(片断)，把查询序列中的所用字编成索引，然后在数据库搜索时查询这些索引，以检索出可能的匹配，这样那些命中的字很快被鉴定出来。

确定参数ktup，在两个序列中查找长度为ktup的、相匹配的片段(增强点)。为了提高速度，可以通过查询表格或hash表来完成，然后在表格中搜索与另一条序列相匹配的、长度为ktup的片段。在同一条对角线中临近的增强点成为一个增强段。每一个增强点都赋予一个正的分值，一个增强段中相邻的两个增强点之间的不匹配区域赋予一定的负值。一个增强段对应于一段相匹配的子序列，分值最高的段被标记为init1。

引入indel。把那些没有重叠(non-overlap)的增强段拼接起来(增强段的分值之和减去空位处罚)。分值最高的区域记为initn。对最有可能的匹配序列进一步评分：以增强段init1所在的对角线为中心，划分出一个较狭窄的对角线带，利用S-W算法，来获得分值最高的局部比对，记作opt。决定采用initn或opt的分值，前者敏感度低但速度快。FASTA对每一个检索到的比对都提供一个统计学显著性的评估，以判断该比对的意义。FASTA对DNA序列搜索的结果要比对蛋白质序列搜索的结果更敏感。它对数据库的每一次搜索都只有一个最佳的比对，一些有意义的比对可能被错过。

启发式算法-BLAST基本思想是：通过产生数量更少的但质量更好的增强点来提高速度。BALST算法是建立在严格的统计学的基础之上的。它集中于发现具有较高的相似性的局部比对，且局部比对中不能含有空位。由于局部比对的限制条件，在大多数情况下比对会被分解为若干个明显的HSP(High-scoreSequencePairs)。首先确定一个终止值S、步长参数w和一个阈值t，S值通常是基于统计学的原理指明一个预期的终止E值，然后软件会在考虑搜索背景性质的基础上计算出合适的S值。使要比对的序列中包含一个分值不小于S的HSP。引入邻近字串的思想：不需要字串确切地匹配，当有一个字串的分值高于t时，BALST就宣称找到了一个选中的字串。为了提高速度，允许较长的字串长度W。W值很少变化，这样，t值就成为权衡速度和敏感度的参数。一个字串选中后，程序会进行没有空位的局部寻优，比对的最低分值是S，当比对延伸时会遇到一些负的分值，使得比对的分值下降，当下降的分值小于S时，命中的延伸就会终止。这样系统会减少消耗于毫无指望的选中延伸的时间，使系统的性能得以改进。BLAST的改进在1997年提出了对BLAST程序的改进算法，提高了搜索速度、敏感度和实用性。可处理间隔(gap)的gappedBLAST算法PSI-BLAST算法对一个选中字串长度标准的延伸利用profile(表头文件)的数据结构来进行搜索算法算法扩大步长，以步长为2w来搜索。允许位于不同的对角线的两个片段拼接在一起。位于不同对角线的两个片段拼接在一起的前提条件是：拼接后片段的分值不小于某一个终止值。执行通常的BLAST算法，使用一种不同的记分方式，根据高度显著比对(HSPs)的最高分值建立一个最初的profile。根据该profile反复利用BLAST算法对数据库进行搜索，这一步实际上是根据表头文件的统计结果扩展局部比对。这一过程是反复进行的，直到再没有发现新的有意义的匹配为止。由于在每一轮都会有新的片段加入，因此在操作过程中profile需要在每一个循环结束之后更新。为适应同源比较的不通情况，BLAST得到了各种各样的补充和发展，形成了一个程序家族。程序名检索序列数据库序列BLASTP蛋白质蛋白质BLASTN核苷酸核苷酸BLASTX核苷酸蛋白质TBLASTN蛋白质核苷酸

SubstitutionMatrixPAM250BLOSUM62在SUN10000上进行BLAST的测试：8CPU8GmemoryVERSION:BLASTN2.2.1[Jul-12-2001]blastall-pblastn-i$query.seq–d/mnt/database/public/NCBI/ftp.ncbi.nlm.nih.gov/BLAST/nt/2001-12-02/nt-e1e-10-o$result.out-a4>rzy_63001.y1.abdCHROMAT_FILE:rzy_63001.y1.abdPHD_FILE:rzy_63001.y1.abd.phd.1CHEM:termDYE:ETTIME:ThuJun2116:03:572001TAAAGGAGAATGACTCTTTTAAGTCATAGCTTGCTGNCATTAGTCGTTTGAAAGGTCGAGTGCTCAATATCTTAACCCTTTTACCTTTTTTGAAAAACACTAGGCTAGATAGATGCCATATCGGATTGAGGGGTTATTATTGGAGATTTGGAGGGGCATTTTTTTTTAGTTATTTTTTTTCTAAATTTTTTGGATGGTGAGGCTATATGAGGTGCTCTTGAAGATGCTCTCCCTTTATGTCCGAGTGATTGAATGGTATAATTTGCAACACACAAGCTTGCTTTTACAGAACTCATGGACTTCTTTACTTATATCGTATAATTCGCAGGAATATTTCGCGAAATGGGGAGAGCAAGGCACGGTGGATCTGAAACGCGAGCTCGACCTCCTCATCCTGACGATCGCGAGCCGCGTCCTCCTCGGCAAGGAGGTTAGGGAGACCATGTTCGCCGACGTCGTGGCATCCTTCCACGAGCTCATGGACAACAGCATGCACCTCATCAGCCTCTGCTTCCCCAACCTCCCGATCCCGCGGCACCGCCGGCGCGACACGGCGTCCGCTAGGCTCAAGGAGCTCTTCTCCCGTGCCATCCAGCTTCGCCGCGGCTCCGGCCGCGCCGAGGACGACGTGCTCCAGCGGTTCTTGGATCCAGGTACAGGGACGGCntdatabase:

mnt/database/public/NCBI/ftp.ncbi.nlm.nih.gov/BLAST/nt/2001-12-02/ntSize:4.8G1´030´000个reads

借助于一些比最初的动态规划算法更快的启发式算法。借助硬件来完成动态规划的算法采用高性能计算系统，把问题有效地分配到多个处理器上运行，最后再把各处理器的运算结果收集起来生物序列比对中的

并行算法两条序列比对的并行算法根据序列的相似性比较，找出两者的最佳匹配找出从一条序列转化到另一条序列的最佳路径核心：动态规划动态规划中的数据相关Lander的处理方法如果已知第i-1和i-2行反对角线上的各项值，那么第i行反对角线上各项的值可以并行计算。△☆△☆△☆△☆☆??????D12D21D31D11D22D13Dm,n●●●处理器01m时间步12m+n-2m+n-1Smith-Waterman的并行处理Row-Wavefront算法和DiagonalWavefront算法

Row-Wavefront算法是让每个处理器顺序处理动态规划矩阵中的每一行，当一个处理器检测到它所需要的上一行矩阵的元素还没有计算出来时，该处理器就阻塞自己，直到所需要的元素计算出来。

处理器1处理器2处理器3序列S序列TSmith-Waterman的并行处理Row-Wavefront算法和DiagonalWavefront算法

DiagonalWavefront算法中所有的处理器同时沿着矩阵的一条反对角线进行计算，当一条反对角线的元素全部计算完，才转到下一条反对角线计算。当一个处理器空闲的时候，它从当前的反对角线上寻找一个还没有计算的元素进行计算。这种算法没有严格的处理器计算顺序序列T序列SHuang的处理方法采用了分而治之的策略对回溯算法进行了改进:按照中间的一条反对角线来分割动态规划矩阵。(0,0)r(