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

第二章雙序列比對(duì)

PairwiseSequenceAlignment

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

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

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

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

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

（hydrophobicmatrix)（iv）PAM矩陣（pointacceptedmatrix,PAM）（v）BLOSUM矩陣（BLOckSUbstitutionMatrix，BLOSUM）其中Rij代表打分矩陣元素i、j分別代表字母表第i和第j個(gè)字符。遺傳密碼矩陣遺傳密碼矩陣通過計(jì)算一個(gè)氨基酸變成另一個(gè)氨基酸所需的密碼子變化的數(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)化的點(diǎn)突變模型

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

但這并不意味100次PAM后，每個(gè)氨基酸都發(fā)生變化，因?yàn)槠渲幸恍┪恢每赡軙?huì)經(jīng)過多次突變，甚至可能會(huì)變回到原來的氨基酸。PAM矩陣的制作步驟構(gòu)建序列相似（大于85％）的比對(duì)計(jì)算氨基酸j的相對(duì)突變率mj（j被其他氨基酸替換的次數(shù)）針對(duì)每個(gè)氨基酸對(duì)i和j,計(jì)算j

被i

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

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

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

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

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

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

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

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

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

取最小值！算法求解:

從起點(diǎn)到終點(diǎn)逐層計(jì)算

計(jì)算過程：

計(jì)算過程：按行計(jì)算其他方式

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

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

時(shí)間復(fù)雜度：O(mn)由于在所考慮的子問題空間中，總共有θ(mn)個(gè)不同的子問題，因此，用動(dòng)態(tài)規(guī)劃算法自底向上地計(jì)算最優(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}}子序列與完整序列的比對(duì)

----AGCT----ATGCAGCTGCTT目標(biāo)： 使S(s,i:t:j)最大序列S:序列t:ij不計(jì)前綴0:t:i的得分,也不計(jì)刪除后綴的j+1:t:|t|得分不計(jì)刪除后綴的j+1:t:|t|得分

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

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

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

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

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

脂質(zhì)運(yùn)載蛋白多結(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ì)運(yùn)載蛋白序列>sp|P31025|LCN1_HUMANLipocalin-1OS=HomosapiensGN=LCN1PE=1SV=1MKPLLLAVSLGLIAALQAHHLLASDEEIQDVSGTWYLKAMTVDREFPEMNLESVTPMTLTTLEGGNLEAKVTMLISGRCQEVKAVLEKTDEPGKYTADGGKHVAYIIRSHVKDHYIFYCEGELHGKPVRGVKLVGRDPKNNLEALEDFEKAAGARGLSTESILIPRQSETCSPGSD使用Blosum62矩陣搜索使用PAM30矩陣搜索FastA簡(jiǎn)介FastA算法是由Lipman和Pearson于1985年發(fā)表的（Lipman和Pearson，1985）。FastA的基本思路是識(shí)別與代查序列相匹配的很短的序列片段，稱為k-tuple。

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

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

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