生物信息學(xué)講義-雙序列比對

第二章雙序列比對

PairwiseSequenceAlignment

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

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

通過點矩陣進(jìn)行序列比較

編輯距離（editdistance)相似性得分第二節(jié)替換記分矩陣（1）核酸打分矩陣設(shè)DNA序列所用的字母表為

={A，C，G，T}a.等價矩陣(unitarymatrix)b.BLAST矩陣c.轉(zhuǎn)換-顛換矩陣（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轉(zhuǎn)移矩陣表3.2BLAST矩陣（2）蛋白質(zhì)打分矩陣（i）等價矩陣（ii）遺傳密碼矩陣（geneticcodematrix,GCM)（iii）疏水性矩陣

（hydrophobicmatrix)（iv）PAM矩陣（pointacceptedmatrix,PAM）（v）BLOSUM矩陣（BLOckSUbstitutionMatrix，BLOSUM）其中Rij代表打分矩陣元素i、j分別代表字母表第i和第j個字符。遺傳密碼矩陣遺傳密碼矩陣通過計算一個氨基酸變成另一個氨基酸所需的密碼子變化的數(shù)目而得到。通常為1或2，只有Met到Tyr為3。ASGLKVTPEDNIQRFYCHMWZBXA01122111112222222222222S10112211221121111221222G11022122112221221221222L21202121222111122111222K22220212121111222212122V12112022112122122212222T11221201221121222212222P11212210222211222122222E12121122012212222222122D12122122101222212122212N21221212210122212122212I21211112221021122212222Q22211221122201222122122R21111211222110221111222F21212122222122011222222Y21222222211222101132212C21122222222221110221222H22212221211211212022212M22211112222121232202222W21112222222221221220222Z22221222122212222222122B22222222211222212122212X22222222222222222222222GCM矩陣疏水矩陣RKDEBZSNQGXTHACMPVLIYFWR1010998866655555433333210K1010998866655555433333210D9910108876665555544433321E9910108876665555544433321B8888101088887777666555443Z8888101088887777666555443S667788101010109999887777664N666688101010109999888777664Q666688101010109999888777664G556688101010109999888877665X555577999910101010998888775T555577999910101010998888775H555577999910101010999888775A555577999910101010999888775C4455668888999910109999885M334466888899991010101099887P33446678888899910101099987V3344557778888891010101010987L33335577778888999101010998I33335577778888999101010998Y2233446666777788999910108F1122446666777788889910109W001133444555556777888910PAM&BLOSOM

這類矩陣列出同源蛋白質(zhì)在進(jìn)化過程中氨基酸變化的可能性。這類矩陣是基于進(jìn)化原理的證據(jù)：編碼相同蛋白質(zhì)的基因隨著進(jìn)化發(fā)生分歧，相似度降低?？茖W(xué)用得多

PAM矩陣（pointacceptedmutaion）

基于氨基酸進(jìn)化的點突變模型

如果兩種氨基酸替換頻繁，說明自然界接受這種替換，那么這對氨基酸替換得分就高一個PAM就是一個進(jìn)化的變異單位,即1%的氨基酸改變

但這并不意味100次PAM后，每個氨基酸都發(fā)生變化，因為其中一些位置可能會經(jīng)過多次突變，甚至可能會變回到原來的氨基酸。PAM矩陣的制作步驟構(gòu)建序列相似（大于85％）的比對計算氨基酸j的相對突變率mj（j被其他氨基酸替換的次數(shù)）針對每個氨基酸對i和j,計算j

被i

替換次數(shù)替換次數(shù)除以相對突變率（mj）利用每個氨基酸出現(xiàn)的頻度對j

進(jìn)行標(biāo)準(zhǔn)化取常用對數(shù)，得到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針對不同的進(jìn)化距離采用PAM矩陣序列相似度=40%50%60%

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

PAM矩陣與BLOSUM矩陣的比較第三節(jié)雙序列比對算法序列的兩兩比對 （PairwiseSequenceAlignment）

按字符位置重組兩個序列，使得兩個序列接近一樣的長度序列兩兩比對基本算法直接方法——生成兩個序列所有可能的比對，分別計算代價函數(shù)，然后挑選一個代價最小的比對作為最終結(jié)果，需要計算2300

次——天文數(shù)字ATTC………CGAAGA

AGTC………GAAGGT假設(shè)比較300個氨基酸長度的兩條序列動態(tài)規(guī)劃方法DynamicProgramming起點終點ATTC………CGAAGA

AGTC………GAAGGTATTC………CGAAGAGTC………GAAGGAT+（1）ATTC………CGAAGAAGTC………GAAGG-T+（2）ATTC………CGAAGAGTC………GAAGGTA-+（3）最短路經(jīng)問題起點終點C1C2W1

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

取最小值！算法求解:

從起點到終點逐層計算

計算過程：

計算過程：按行計算其他方式

計算過程：（3）求最佳路徑算法分析： 數(shù)據(jù)結(jié)構(gòu)di,j

空間復(fù)雜度：O(mn)

時間復(fù)雜度：O(mn)由于在所考慮的子問題空間中，總共有θ(mn)個不同的子問題，因此，用動態(tài)規(guī)劃算法自底向上地計算最優(yōu)值能提高算法的效率。矩陣賦值算法

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)}

算法程序反向構(gòu)造匹配序列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目標(biāo)： 使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的后面出現(xiàn)時不會被罰分影響三、比對的統(tǒng)計學(xué)顯著性（1）典型方法：將兩條待比較的序列分別隨機(jī)打亂使用相同的程序與打分函數(shù)（或打分矩陣）進(jìn)行比對計算這些隨機(jī)序列的相似性得分重復(fù)這一過程（50～100次）用和分別表示其平均值與標(biāo)準(zhǔn)差。設(shè)原來兩條序列的比對得分為x，利用下式計算大于或等于x的比對得分概率：z=(x-)/根據(jù)z值判斷兩個序列相似得分的顯著性，當(dāng)z值是3.1、4.3、5.2時，x出現(xiàn)的概率為10-3、10-5、10-7Z>5，同源；Z<3,不同源；Z=3~5,可能同源經(jīng)驗法則（針對蛋白質(zhì)序列）：①如果兩個序列的長度都大于100，在適當(dāng)?shù)丶尤肟瘴恢螅鼈兣鋵Φ南嗤蔬_(dá)到25%以上，則兩個序列相關(guān)；②如果配對的相同率小于15%，則不管兩個序列的長度如何，它們都不可能相關(guān)；③如果兩個序列的相同率在15%25%之間，它們可能是相關(guān)的。數(shù)據(jù)庫的搜索簡介第四節(jié)雙序列比對工具

數(shù)據(jù)庫查詢?yōu)樯飳W(xué)研究提供了一個重要工具，在實際工作中經(jīng)常使用。然而，在分子生物學(xué)研究中，對于新測定的堿基序列或由此翻譯得到的氨基酸序列，往往需要通過數(shù)據(jù)庫搜索，找出具有一定相似性的同源序列，以推測該未知序列可能屬于哪個基因家族，具有哪些生物學(xué)功能。對于氨基酸序列來說，有可能找到已知三維結(jié)構(gòu)的同源蛋白質(zhì)而推測其可能的空間結(jié)構(gòu)。因此，數(shù)據(jù)庫搜索與數(shù)據(jù)庫查詢一樣，是生物信息學(xué)研究中的一個重要工具。

數(shù)據(jù)庫搜索的基礎(chǔ)是序列的相似性比對，即雙序列比對(pairwisealignment)。新測定的、希望通過數(shù)據(jù)庫搜索確定其性質(zhì)或功能的序列稱作檢測序列(probesequence)；通過數(shù)據(jù)庫搜索得到的和檢測序列具有一定相似性的序列稱目標(biāo)序列(subjectsequence)。為了確定檢測序列和一個已知基因家族之間的進(jìn)化關(guān)系，在通過數(shù)據(jù)庫搜索得到某些相似序列后，還需要判斷其序列相似性程度。如果檢測序列和目標(biāo)序列的相似性程度很低，還必須通過其他方法或?qū)嶒炇侄尾拍艽_定其是否屬于同一基因家族。一、BLAST簡介BLAST程序是目前最常用的基于局部相似性的數(shù)據(jù)庫搜索程序，它們都基于查找完全匹配的短小序列片段，并將它們延伸得到較長的相似性匹配。它們的優(yōu)勢在于可以在普通的計算機(jī)系統(tǒng)上運行，而不必依賴計算機(jī)硬件系統(tǒng)而解決運行速度問題。BLAST數(shù)據(jù)庫搜索策略BLAST僅通過部分而不是全部序列計算最適聯(lián)配值——贏得搜索速度

比對統(tǒng)計學(xué)意義的評價——E值(E-Value)P值(P-Value)(概率值)BLAST程序中使用了E值而非P值，這主要是從直觀和便于理解的角度考慮。比如E值等于5和10，總比P值等于0.993和0.99995更直觀。但是當(dāng)E<0.01時，P值與E值接近相同參數(shù)K和λ可分別被簡單地視為搜索步長(searchspacesize)和計分系統(tǒng)(scoringsystem)的特征數(shù)BLAST軟件包實際上是綜合在一起的一組程序，不僅可用于直接對蛋白質(zhì)序列數(shù)據(jù)庫和核酸序列數(shù)據(jù)庫進(jìn)行搜索，而且可以將檢測序列翻譯成蛋白質(zhì)或?qū)?shù)據(jù)庫翻譯成蛋白質(zhì)后再進(jìn)行搜索，以提高搜索結(jié)果的靈敏度。BLAST程序檢測序列和數(shù)據(jù)庫類型程序名檢測序列數(shù)據(jù)庫類型方法Blastp蛋白質(zhì)蛋白質(zhì)用檢測序列蛋白質(zhì)搜索蛋白質(zhì)序列數(shù)據(jù)庫Blastn核酸核酸用檢測序列核酸搜索核酸序列數(shù)據(jù)庫Blastx核酸蛋白質(zhì)將核酸序列按6條鏈翻譯成蛋白質(zhì)序列后搜索蛋白質(zhì)序列數(shù)據(jù)庫Tblastn蛋白質(zhì)核酸用檢測序列蛋白質(zhì)搜索由核酸序列數(shù)據(jù)庫按6條鏈翻譯成的蛋白質(zhì)序列數(shù)據(jù)庫Tblastx核酸核酸將核酸序列按6條鏈翻譯成蛋白質(zhì)序列后搜索由核酸序列數(shù)據(jù)庫按6條鏈翻譯成的蛋白質(zhì)序列數(shù)據(jù)庫

對一般用戶來說，目前常用的辦法是通過NCBI、EBI等國際著名生物信息中心的BLAST服務(wù)器進(jìn)行搜索。需要說明的是，各生物信息中心BLAST用戶界面有所不同，所提供的數(shù)據(jù)庫也可能不完全相同，使用前最好先進(jìn)行適當(dāng)?shù)倪x擇。BLAST應(yīng)用實例多結(jié)構(gòu)域蛋白(H1N1)

脂質(zhì)運載蛋白多結(jié)構(gòu)域蛋白(H1N1)的BLAST檢索>gi|224983683|pdb|3GBN|BChainB,CrystalStructureOfFabCr6261InComplexWithThe1918H1n1InfluenzaVirusHemagglutininGLFGAIAGFIEGGWTGMIDGWYGYHHQNEQGSGYAADQKSTQNAIDGITNKVNSVIEKMNTQFTAVGKEFNNLERRIENLNKKVDDGFLDIWTYNAELLVLLENERTLDFHDSNVRNLYEKVKSQLKNNAKEIGNGCFEFYHKCDDACMESVRNGTYDYPKYSEESKLNREEIDGVSGR多結(jié)構(gòu)域蛋白(H1N1)的BLAST檢索H1N1聚合酶序列BLAST結(jié)果綜述BLAST結(jié)果表述BLAST結(jié)果逐條顯示BLAST結(jié)果逐條顯示BLAST：改變打分矩陣的作用脂質(zhì)運載蛋白序列>sp|P31025|LCN1_HUMANLipocalin-1OS=HomosapiensGN=LCN1PE=1SV=1MKPLLLAVSLGLIAALQAHHLLASDEEIQDVSGTWYLKAMTVDREFPEMNLESVTPMTLTTLEGGNLEAKVTMLISGRCQEVKAVLEKTDEPGKYTADGGKHVAYIIRSHVKDHYIFYCEGELHGKPVRGVKLVGRDPKNNLEALEDFEKAAGARGLSTESILIPRQSETCSPGSD使用Blosum62矩陣搜索使用PAM30矩陣搜索FastA簡介FastA算法是由Lipman和Pearson于1985年發(fā)表的（Lipman和Pearson，1985）。FastA的基本思路是識別與代查序列相匹配的很短的序列片段，稱為k-tuple。

蛋白質(zhì)序列數(shù)據(jù)庫搜索時，短片段的長度一般是1～2個殘基長；DNA序列數(shù)據(jù)庫搜索時，通常采用稍大點的值，最多為6個堿基。通過比較兩個序列中的短片段及其相對位置，可以構(gòu)成一個動態(tài)規(guī)劃矩陣的對角線方向上的一些匹配片段。

FastA程序采用漸進(jìn)（heuristicapproach）算法將位于同一對角線上相互接近的短片段連接起來。也就是說，通過不匹配的殘基將這些匹配殘基片段連接起來，以便得到較長的相似性片段。這就意味著，F(xiàn)astA輸出結(jié)果中允許出現(xiàn)不匹配殘基。這和BLAST程序中的成對片段類似。如果匹配區(qū)域很多，F(xiàn)astA利用動態(tài)規(guī)劃算法在這些匹配區(qū)域間插入空位。

由FastA搜索產(chǎn)生的典型輸出結(jié)果的第一行列出程序名稱和版本號，以及該程序發(fā)表的雜志。接下來列出所提交的序列，然后是所用參數(shù)和運行時間，緊跟這些一般信息的是數(shù)據(jù)庫搜索結(jié)果。首先列出搜索得到的目標(biāo)序列簡單說明，其數(shù)目可由用戶定義。所列出的目標(biāo)序列的信息包括：序列所在數(shù)據(jù)庫名稱的縮寫，目標(biāo)序列的標(biāo)識碼、序列號和序列名等部分信息。括號中標(biāo)明匹配部分的殘基數(shù)。緊接著是由程序計算得到的初始化和優(yōu)化后的分?jǐn)?shù)值。最后一列是期望值即E值，用來判斷比對結(jié)果的置信度。接近于0的E值表明兩序列的匹配不大可能是由隨機(jī)因素造成的。以兩條氨基酸序列的比較為例介紹算法的基本思路，算法可以分為4步：第一步：FASTA首先找出進(jìn)行比較的兩條序列所有長度為k-tuple的連續(xù)的一致序列片段。例如以下兩條蛋白質(zhì)序列：··設(shè)k-tupl