生物信息学讲义-双序列比对

第二章双序列比对

PairwiseSequenceAlignment

第一节引言同源（homology）-具有共同的祖先垂直同源（ortholog）水平同源（paralog）相似（similarity）

同源序列一般是相似的，相似序列不一定是同源的

通过点矩阵进行序列比较

编辑距离（editdistance)相似性得分第二节替换记分矩阵（1）核酸打分矩阵设DNA序列所用的字母表为

={A，C，G，T}a.等价矩阵(unitarymatrix)b.BLAST矩阵c.转换-颠换矩阵（transition-transversionmatrix）（嘌呤：腺嘌呤A，鸟嘌呤G；嘧啶：胞嘧啶C，胸腺嘧啶T）ATCGA1000T0100C0010G0001ATCGA5-4-4-4T-45-4-4C-4-45-4G-4-4-45ATCGA1-5-5-1T-51-1-5C-5-11-5G-1-5-51表3.1等价矩阵表表3.3转移矩阵表3.2BLAST矩阵（2）蛋白质打分矩阵（i）等价矩阵（ii）遗传密码矩阵（geneticcodematrix,GCM)（iii）疏水性矩阵

（hydrophobicmatrix)（iv）PAM矩阵（pointacceptedmatrix,PAM）（v）BLOSUM矩阵（BLOckSUbstitutionMatrix，BLOSUM）其中Rij代表打分矩阵元素i、j分别代表字母表第i和第j个字符。遗传密码矩阵遗传密码矩阵通过计算一个氨基酸变成另一个氨基酸所需的密码子变化的数目而得到。通常为1或2，只有Met到Tyr为3。ASGLKVTPEDNIQRFYCHMWZBXA01122111112222222222222S10112211221121111221222G11022122112221221221222L21202121222111122111222K22220212121111222212122V12112022112122122212222T11221201221121222212222P11212210222211222122222E12121122012212222222122D12122122101222212122212N21221212210122212122212I21211112221021122212222Q22211221122201222122122R21111211222110221111222F21212122222122011222222Y21222222211222101132212C21122222222221110221222H22212221211211212022212M22211112222121232202222W21112222222221221220222Z22221222122212222222122B22222222211222212122212X22222222222222222222222GCM矩阵疏水矩阵RKDEBZSNQGXTHACMPVLIYFWR1010998866655555433333210K1010998866655555433333210D9910108876665555544433321E9910108876665555544433321B8888101088887777666555443Z8888101088887777666555443S667788101010109999887777664N666688101010109999888777664Q666688101010109999888777664G556688101010109999888877665X555577999910101010998888775T555577999910101010998888775H555577999910101010999888775A555577999910101010999888775C4455668888999910109999885M334466888899991010101099887P33446678888899910101099987V3344557778888891010101010987L33335577778888999101010998I33335577778888999101010998Y2233446666777788999910108F1122446666777788889910109W001133444555556777888910PAM&BLOSOM

这类矩阵列出同源蛋白质在进化过程中氨基酸变化的可能性。这类矩阵是基于进化原理的证据：编码相同蛋白质的基因随着进化发生分歧，相似度降低。科学用得多

PAM矩阵（pointacceptedmutaion）

基于氨基酸进化的点突变模型

如果两种氨基酸替换频繁，说明自然界接受这种替换，那么这对氨基酸替换得分就高一个PAM就是一个进化的变异单位,即1%的氨基酸改变

但这并不意味100次PAM后，每个氨基酸都发生变化，因为其中一些位置可能会经过多次突变，甚至可能会变回到原来的氨基酸。PAM矩阵的制作步骤构建序列相似（大于85％）的比对计算氨基酸j的相对突变率mj（j被其他氨基酸替换的次数）针对每个氨基酸对i和j,计算j

被i

替换次数替换次数除以相对突变率（mj）利用每个氨基酸出现的频度对j

进行标准化取常用对数，得到PAM-1(i,j)将PAM-1自乘N次，可以得到PAM-nPAMMatrices

MutationsacceptedbynaturalselectionConstructingPAMMatrix:TrainingDataPAM:PhylogeneticTreePAM:AcceptedPointMutationMutabilityofResiduejTotalMutationRateisthetotalmutationrateofallaminoacidsNormalizeTotalMutationRateto1%Thisdefinesanevolutionaryperiod:theperiodduringwhichthe1%ofallsequencesaremutated(acceptedofcourse)MutationProbabilityMatrixNormalizedSuchthattheTotalMutationRateis1%MutationProbabilityMatrix(transposed)M*10000elementsareshownmultipliedby10,000From:http://www.icp.ucl.ac.be/~opperd/private/pam1.htmlPAM-250PAM60—60%,PAM80—50%,PAM120—40%PAM-250matrixprovidesabetterscoringalignmentthanlower-numberedPAMmatricesforproteinsof14-27%similarityPAMMatrix:Assumptions

PAM=%AcceptedMutations:

1500changesin71groups>85%similarity

BLOSUM=BlocksSubstitutionMatrix:

2000“blocks”from500families

TwoclassesofwidelyusedproteinscoringmatricesBLOSUM62ChoiceofScoringMatrix针对不同的进化距离采用PAM矩阵序列相似度=40%50%60%

|||打分矩阵=PAM120PAM80PAM60PAM250→14%-27%

PAM矩阵与BLOSUM矩阵的比较第三节双序列比对算法序列的两两比对 （PairwiseSequenceAlignment）

按字符位置重组两个序列，使得两个序列接近一样的长度序列两两比对基本算法直接方法——生成两个序列所有可能的比对，分别计算代价函数，然后挑选一个代价最小的比对作为最终结果，需要计算2300

次——天文数字ATTC………CGAAGA

AGTC………GAAGGT假设比较300个氨基酸长度的两条序列动态规划方法DynamicProgramming起点终点ATTC………CGAAGA

AGTC………GAAGGTATTC………CGAAGAGTC………GAAGGAT+（1）ATTC………CGAAGAAGTC………GAAGG-T+（2）ATTC………CGAAGAGTC………GAAGGTA-+（3）最短路经问题起点终点C1C2W1

W2路径1：C1+w1？路径2：C2+w2？

取最小值！算法求解:

从起点到终点逐层计算

计算过程：

计算过程：按行计算其他方式

计算过程：（3）求最佳路径算法分析： 数据结构di,j

空间复杂度：O(mn)

时间复杂度：O(mn)由于在所考虑的子问题空间中，总共有θ(mn)个不同的子问题，因此，用动态规划算法自底向上地计算最优值能提高算法的效率。矩阵赋值算法

fori=0to

length(A)F(i,0)←0

forj=0to

length(B)F(0,j)←0

fori=1to

length(A)

forj=1to

length(B){Choice1←F(i-1,j-1)+S(A(i),B(j))Choice2←F(i-1,j)+dChoice3←F(i,j-1)+dF(i,j)←max(Choice1,Choice2,Choice3)}

算法程序反向构造匹配序列AlignmentA←""AlignmentB←""i←length(A)j←length(B)

while(i>0andj>0){Score←F(i,j)ScoreDiag←F(i-1,j-1)ScoreUp←F(i,j-1)ScoreLeft←F(i-1,j)

if(Score==ScoreDiag+S(A(i-1),B(j-1))){AlignmentA←A(i-1)+AlignmentAAlignmentB←B(j-1)+AlignmentBi←i-1j←j-1}

else

if(Score==ScoreLeft+d){AlignmentA←A(i-1)+AlignmentAAlignmentB←"-"+AlignmentBi←i-1}

otherwise(Score==ScoreUp+d){AlignmentA←"-"+AlignmentAAlignmentB←B(j-1)+AlignmentBj←j-1}}子序列与完整序列的比对

----AGCT----ATGCAGCTGCTT目标： 使S(s,i:t:j)最大序列S:序列t:ij不计前缀0:t:i的得分,也不计删除后缀的j+1:t:|t|得分不计删除后缀的j+1:t:|t|得分

——处理最后一行+p(-,tj)不计前缀0:t:i的得分——处理第一行tsACACACTA000000000C-101010100A-200212110C-3-11132321A-4-20224444最后一行不计代价子序列s在全序列t的后面出现时不会被罚分影响三、比对的统计学显著性（1）典型方法：将两条待比较的序列分别随机打乱使用相同的程序与打分函数（或打分矩阵）进行比对计算这些随机序列的相似性得分重复这一过程（50～100次）用和分别表示其平均值与标准差。设原来两条序列的比对得分为x，利用下式计算大于或等于x的比对得分概率：z=(x-)/根据z值判断两个序列相似得分的显著性，当z值是3.1、4.3、5.2时，x出现的概率为10-3、10-5、10-7Z>5，同源；Z<3,不同源；Z=3~5,可能同源经验法则（针对蛋白质序列）：①如果两个序列的长度都大于100，在适当地加入空位之后，它们配对的相同率达到25%以上，则两个序列相关；②如果配对的相同率小于15%，则不管两个序列的长度如何，它们都不可能相关；③如果两个序列的相同率在15%25%之间，它们可能是相关的。数据库的搜索简介第四节双序列比对工具

数据库查询为生物学研究提供了一个重要工具，在实际工作中经常使用。然而，在分子生物学研究中，对于新测定的碱基序列或由此翻译得到的氨基酸序列，往往需要通过数据库搜索，找出具有一定相似性的同源序列，以推测该未知序列可能属于哪个基因家族，具有哪些生物学功能。对于氨基酸序列来说，有可能找到已知三维结构的同源蛋白质而推测其可能的空间结构。因此，数据库搜索与数据库查询一样，是生物信息学研究中的一个重要工具。

数据库搜索的基础是序列的相似性比对，即双序列比对(pairwisealignment)。新测定的、希望通过数据库搜索确定其性质或功能的序列称作检测序列(probesequence)；通过数据库搜索得到的和检测序列具有一定相似性的序列称目标序列(subjectsequence)。为了确定检测序列和一个已知基因家族之间的进化关系，在通过数据库搜索得到某些相似序列后，还需要判断其序列相似性程度。如果检测序列和目标序列的相似性程度很低，还必须通过其他方法或实验手段才能确定其是否属于同一基因家族。一、BLAST简介BLAST程序是目前最常用的基于局部相似性的数据库搜索程序，它们都基于查找完全匹配的短小序列片段，并将它们延伸得到较长的相似性匹配。它们的优势在于可以在普通的计算机系统上运行，而不必依赖计算机硬件系统而解决运行速度问题。BLAST数据库搜索策略BLAST仅通过部分而不是全部序列计算最适联配值——赢得搜索速度

比对统计学意义的评价——E值(E-Value)P值(P-Value)(概率值)BLAST程序中使用了E值而非P值，这主要是从直观和便于理解的角度考虑。比如E值等于5和10，总比P值等于0.993和0.99995更直观。但是当E<0.01时，P值与E值接近相同参数K和λ可分别被简单地视为搜索步长(searchspacesize)和计分系统(scoringsystem)的特征数BLAST软件包实际上是综合在一起的一组程序，不仅可用于直接对蛋白质序列数据库和核酸序列数据库进行搜索，而且可以将检测序列翻译成蛋白质或将数据库翻译成蛋白质后再进行搜索，以提高搜索结果的灵敏度。BLAST程序检测序列和数据库类型程序名检测序列数据库类型方法Blastp蛋白质蛋白质用检测序列蛋白质搜索蛋白质序列数据库Blastn核酸核酸用检测序列核酸搜索核酸序列数据库Blastx核酸蛋白质将核酸序列按6条链翻译成蛋白质序列后搜索蛋白质序列数据库Tblastn蛋白质核酸用检测序列蛋白质搜索由核酸序列数据库按6条链翻译成的蛋白质序列数据库Tblastx核酸核酸将核酸序列按6条链翻译成蛋白质序列后搜索由核酸序列数据库按6条链翻译成的蛋白质序列数据库

对一般用户来说，目前常用的办法是通过NCBI、EBI等国际著名生物信息中心的BLAST服务器进行搜索。需要说明的是，各生物信息中心BLAST用户界面有所不同，所提供的数据库也可能不完全相同，使用前最好先进行适当的选择。BLAST应用实例多结构域蛋白(H1N1)

脂质运载蛋白多结构域蛋白(H1N1)的BLAST检索>gi|224983683|pdb|3GBN|BChainB,CrystalStructureOfFabCr6261InComplexWithThe1918H1n1InfluenzaVirusHemagglutininGLFGAIAGFIEGGWTGMIDGWYGYHHQNEQGSGYAADQKSTQNAIDGITNKVNSVIEKMNTQFTAVGKEFNNLERRIENLNKKVDDGFLDIWTYNAELLVLLENERTLDFHDSNVRNLYEKVKSQLKNNAKEIGNGCFEFYHKCDDACMESVRNGTYDYPKYSEESKLNREEIDGVSGR多结构域蛋白(H1N1)的BLAST检索H1N1聚合酶序列BLAST结果综述BLAST结果表述BLAST结果逐条显示BLAST结果逐条显示BLAST：改变打分矩阵的作用脂质运载蛋白序列>sp|P31025|LCN1_HUMANLipocalin-1OS=HomosapiensGN=LCN1PE=1SV=1MKPLLLAVSLGLIAALQAHHLLASDEEIQDVSGTWYLKAMTVDREFPEMNLESVTPMTLTTLEGGNLEAKVTMLISGRCQEVKAVLEKTDEPGKYTADGGKHVAYIIRSHVKDHYIFYCEGELHGKPVRGVKLVGRDPKNNLEALEDFEKAAGARGLSTESILIPRQSETCSPGSD使用Blosum62矩阵搜索使用PAM30矩阵搜索FastA简介FastA算法是由Lipman和Pearson于1985年发表的（Lipman和Pearson，1985）。FastA的基本思路是识别与代查序列相匹配的很短的序列片段，称为k-tuple。

蛋白质序列数据库搜索时，短片段的长度一般是1～2个残基长；DNA序列数据库搜索时，通常采用稍大点的值，最多为6个碱基。通过比较两个序列中的短片段及其相对位置，可以构成一个动态规划矩阵的对角线方向上的一些匹配片段。

FastA程序采用渐进（heuristicapproach）算法将位于同一对角线上相互接近的短片段连接起来。也就是说，通过不匹配的残基将这些匹配残基片段连接起来，以便得到较长的相似性片段。这就意味着，FastA输出结果中允许出现不匹配残基。这和BLAST程序中的成对片段类似。如果匹配区域很多，FastA利用动态规划算法在这些匹配区域间插入空位。

由FastA搜索产生的典型输出结果的第一行列出程序名称和版本号，以及该程序发表的杂志。接下来列出所提交的序列，然后是所用参数和运行时间，紧跟这些一般信息的是数据库搜索结果。首先列出搜索得到的目标序列简单说明，其数目可由用户定义。所列出的目标序列的信息包括：序列所在数据库名称的缩写，目标序列的标识码、序列号和序列名等部分信息。括号中标明匹配部分的残基数。紧接着是由程序计算得到的初始化和优化后的分数值。最后一列是期望值即E值，用来判断比对结果的置信度。接近于0的E值表明两序列的匹配不大可能是由随机因素造成的。以两条氨基酸序列的比较为例介绍算法的基本思路，算法可以分为4步：第一步：FASTA首先找出进行比较的两条序列所有长度为k-tuple的连续的一致序列片段。例如以下两条蛋白质序列：··设k-tupl