第三章序列比较

第三章序列比较序列比较是生物信息学中最基本、最重要的操作。序列比较的根本任务是：通过比较生物分子序列，发现它们的相似性，找出序列之间共同的区域，同时辨别序列之间的差异。在分子生物学中，DNA或蛋白质的相似性是多方面的，可能是核酸或氨基酸序列的相似，可能是结构的相似，也可能是功能的相似。一个普遍的规律是序列决左结构，结构决左功能°研究序列相似性的目的之一是，通过相似的序列得到相似的结构或相似的功能。这种方法在大多数情况下是成功的，当然也存在着这样的情况，即两个序列几乎没有相似之处，但分子却折叠成相同的空间形状，并具有相同的功能。这里先不考虑空间结构或功能的相似性，仅研究序列的相似性。研究序列相似性的另一个目的是通过序列的相似性，判别序列之间的同源性，推测序列之间的进化关系。这里将序列看成由基本字符组成的字符串，无论是核酸序列，还是蛋白质序列，都是特殊的字符串•本章着重介绍通用的序列比较方法。第一节序列的相似性序列的相似性可以是泄量的数值，也可以是左性的描述。相似度是一个数值，反应两个序列的相似程度。关于两条序列之间的关系，有许多家词，如相同.相似、同源、同功、直向同源、共生同源等。在进行序列比较时经常使用“同源”(homology)和“相似”(similarity)这两个概念，这是经常容易被混淆的两个不同的概念。两个序列同源是指它们具有共同的祖先，在这个意义上无所谓同源的程度，两个序列要么同源，要么不同源。而相似则是有程度的差别，如两个序列的相似程度达到30%或60%。一般来说，相似性很高的两个序列往往具有同源关系。但也有例外，即两个序列的相似性程度很髙，但它们可能并不是同源序列，这两个序列的相似性可能是由随机因素所产生的，这在进化上称为“趋同”(convergence),这样一对序列可称为同功序列。直向同源序列来自于不同的种属，而共生同源序列则是来自于同一种属序列，其产生是由于进化过程中的序列复制。序列比较的基本操作是比对(Align)°两个序列的比对(Alignment)是指这两个序列中各个字符的一种一一对应关系，或字符对比排列。序列的比对是一种关于序列相似性的泄性描述，它反应在什么部位两个序列相似，在什么部位两个序列存在差别。最优比对揭示两个序列的最大相似程度，指出序列之间的根本差异。K字母表和序列在生物分子信息处理过程中，将生物分子序列抽象为字符串，英中的字符取自特左的字母表。字母表是一组符号或字符，字母表中的元素组成序列。一些重要的字母表有：4字符DNA字母表：{A,C,G,T}°扩展的遗传学字母表或IUPAC编码(见表2・3)°单字母氨基酸编码(见2-1)o上述字母表形成的子集。下而所讨论的内容独立于特左的字母表。首先规泄一些特立的符号：•A-代表字母表A*—代表由字母表A中字符所形成的一系列有限长度序列或字符串或序列的集合•a.b、c—代表单独的字符s、I、11、v—代表A*中的序列Isl—代表序列s的长度为了说明序列s子序列和s中单个字符，在s中各字符之间用数字标明分割边界。例如，设s二ACCACGTA,贝ljs可表示为oA1C2C3A4C5G6T7Agi：s：j指明第i位或第j位之间的子序列。当然，0<i<j<lsL子序列:s:i称为前缀，即prefix(sj),而子序列i：s:©称为后缀suffix(sJs|.i+r)o有两种特殊的情况，即冃或i=j-K

•i：s：i表示空序列表示S中的第j个字符,简记为Sj一般认为子序列与讣算机算法中子串的槪念相当，但是严格地讲，子序列与子串的概念是有区别的0子串是子序列，而子序列不一左是子串。可以通过选取S中的某些字符(或删除S中的某些字符)而形成s的子序列，例如TTT是ATATAT的子序列。而s的子串则是由s中相继的字符所组成，例如TAC是AGTACA的子串，但不是TTGAC的子串。如果t是s的子串，则称s是t的超串。两个序列s和(的连接用s++t来表示，如：ACC++CTA=ACCCTA字符串操作除连接操作之外，另有一个k操作，即删除一个字符串两端的字符。其左义如下:prefix(s,l)=sklsUsuffix(sj)=kkklsi：S：j=时殊円序列比较可以分为四种基本情况，具体任务和应用说明如下：假设有两条长度相近的来自同一个字母表的序列，它们之间非常相似，仅仅是有一些细微的差别，例如字符的插入、字符的删除和字符替换，要求找岀这两条序列的差别。这种操作实际应用比较多，例如，有两个实验室同时测宦某个基因的DNA序列，其结果可能不一样，需要通过序列比较来比较实验结果；(2)(3)(4)假设有两条序列.取出前缀和后缀。假设有两条序列,特左的序列模式;假设有两条序列,序列。要求判断是否有一条序列的前缀与另一条序列的后缀相似，如果是，则分别该操作用于大规模DNA(2)(3)(4)假设有两条序列.取出前缀和后缀。假设有两条序列,特左的序列模式;假设有两条序列,序列。要求判断是否有一条序列的前缀与另一条序列的后缀相似，如果是，则分别该操作用于大规模DNA测序中序列片段的组装：要求判断其中的一条序列是否是另一条序列的子序列。这种操作常用于搜索要求判断这两条序列中是否有非常相似的子序列。这种操作可用于分析保守当然，进行序列比较时，往往还需要说明是采取全局比较，还是采取局部比较。全局比较是比较两个完整的序列，而局部比较是找出最大相似的子序列。2、编辑距离(EditDistance)观察这样两个DNA序列：GCATGACGAATCAG和TATGACAAACAGC*—眼看上去，这两个序列并没有什么相似之处，然而如果将第二条序列错移一位，并对比排列起来以后，就可以发现它们的相似性。GCATGACGAATCAGHillIITATGACAAACAGC如果进一步在第2条序列中加上一条横线，就会发现原来这两条序列有很多相似之处。GCATGACGAATCAGHillIIIIITATGAC-AAACAGC上而是两条序列相似性的一种左性表示方法，为了说明两条序列的相似程度，还需要左疑计算(Sellers1974：SmithandWaterman1981)。有两种方法可用于虽:化两个序列的相似程度，一为相似度，它是两个序列的函数，其值越大，表示两个序列越相似。与相似度对应的另一个概念是两个序列之间的距离。距离越大，则两个序列的相似度就越小。在大多数情况下，相似度和距离可以交互使用，并且距离越大,相似度越小，反之亦然。但一般而言，相似度使用得较多，并且灵活多变。最简单的距离就是海明(Hamming)距离。对于两个长度相等的序列，海明距离等于对应位字符不同的个数。例如，图3.1是3组序列的海明距藹。s= AATAGCAAAGCACACAt= TAAACATAACACACTAHammingDistance(s,t)=2 3 6图3.1海明距藹使用距离来计算不够灵活，这是因为序列可能具有不同的长度，两个序列中各位置上的字符并不一左是真正的对应关系。例如，在DNA复制的过程中，可能会发生象删除或插入一个碱基的错误，虽然两个序列的其他部分相同，但由于位置的移动导致海明距离的失真。就上述例子最右边的情况，海明距禽为6,简单地从海明距离来看，两个序列差别很大(整个序列的长度只有8bp),但如果从s中删除G,从t中删除T,则两个序列都成为ACACACA,这说明两个序列仅仅相差两个字符。实际上，在许多情况下，直接运用海明距离来衡量两个序列的相似程度是不合理的。为了解决字符插入和删除问题，引入字符编辑操作(EditOpcndion)的概念，通过编辑操作将一个序列转化为一个新序列。用一个新的字符“/代表空缺(Space),并泄义下述字符编辑操作：Match(a»a)Delete(a,・)Replace(a,b)Insert b)—表示字符匹配-从第一个序列删除一个字符，或在第二个序列相应的位置插入空白字符-以第二序列中的字符b替换第一个序列中的字符a,a=b—在第一个序列插入空白字符，或删除第一个序列中的对应字符很显然，在比较两个序列s和t时，在s中的一个删除操作等价于在(中对应位置上的一个插入操作，反之亦然。需要注意的是两个空缺字符不能匹配，因为这样的操作没有意义。引入上述编辑操作后，重新计算两个序列的距离，就成为编借距离。以上的操作仅仅是关于序列的常用操作，在实际应用中还可以引入复杂的序列操作。下而是两条序列的一种比对：ACCGACAATATGCATAIlliIATAGGTATAACAGTCA如果将第2条序列头尾倒置，可以发现两条序列惊人的相似：ACCGACAATATGCATAIIIIIIIIIIIIIIACTGACAATATGGATA下而两条序列有什么关系？如果将其中一条序列中的碱基替换为其互补碱基，就会发现其中的关系。CTAGTCGAGGCAATCTGAACAGCTTCGTTAGT3、通过点矩阵分析两条序列的相似之处进行序列比较的一个简单的方法是“矩阵作图法”或“对角线作图”，这种方法是由Gibb首先提出(Gibb“ndMcIntyre,1970)。将两条待比较的序列分别放在矩阵的两个轴上，一条在X轴，从左到右，一条在Y轴，从下往上，如图3.2所示。当对应的行与列的序列字符匹配时，则在矩阵对应的位置作出“点”标记。

ATCGATCGIlliATCG图3.2序列比较矩阵标记图显然，如果两个序列完全相同，则在矩阵主对角线的位置都有标记。如果两个序列存在相同的子序列，则对于每一个相同的子序列对，有一条与对角线平行的由标记点所组成的斜线，如图3.3所示。而对于两条互为反向的序列，则在反对角线方向上有标记点组成的斜线，如图3.4所示。图3.3相同子序列矩阵标记图图3.4反向序列矩阵标记图对于序列标记图中非重叠的与对角线平行斜线，可以组合起来，形成与两个序列的一种比对。在两个子序列的中间插入符号，表示插入空白字符。在这种对比之下分析两条序列的相似性，如图3.5所示。找两条序列的最佳比对（对应位置等同序列字符最多），实际上就是在序列标记图中找非重叠平行斜线最长的组合。TCCGTCCGTATTAT——GCCTIIIIlliTATCCGCCTC图3.5多个相同子序列矩阵标记图4、序列的两两比对序列的两两比对(PairwiseSequenceAlignment)就是按字符位置重组两个序列，使得两个序列达到一样的长度。设两个序列分别是s和(，在s或t中插入空缺符号，使s和t达到一样的长度❾图3.6是对序列AGCACACA和ACACACTA的两种比对结果。AG-CACACAACACACT-As: AG-CACACAACACACT-At: A-CACACTAMatch(A,A)Delete(G,・Match(A,A)Delete(G,・)Match(C,C)Match(A,A)Match(C,C)Match(A.A)Match(C,C)Insert(T)Match(A.A)Match(A.A)Replace(G,C)Insert(A)Match(C,C)Match(A.A)Match(C,C)Rcplace(A,T)Delete(C,-)Match(ArA)图3.6序列AGCACACA和ACACACTA的两种比对结果下而就不同类型的编辑操作左义函数w,它表示“代价(cost)”或“权重(weight)”o对字母表A中的任意字符也b,定义：w(a,a)=0w(a,b)=1a*b (3-1)w(a,・)=w(・，b)=1这是一种简单的代价泄义，在实际应用中还需使用更复杂的代价模型。一方面，可以改变各编辑操作的代价值。例如，在蛋白质序列比较时，用理化性质相近的姐基酸进行替换的代价应该比完全不同的氨基酸替换代价小。另一方而，也可以使用得分(score)函数来评价编借操作(Altschul,1993)。这里给出一种基本的得分函数：p(a,a)=1p(a,b)=-1a*b (3-2)p(a「)=w(b)=-2卜•而给出在进行序列比对时常用的概念:

（1） 两个序列S和t的比对代价等于将S转化为I所用的所有编辑操作的代价和:（2） s和t的最优比对是所有可能的比对中代价最小的一个比对：（3） s和t的真实距离应该是在代价函数w值最优时的距离，记为dw（s.t）o使用前面w的左义，可以得到下列代价：s:AGCACAC-At:A-CACACTAcost=2进行序列比对的目的是寻找一个代价最小的比对。5.用于序列相似性打分的权值矩阵（Wei曲Matrices）无论是3-1式还是3・2式，都是简单相似性评价模型，在计算距离或打分时对字符替换或插入、删除进行统一的处理，没有考虑“同类字符”替换与“非同类字符”替换的差别。实际上，不同类型的字符替换，其代价或得分是不一样的，特別是对于蛋白质序列。某些氨基酸可以很容易地相互取代而不用改变它们的理化性质，例如体积小并且是疏水残基的异亮氨酸和颉氨酸互换，极性残基线姐酸和苏氨酸的互换。理化性质相近的氨基酸残基之间替换的代价显然应该比理化性质相差甚远的氨基酸残基替换代价小，或者得分髙。同样，保守的氨基酸替换打分高于非保守的氨基酸替换。这样的打分方法在比对非常相近的序列以及差异极大的序列得出不同的分值。这就是提出打分矩阵（或者称为取代矩阵）的原由。在打分矩阵中.详细地列出齐种字符替换的得分，从而使得计算序列之间的相似度更为合理。在比较蛋白质时，我们可以用打分矩阵来增强序列比对的敏感性。打分矩阵是序列比较的基础，选择不同的打分矩阵将得到不同的比较结果，而了解打分矩阵的理论依据将有助于在实际应用选择合适的打分矩阵。以下介绍一些常用的打分矩阵或代价矩阵。（1）核酸打分矩阵设DNA序列所用的字母表为A={A,C,G,T}o等价矩阵（表3.1）ATcGA1ATcGA1000T0100C0010G0001表3.1等价矩阵表表3.2BLAST矩阵ATCGA5-4-4-4T-45-4-4C-4-45-4G-4-4-45ATCGA0551T5015C5105G1550表3.3转移矩阵BLAST矩阵（表3.2）BLAST是目前最流行的核酸序列比较程序，表3.2是其打分矩阵。转移矩阵（Transition/Transversion）（表3.3）核酸的碱基按照环结构分为两类，一类是噪吟（腺噪吟A,鸟瞟吟G），它们有两个环，另一类是嚅咙（胞喀咙C,胸腺«T）,它们的碱基只有一个坏。如果DNA碱基的变化（碱基替换）保持环数不变，则称为转换（transition）,如AtG、CtT。如果环数发生变化，则称为颠换（transversion）,如atC、AtT等。转换发生的频率远比颠换髙，而表3・3所示的矩阵正是反应了这种情况。注意，这里矩阵元素的值与得分相反.代表字符替换代价。（2）蛋白质打分矩阵设蛋白质的字母表为表2.1。（i）等价矩阵(3-3)其中陶代表打分矩阵元素，i、j分别代表字母表第i和第j个字符。氨基酸突变代价矩阵GCM(表3.4)GCM(GeneticCodeMatrix,HaigandHurst,1991)矩阵通过计算一个氨基酸残基到转变到另一个氨基酸残基所需的密码子变化数目而得到，矩阵元素的值对应于代价。如果变化一个碱基使某些氨基酸的密码子改变为另一些氨基酸的密码子，其替换代价为1：如果需要2个碱基的改变，则替换代价为2；以此类推。注意，Met到Tyr的转变是仅有的密码子三个位苣都发生变化的转换。其中Glx代表Gly.Gin或Glu,而Asx则代表Asn或Asp,X代表任意氨基酸。GCM常用于进化距离的计算，其优点是计算结果可以直接用于绘制进化树,但是它在蛋白质序列对准尤其是类似程度很低的序列对准中很少被使用。表3.4氨基酸突变代价矩阵ASGLKVTPEDNIQRFYCHMWZBXAla=A01122111112222222222222Ser=S10112211221121111221222Gly=G11022122112221221221222Leu=L21202121222111122111222Lys=K22220212121111222212122VaI=V12112022112122122212222Thr=T1I221201221121222212222Pro=P112122102222I1222122222Glu-E12121122012212222222122Asp=D12122122101222212122212Asn=N21221212210122212122212Ile=I21211112221021122212222Gln=Q22211221122201222122122Arg=R2I111211222110221111222Phe=F21212122222122011222222Tyr=Y21222222211222101132212Cys=C21122222222221110221222His=H22212221211211212022212Met=M22211112222121232202222Trp=W21112222222221221220222Glx=Z22221222122212222222122Asx=B22222222211222212122212???=X22222222222222222222222疏水矩阵(表3.5)该矩阵是根据氨基酸残基替换前后疏水性的变化而得到得分矩阵。若一次氨基酸替换疏水特性不发生太大的变化，则这种替换得分髙，否则替换得分低。表3.5蛋白质疏水矩阵RKI)EBZSNQGXTHACMPVL1YFWArg=R1()10998866655555433333210Lys=K1()10998866655555433333210Asp=D991()108876665555544433321Glu=E991()108876665555544433321Asx=B88881()1088887777666555443G)x=Z88881()1088881777666555443Scr=S66778810101()10999988777664Asn=N666688101010109999888177664G)n=Q6666881()1010109999888777664Gly=G556688101010109999888877665?T?=X555577999910101()10998888775ThuT555577999910101010998888775His=H55557799991()101010999888775Ala=/\55557799991()101010999888775Cys=C445566888899991()109999885Mcl=M334466888899991010101099887Pro=P33446678888S99910101099987Val=V3344557778888891010101010987Lcu=L33335577778888999101010998llc=l33335577778888999101()10998Tyr=Y?Av23344666677778899991()108Phc=F112244666677778888991()109l*rp=\V001133444555556777888910PAM矩阵PAM矩阵是第一个广泛使用的最优矩阵，它建立在进化的点突变模型PAM(PointAcceptedMutation.Dayhoff1969；Dayhoffetal.J978)基础上。Dayhoff和她的同事们研究了71个相关蛋白质家族的1572个突变，他们发现家族中氨基酸的替换并不是随机的，由此断言一些氨基酸的替换比其它替换更容易发生.英主要原因是这些替换不会对蛋白质的结构和功能产生太大的影响。如果氨基酸的替换是随机的，那么每一种可能的取代频率仅仅取决于不同氨基酸的出现的背景频率。然而在相关蛋白中，存在取代频率大大地倾向于那些不影响蛋白质功能的取代，换句话说，这些点突变已经被进化所接受。这意味着在进化历程上相关的蛋白质在某些位置上可以岀现不同的氨基酸。一个PAM就是一个进化的变异单位，即1%的氨基酸改变，这并不意味着经过100次PAM后，每个氨基酸都发生变化，因为其中一些位置可能会经过多次改变，甚至可能变回到原先的氨基酸，因此另外一些氨基酸可能不发生改变。PAM有一系列的替换矩阵，每个矩阵用于比较具有特泄进化距禽的两个序列。例如，PAM120矩阵用于比较相距120个PAM单位的序列。一个PAM-N矩阵元素(i,j)的值反应两个相距N个PAM单位的序列中第i种氨基酸替换第j种氨基酸的频率。从理论上讲，PAM0是一个单位矩阵，主对角线上的元素值为1,其它矩阵元素的值为0。其它PAM-N矩阵可以通过统汁计算而得到。首先针对那些确信是相距一个PAM单位的序列进行统汁分析，得到PAM1矩阵。PAM1对角线上的元素值接近于1，而其它矩阵元素值接近于0。将PAM1自乘N次，可以得到PAM-NODayhoff等第一次使用了log-odd处理，在这种处理中，矩阵中的取代分值同目标频率与背景频率的比值的自然对数成比例。为了评估目标频率，人们用进化历程上相关的序列来收集对应于一个PAM的突变频率，然后将数据外推至250个PAMo虽然Dayhoff等人只发表了PAM250,但潜在的突变数据可以外推至其它PAM值，产生一组矩阵，在比较差异极大的序列时，通常在较髙的PAM值处得到最佳结果，比如在PAM200到250之间，较低值的PAM矩阵一般使用于髙度相似的序列(AltschulJ991)<>BLOSUM矩阵BL0SUM矩阵是由Henikoff首先提出的另一种氨基酸替换矩阵(Henikoff,1992),采用与PAM冋样的方式可以建立BLOSUM替换矩阵。PAM矩阵是从蛋白质序列的全局比对结果推导岀来的，而BLOSUM矩阵则是从蛋白质序列块(短序列)比对而推导岀来的。但在评估氨基酸替换频率时，应用了不同的策略.基本数据来源于BLOCKS数据库，其中包括了局部多重比对(包含较远的相关序列，同在PAM中使用较近的相关序列相反)。虽然在这种情况下，没有进化模型，但它的优点在于可以通过直接观察获得数据而不是通过外推获得。同PAM模型一样，也有许多编号的BLOSUM矩阵，这里的编号指的是序列可能相同的最髙水平，并且同模型保持独立性(李衍达，等，2000)0

第二节两两比对算法进行序列两两比对最直接方法就是生成两个序列所有可能的比对，分别计算代价函数，然后挑选一个代价最小的比对作为最终结果。但是，两个序列可能的比对数非常多，是序列长度的指数函数。从算法时间复杂性的角度来看，这种比较方法显然不合适。用上一节中所介绍的矩阵分析方法，在寻找斜线及斜线组合时，仍然需要较大的运算量。因此，必须设计髙效的算法以找出最优的比对。著名的Needleman-Wunsch算法就是针对寻求最佳序列比对这一问题所设计的动态规划寻优策略（NeedlemanandWunsch,1970）o下而首先介绍基本的序列比较算法，即动态规划算法（DynamicProgramming;GalilabdGiancarlo,1989;Giegerich,2000;Huang,1993）o动态规划是一种常用的规划方法，往往用于在一个复杂的空间中寻找一条最优路径。对于一个具体的问题，如果该问题可以被抽象为一个对应的图论问题，并且问题的解对应于图中从起点到终点的最短距禽，那么就可以通过动态规划算法解决这个问题。在运用动态规划时，有以下几个要求：（1）首先,搜索问题能够划分成一系列相继的阶段：（2）起始阶段包含基本子问题的解：（3）在后续阶段中，能够按递归方式根据前面阶段的部分结果计算每个部分解：（4）最后阶段包含全局解“K序列两两比对基本算法设序列s、（的长度分别为m和叽考虑两个前缀，o：s：i和o：t：j,i,jni°假如已知序列o：s：i和所有较短子列的最优比对，即已知：（1） O：s：（i.|）和 的最优比对（2） o：s:“・n和o-t'j的最优比对（3） o：s：i和的最优比对则O：S：j和o：t：j的最优比对一泄是上述三种情况之一的扩展，（1） 替换（Si,tj）或匹配（Si,tj）,这取决于Si是否等于tj:（2） 删除（s“ -）:（3） 插入（tj）o可根拯下列递归算式讣算最小值：Co：S：（f0:f：（丿_1J+W（耳心）(3-4)«0:S：门0：/：j）=nin血（（）：s:G_lpo:r：y）+叽»—）心（o：s：e：/：（iJ+w（-心）(3-4)其初值为：血（o：s：o,o：『：o）=O血（o：$:m:r：o）=＜（o:s：（r_1）w汀：。）+w（»—） fori=l,2‘・・・・・・，m （3-5）〃w（o：s：o,o：『：_/）=d*（o：s:o»o•，•（；-!））■*■w（—，°） forj=l,2 n按照这种方法，对于给泄的代价函数w,两个序列所有前缀的编辑距离泄义了一个（m+l）x（n+l）的距离矩阵D=（dij）其中di.j=dw（o：s：i,o：t：j）。对于一个长度为n的序列，有n+1个前缀（包括一个空序列），所以距离矩阵的大小为（m+l）x（n+l）odj.j的计算公式如下：闩+w(s“)闩+w(s“)(3-6)(3-6)dij最小值的三种选择决左了各矩阵元素之间的关系，用图3.7表示:图3.7距离矩阵元素d-j的计算矩阵右卜角元素即为期望的结果：dmn=dw(0：S：m,0：t：n)=dw(SJ)o计算过程从do』开始，可以是按行计算，每行从左到右，也可以是按列计算，每列从上到下。当然，任何计算过程，只要满足在计算d-j时di-JJ.di-JJ-H和di.jj都已经被计算这个条件即可。在讣算di.j后，需要保存di・j是从dx.j、dij.j」、或di・jj中的哪一个推进的，或保存计算的路径，以便于后续处理。上述计算过程到dm.n结束。与计算过程相反，求最优路径或最优比对时，从din开始，反向前推。假设在反推时到达di.j,现在要根据保存的计算路径判断di.j究竟是根据di.u>dMj.H和中的那一个汁算而得到的。找到这个点以后，再从此点岀发，一直到du・o为止。走过的这条路径就是最优路径(即代价最小路径),其对应于两个序列的最优比对。设s=AGCACACA・t=ACACACTA,则距离矩阵和最优路径如图3.8所示。在图3.8中，画出了一条穿越距禽矩阵的路径，该路径表明了如何通过合理的比对达到最小的代价。其中斜线表示匹配或替换，垂直线表示删除，而水平线则表示插入。由该路径可以得到比对s: AGCACAC-AIIIIIIIt: A-CACACTA从图3-8可知dw(s,t)=20值得注意的一点是，在有些情况下，最优比对并非唯一，亦即存在几条最优路径。以上计算是在代价函数的基础上进行的，在实际应用中也常常用“得分(scores)”表示相似程度,得分越髙，序列就越相似，因为计算方式刚好与代价函数相反，所以要求出最大得分的路径。一般来讲,距禽总是正的，代价函数亦如此，而得分则可正可负，因此用起来就更加灵活，例如后而要介绍的局部相似性计算方法就使用了得分的概念。在计算序列相似程度时还应该考虑序列长度的影响。令Cw(s,t)表示两个长度分別为m和n的序列的相似性得分，设cw(s,t)=99,即两个序列有99个字符一致。如果m=n=100,则可以说这两个序列非常相似，几乎一样。但如果m=n=1000,则仅有10%相同。所以在实际序列比较时，使用相对长度的得分就更有意义,即用sim(s,t)=2*cw(sj)/(m+n)而非cw(s,t)作为衡量序列相似性的指标。

t—>ACACACTA012345678A1•♦••O 1 2 3 4 5 6 7G2 1...1 2 3 4 5 6 7SC3 2工2 2 3 4 5 6JA4 3 2'X.2 2 3 4 5c5 4 3 2耳•・2 2 3 4A6 5 4 3 2乜..2 3 3c7 6 5 4 3 2*1…•…2. 3A87654322X-2图3.8距离矩阵和最优路径从数据结构di.j本身及其讣算过程来看，序列两两比对基本算法的空间复杂度和时间复杂度都是0(mn).如果两个序列的长度相近，空间和时间复杂度为0(2)。2>子序列与完整序列的比对序列比较的基本动态规划算法找出两条序列的最佳比对.而不在乎它是否具有生物学意义。有些同源序列虽然全序列的相似性很小，但是存在髙度相似的局部区域。如果在进行序列比较时注重序列的局部相似性，则可能会发现重要的比对，因此，不能仅仅只关注最佳的一个。Smith和Waterman在Needleman-Wunsch算法的基础上进行改进，提岀的局部比对算法(SmithandWaterman,1981)。后来其他人又进一步改进，形成改良Smith-Waterman算法(AltschulandErickson」986;WatermanandEggert.1987),算法将寻找K种最好的但不相互交叉的比对方式作为目标，这些思想后来都在SIM算法(Huangetal.J990)的发展中得以体现。FASTA程序包中的LALIGN程序提供了有用的SIM工具(PearsonJ996)o对于比对多模块的蛋白质而言，寻找次优比对尤为重要。LALIGN程序可以被用来得到多个最好的局部比对，而一个标准的Smith-waterman算法只会报告出最好的一个比对，改良的算法会报告出第二和第三的比对方式，从而显示岀功能结构域。下而分别介绍各种局部比对方法。在有些情况下，我们需要将一个较短的序列(或探测序列，或模式序列)与一个较长的完整序列比较，试图找出局部的最优匹配。这是一种局部的比对(localalignment),其左义如下：给左两个序列U：s：m和o：t：n,从I中寻找一个子列i：t：j使得dw(S,i：t：j)最小，0<i<j<no已有许多高效的算法可解决局部比对问题，而动态规划算法也只要作一点小小的修改就可以用于局部比对。如果a,b相似如果如果a,b相似如果a,b不相似(3-7)w(a,b)>0w(a,b)<0w(a,-)<0w(-.b)<0以“(T作为子序列相似与否的阈值。计算的目标是寻求最大的相似度(或得分)，而不是最小距离,因而公式(3-6)中的min改变为max。局部比对意味着不计删除前缀o：t：i的得分，也不计删除后缀的得分，即(3-8)

(3-(3-#)d：d：:=max<1•J+w(-心)阈值“0”意味着矩阵中的“0”元素分布区域对应于不相似的子序列，而正数区域则是局部相似的区域。最后，在矩阵中找最大值，该值就是最优的局部比对得分，该值对应的点为序列局部比对的末点，然后反向推演前而的最优路径，直到局部比对的起点。4、半全局比较所谓半全局比较就是在评价序列比对时不计终端“空缺”(endspace,或空位)的得分或代价(SetubalandMeidanis,1997)。在下而的讨论中，以“空位”代表字符，表示插入或删除操作。我们仍然采用基本的比对算法，而改进对两端空位的打分函数。终端空位是那些出现在序列第一个字符前或最后一个字符之后的空位。例如，图3.10(a)第二个序列中所有的空位都是终端空位。这两个序列长度相差较大，一个序列长度8,而另一个序列的长度18.如果两个序列的长度相差比较大，那么，在最终的比对中一淀存在很多空位，从而引起得分下降或代价增大。然而，如果忽略这些终端空位，则比对相当好，6个匹配，1个失配，1个空位。但是根据基本的序列比对算法，图3.10(a)的比对不是最优的，最优的比对似乎应该是图3.10(b)所示的比对，因为其得分更高(同样多的空位，但除了空位以外全是匹配)。尽管具有更髙的得分，第二个序列的字符全部匹配，但从寻找相似区域的观点来看，图3.10(b)所示的并不是我们所感兴趣的比对。为了匹配更多的字符，第二个序列被严重地割裂开。如果我们在一段长序列寻找一段与另一个短序列区域相似的区威，无疑图3」0(a)所示的排列更好。CAGCA-CTTGGATTCTCGG——CAGCGTGG (a)CAGCACTTGGATTCTCGGCAGC G-T GG(b)图3・10终端空白序列假设有两个序列S和toS后而的空位与t的后缀匹配，去掉这部分比对，仅保留S与t前缀的比对,被去掉的比对不再计分。因此，为了得到不计S右端空位得分的最优比对，我们需要寻找S和t前缀的最优比对。由于在基本的序列比对算法中，矩阵元素d订代表o：s：j和o：t：j的相似性，因而可以取最后一行的最大值，即：n=max； (3-14)取最后一行的最大值作为比对的得分，就相当于不计序列S右端的空位得分(负分)。注意，这里矩阵元素尙代表o：s：i和u：t：j前缀的最大相似性，并且计算公式就是(3-6)式(min改为max)。按照上述方法，我们可以同样取最后一列的最大值作为比对的得分，从而达到不计序列I后而空位的目的。进一步将两者组合起来，形成一种不计两个序列末尾空位的最优比对方法。现在我们再讨论序列的前端空白。情况与末端空白相似，只要在设置初始条件时将矩阵第一行和第一列元素的值为“0”,即达到不计序列前端空位得分的目的。综合上述分析结果，对两个序列进行半全局比较的算法如下：将矩阵的第一行和第一列元素置为“0”，按行或按列优先的次序计算矩阵D所有元素di・j(0<i<m,0<j<n),取最后一行和最后一列的最大值，而在搜索最优路径时，从矩阵最后一行和最后一列的最大值岀发，按照基本算法中的方法返回。半全局比较算法与基本算法在计算dM时的区别归纳为下列四个方而：第一行初始值为“0”，表示不计第一个序列的前端空位；寻找最后一行的最大值，表示不计第一个序列的末端空位：第一列初始值为“0”，表示不计第二个序列的前端空位；寻找最后一列的最大值，表示不计第二个序列的末端空位。5、关于连续空位的问题左义：“k阶间隙(Gap)"是一个连续的“空位”区域，其空缺字符的数目k>l(SetubalandMeidanis,1997)o在研究序列突变时，k阶间隙出现的可能性要大于k个不连续的空位，这是因为一个间隙对应于一串字符的插入或者删除，对应于一个突变事件，不连续的空位可能对应于多个不同的突变事件，而一个突变的发生比多个突变同时发生的可能性要大。从这个意义上说，我们希望将k连续空位与k个孤立空位在打分方而区別开来，k连续空位的得分(或代价)应该比k个孤立空位的得分高(或代价低)。妇外，在将cDNA与基因DNA比较时，通过改进的空位打分方式可以正确处理内含子。真核生物的基因是非连续的，在编码区域的外显子之间插入了内含子。在基因转录过程中，内含子被剔除，外显子连接起来形成完整的编码区域，并表现为mRNAo通过反转录，根据mRNA形成cDNA“因此，同一个基因的cDNA序列与DNA序列的差别主要反映在内含子的删除。这样，在根据基因的cDNA寻找DNA时，我们也希望对连续空位的处理与孤立空位的处理有所区别*然而直到现在，我们在算法中还没有区分单个空位及连续的多个空位，这意味着对于间隙的代价是一个线性的函数。设w(k)代表空位得分函数，其中k是连续空位的个数，则：w(k)=-bk (3-15)这里b(>0)是单个“空位”得分的绝对值。在本节中，将讨论解决序列比对连续空位问题的算法，该算法的时间复杂性是O(2)。该算法与基本算法类似.主要差别在于得分函数不具有加和性，即不能将一个比对分成两个部分.而希望总的得分是两个部分的和。然而，当比对的分割点跨越“块(Block)”边界时，得分的加和性依然成立•任何一个比对可以被唯一地分为若干个相继的块。有三类块：两个字符的比对：与序列s空位进行比对的t的最大连续字符序列：与序列t空位进行比对的s的最大连续字符序列。“最大”意味着不能再扩展。图3.11显示了一个比对及块的分割。上述第一类块的得分是p(a.b),a>b是两个比对的字符；其它两类的得分是-w(k),k为连续的空位数。sAAC——AATTCCGACTACtACTACCT CGC--saIaIcI——IaIattccgIaIcltIactaIcltIaccItI 1clgIcl—图3.11比对及块分割现在，对于一个比对的打分不再是按照单个比对的列进行，而是按块进行。只有当打破块的边界时,得分的加和性才成立。对基本的动态规划方法进行修改，优化过程也不再是针对每个字符，而是针对每个块。另外，块和块的连接不是任意的，第二类和第三类的块不能和相同类型的块连接匚对于每一个点

对（ij）•我们不仅要保存o：s：i和0心的最优比对的得分，还特别要保存那些终结于特泄块的前缀的最佳得分。为比较序列S（长度为m）和序列t（长度为n）,我们使用三个（m+1）x（n+1）的矩阵，每一个矩阵对应于一种块的类型。矩阵A对应结朿于字符对比的块，a»用于保存°:s：i和。心的最优比对的得分•其所对应的比对的最后一个比对操作是字符匹配或替换。矩阵B对应结束于s空位的块，切保存的也是o：s：i和庶:j的最优比对的得分，但是其所对应的比对的最后一个比对操作是在U：S：i末尾插入空位。矩阵C对应结束于t空位的块，5保存的也是u：s：i和。:冷的最优比对的得分，但是其所对应的比对的最后一个比对操作是在庶:j末尾插入空位。各矩阵第一行和第一列初始值的设左如下：(3-16)(3-17)(3-18)^<).0=°(3-16)(3-17)(3-18)Cj()=所有其它元素初始值为（这对应于取最大值的操作）o递归计算过程如下：(3-19)%=吨J)+max］侏|jtfiTJT(3-19)bt；=max<；S+'KIforl<k<jci.i-k+W)for1<k<j(3-20)G—j+'V(k)forl<k<iCii=max<b*+w伙)forl<k<i(3-21)注意，由于块可能有不同的长度，所以矩阵B和矩阵C中的元素依赖于多个以前的值。另外计算矩阵B时，不需要用到B以前的值，这因为第二种类型的块不可能连续出现。对于C也一样。至于最终的结果Sim（S,t）,取Hnunxbm.n和Cm.n的最大值q在序列长度相近的条件下，上述算法的时间复杂度为0（2）。比起标准算法，其多花的时间主要用于处理连续的空位（间隙几那么是否可以改进间隙的代价函数，而使得算法的时间复杂度降低为0（2）呢？下而讨论这个问题。如果认为k个连续空位比k个孤立空位岀现的可能更大，则w(k)>kw(l) (3-22)或更一般地，w(ki+k：+..・+kn)Aw(ki)+w(k2)+...+w(kn) (3-23)可以用下式重新计算连续“空位”的得分：w(0)=0 (3-24)w(k)=-h-g(k-l), k>l (3-25)1】、g为参数,"・h”代表第一个“空位”的代价，而则代表后续“空位”的代价，h、g>0,h>ge

选用这样的函数作为间隙的打分函数，其含义是：一段连续空位的第一个代价较大，而其它后续空位的代价较