生物序列聯(lián)配中的算法-中科院計(jì)算所生物信息學(xué)實(shí)驗(yàn)室

1、生物序列聯(lián)配中的算法,張 法,提 綱,背景知識(shí)序列相似性的比較兩條序列的聯(lián)配問題多序列的聯(lián)配問題一些啟發(fā)式的算法生物序列聯(lián)配中的并行算法,DNA(1) 脫氧核糖核酸,DNA的分子組成核甘(nucleotides)磷酸鹽(phosphate)糖(sugar)一種堿基腺嘌呤(Adenine)鳥嘌呤(Guanine)胞嘧啶(Cytosine)胸腺嘧啶(Thymine),DNA(2),堿基的配對(duì)原則A(腺嘌呤)—T(

2、胸腺嘧啶)C(鳥嘌呤)—G(胞嘧啶) 一個(gè)嘌呤基與一個(gè)嘧啶基通 過氫鍵聯(lián)結(jié)成一個(gè)堿基對(duì)。DNA分子的方向性5'→3',DNA(3),DNA的雙螺旋結(jié)構(gòu) 堿基對(duì)之間的互補(bǔ)能力,DNA(4),DNA的復(fù)制在DNA解旋酶的作用 下兩條鏈分離開，分 別作為一個(gè)模板，在 聚合酶的作用下合成 一條新鏈。,RNA、轉(zhuǎn)錄和翻譯,RNA(核糖核酸)：單鏈結(jié)構(gòu)、尿嘧啶U代替胸腺嘧啶T

3、、位于細(xì)胞核和細(xì)胞質(zhì)中。轉(zhuǎn)錄：DNA鏈 → RNA鏈 信使RNA(mRNA)，啟動(dòng)子。翻譯： mRNA上攜帶遺傳信息在核糖體中合成蛋白質(zhì)的過程。,變異,進(jìn)化過程中由于不正確的復(fù)制，使DNA內(nèi)容發(fā)生局部的改變。 變異的種類主要有以下三種： 替代(substitution)插入或刪除(insertion or deletion) indel重排(rearrangement),蛋白質(zhì),由氨基酸依次

4、鏈接形成在生物體中總共有20種氨基酸。蛋白有十分復(fù)雜的三維結(jié)構(gòu)。其三維機(jī)構(gòu)決定了蛋白質(zhì)的功能。,基 因,什么是基因？DNA上具有特定功能的一個(gè)片斷，負(fù)責(zé)一種特定性狀的表達(dá)。一般來講，一個(gè)基因只編碼一個(gè)蛋白質(zhì)。,基因組,任何一條染色體上都帶有許多基因，一條高等生物的染色體上可能帶有成千上萬個(gè)基因，一個(gè)細(xì)胞中的全部基因序列及其間隔序列統(tǒng)稱為genomes（基因組）。,DNA上的基因,基因,,,,,基因的編碼,基因編碼是一個(gè)邏輯的映射，表

5、明存儲(chǔ)在DNA和mRNA中的基因信息決定什么樣的蛋白質(zhì)序列。每個(gè)堿基三元組稱為一個(gè)密碼子(codon)堿基組成的三元組的排列共有43＝64種，而氨基酸共有20種類型，所以不同的密碼子可能表示同一種氨基酸。,帶來的問題,序列排列問題 基因組的重排問題 蛋白質(zhì)結(jié)構(gòu)和功能的預(yù)測(cè) 基因(外顯子、內(nèi)含子)查找問題 序列裝配(Sequence Assembly)問題,生物序列相似性的比較,動(dòng)機(jī),在生物學(xué)的研究中，將未知序列同已知序列進(jìn)行

6、比較分析已經(jīng)成為一種強(qiáng)有力的研究手段 ，生物學(xué)領(lǐng)域中絕大部分的問題在計(jì)算機(jī)科學(xué)領(lǐng)域中主要體現(xiàn)為序列或字符串的問題 。,序列聯(lián)配問題的分類,如果兩個(gè)序列具有足夠的相似性，則認(rèn)為兩者具有同源性。－序列相似性的比較 (兩條序列的聯(lián)配)序列的分類序列的排列多序列的聯(lián)配,,,兩條序列聯(lián)配問題的分類,全局聯(lián)配(Global Alignment)局部聯(lián)配(Local Alignment)空位處罰(Gap Penalty),全局聯(lián)配(1)

7、－定義,定義1：兩個(gè)任意的字符 x和y，?(x,y)表示表x和y比較時(shí)的分值。 ?(x,x)=2, ?(x,y)= ?(x,-)= ?(-,y)=-1定義2：S= s1…sn和T=t1…tm，其全局聯(lián)配A可以用序列S’和T’來表示，其中： (1) | S’ | = | T’ |； (2) 將S’和T’中的空字符除去后所得到的序列分別為S和T;聯(lián)配A的分值Score為：,,全局聯(lián)配(2)－原始算法,輸入：序

8、列S和T，其中 | S | = | T | = n 輸出：S和T的最優(yōu)聯(lián)配 for i=0 to n do for (S的所有的子序列A，其中| A | = i ) do for (T的所有的子序列B，其中| B | = i ) do ……,全局聯(lián)配(3),動(dòng)態(tài)規(guī)劃DP(Dynamic Programming)Smith-Waterman 算法

9、計(jì)算出兩個(gè)序列的相似分值，存于一個(gè)矩陣中。(相似度矩陣、DP矩陣)根據(jù)此矩陣，按照動(dòng)態(tài)規(guī)劃的方法尋找最優(yōu)的聯(lián)配序列。,全局聯(lián)配(4),前提條件遞歸關(guān)系,,,全局聯(lián)配(5),在得到相似度矩陣后，通過動(dòng)態(tài)規(guī)劃回溯（Traceback）的方法可獲得序列的最優(yōu)聯(lián)配序列 。例： S = “a c g c t g”和T = “c a t g t” ?(x,x)=2, ?(x,y)= ?(x,-)= ?(-,y)=-

10、1,,三種可能的最優(yōu)聯(lián)配序列：S: a c g c t g - T: - c – a t g t S: a c g c t g - T: - c a – t g t S: - a c g c t g T: c a t g - t -,局部聯(lián)配(1),兩條序列在一些局部的區(qū)域內(nèi)具有很高的相似度。在生物學(xué)中局部

11、聯(lián)配比全局聯(lián)配更具有實(shí)際的意義。兩條DNA長序列，可能只在很小的區(qū)域內(nèi)(密碼區(qū))存在關(guān)系。不同家族的蛋白質(zhì)往往具有功能和結(jié)構(gòu)上的相同的一些區(qū)域。,局部聯(lián)配(2),前提條件： V(i, 0) = 0; V(0, j) = 0;遞歸關(guān)系： 找出i*和j*，使得：,,,局部聯(lián)配(3),對(duì)全局聯(lián)配策略稍作修改可得到局部最優(yōu)聯(lián)配算法。聯(lián)配的路徑不需要到達(dá)搜索圖的盡頭

12、 ,如果某種聯(lián)配的分值不會(huì)因?yàn)樵黾勇?lián)配的數(shù)量而增加時(shí)，這種聯(lián)配就是最佳的。依賴于記分系統(tǒng)的性質(zhì)：因?yàn)槟撤N路徑的記分會(huì)在不匹配的序列段減少 ,當(dāng)分值降為零時(shí)，路徑的延展將會(huì)終止，一個(gè)新的路徑就會(huì)產(chǎn)生。,局部聯(lián)配(4),S = “ a b c x d e x ”，T= “ x x x c d e ” 記分函數(shù)?(x,y)： ?(x,x)=2, ?(x,y)= ?(x,-)= ?(-,y)=-1,,,S = “ a b

13、 c x d e x ”，T= “ x x x c d e ” 局部最優(yōu)聯(lián)配是： c x d e c - d e或 x - d e x c d e,,空位處罰(1),空位:序列中任意連續(xù)的盡可能長的空格。 空位的引入是為了補(bǔ)償那些插入或缺失，但是在序列的聯(lián)配中引入的空位不能太多，否則會(huì)使序列的排列變得面目全非。每

14、引入一個(gè)空位，聯(lián)配的分值都會(huì)有所扣除，常見的罰分規(guī)則主要有兩種：空位權(quán)值恒定模型仿射空位處罰模型,空位處罰(2),空位權(quán)值恒定模型： 在每個(gè)空位中的空格的分值為零， 即：?(x,-)= ?(-,y) = 0 聯(lián)配的分值為：其中，Wg為開放一個(gè)空位的罰分,,空位處罰(3),仿射空位處罰模型(Affine Gap Model) 用一個(gè)附加的罰分比例去乘空位的長度

15、 Wg+q×Ws 聯(lián)配的分值為：,,空位處罰(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[

16、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}.,利用動(dòng)態(tài)規(guī)劃計(jì)算序列最優(yōu)聯(lián)配的算法的復(fù)雜度分析：時(shí)間復(fù)雜度；O(nm) 空間復(fù)雜度：O(n+m),,最新的相關(guān)的研究,在動(dòng)態(tài)規(guī)劃算法的基礎(chǔ)上提出了一種新的方法，解決在序列局部聯(lián)配的最優(yōu)排列中經(jīng)常出現(xiàn)的

17、馬賽克問題(在最優(yōu)排列中間經(jīng)常出現(xiàn)的相似度很低的保守區(qū)域)把hash表方法引入到基因組序列的局部聯(lián)配問題中，同時(shí)提高了原有算法的效率和質(zhì)量,多序列聯(lián)配(1),兩條序列聯(lián)配問題的一般化推廣動(dòng)機(jī)：攜帶更多的消息、生物學(xué)家的一種重要手段。 DNA或蛋白質(zhì)數(shù)據(jù)庫容量的指數(shù)級(jí)的增長，即使使用很簡單的記分函數(shù)也很難找到一種在有效時(shí)間內(nèi)的解決方法，而幾乎所有的近似算法和許多的啟發(fā)式算法，實(shí)際上都是在算法的計(jì)算速度和獲得最佳聯(lián)配結(jié)果的敏感性之間

18、尋找一種權(quán)衡（tradeoff）策略。,多序列聯(lián)配(2),rigorous算法兩條序列 多條序列， NP hardtree_based算法和 iterative算法首先從序列中找出兩條相似度最大的聯(lián)配，然后按照相似度遞減的順序連續(xù)添加一些序列到最優(yōu)的聯(lián)配序列中。序列非常接近或?qū)儆谝粋€(gè)同源的家族 原始算法基礎(chǔ)上的近似算法，但是它們也是非常耗時(shí)的。,,多序列聯(lián)配(3),記分方法：用函數(shù)δ(x, y)來計(jì)算字符x

19、和y之間的距離，兩個(gè)序列的聯(lián)配的距離分值我們用V來表示： k條序列聯(lián)配的分值：為k條序列中任意兩條序列(共有Ck2條)的分值(距離)V之和，用SP (Sum_of_Pairs)來表示:,,,多序列聯(lián)配(4),中心星算法輸入：一組含有k條序列的集合? 找出序列St，St∈?，使得∑i≠tD( Si, St)的值最小，令A(yù) = { St } 逐次地添加Si∈?-{St}到A中，并使Si與St的聯(lián)配的值最小；假設(shè)S1，S

20、2，…Si-1已添加到A中，由于在分別和St進(jìn)行聯(lián)配的過程中需要加入一些空格，故此時(shí)A為：A = {S1’，S2’，…Si-1’，S’t}。添加Si到A中，按照兩條序列聯(lián)配的動(dòng)態(tài)規(guī)劃算法比較S’t和Si，分別產(chǎn)生新的序列St’’和Si’，再按照St’’中添加空格的位置調(diào)節(jié)序列S1’，S2’，…Si+1’為S1’’，S2’’，…Si-1’’，并用St’’替換St。,借助硬件來完成動(dòng)態(tài)規(guī)劃的算法 采用并行的固件，把問題有效地分配到多個(gè)處理

21、器上運(yùn)行，最后再把各處理器的運(yùn)算結(jié)果收集起來 借助于一些比最初的動(dòng)態(tài)規(guī)劃算法更快的啟發(fā)式算法。以降低運(yùn)算結(jié)果的質(zhì)量為代價(jià)來提高計(jì)算速度的。,啟發(fā)式算法-FASTA,基本思想是：一個(gè)能夠揭示出真實(shí)的序列關(guān)系的聯(lián)配至少包含一個(gè)兩個(gè)序列都擁有的字(片斷)，把查詢序列中的所用字編成索引，然后在數(shù)據(jù)庫搜索時(shí)查詢這些索引，以檢索出可能的匹配，這樣那些命中的字很快被鑒定出來。,確定參數(shù)ktup，在兩個(gè)序列中查找長度為ktup的、相匹配的片段(增強(qiáng)點(diǎn)

22、)。連續(xù)的增強(qiáng)點(diǎn)在動(dòng)態(tài)規(guī)劃矩陣中主對(duì)角線附近。為了提高速度，可以通過查詢表格或hash表來完成，然后在表格中搜索與另一條序列相匹配的、長度為ktup的片段。在同一條對(duì)角線中臨近的增強(qiáng)點(diǎn)成為一個(gè)增強(qiáng)段。每一個(gè)增強(qiáng)點(diǎn)都賦予一個(gè)正的分值，一個(gè)增強(qiáng)段中相鄰的兩個(gè)增強(qiáng)點(diǎn)之間的不匹配區(qū)域賦予一定的負(fù)值。一個(gè)增強(qiáng)段對(duì)應(yīng)于一段相匹配的子序列，分值最高的段被標(biāo)記為init1。,引入indel。把那些沒有重疊(non-overlap)的增強(qiáng)段拼接起來(增

23、強(qiáng)段的分值之和減去空位處罰)。分值最高的區(qū)域記為initn。對(duì)最有可能的匹配序列進(jìn)一步評(píng)分：以增強(qiáng)段init1所在的對(duì)角線為中心，劃分出一個(gè)較狹窄的對(duì)角線帶，利用S-W算法，來獲得分值最高的局部聯(lián)配，記作opt。按照initn和opt的分值，對(duì)分值最高的序列再進(jìn)行一次S-W聯(lián)配。FASTA對(duì)每一個(gè)檢索到的聯(lián)配都提供一個(gè)統(tǒng)計(jì)學(xué)顯著性的評(píng)估，以判斷該聯(lián)配的意義。,FASTA對(duì)DNA序列搜索的結(jié)果要比對(duì)蛋白質(zhì)序列搜索的結(jié)果更敏感，因?yàn)樗鼘?duì)

24、數(shù)據(jù)庫的每一次搜索都只有一個(gè)最佳的聯(lián)配，這樣對(duì)于蛋白質(zhì)序列而言，一些有意義的聯(lián)配可能被錯(cuò)過。,啟發(fā)式算法-BLAST,基本思想是：通過產(chǎn)生數(shù)量更少的但質(zhì)量更好的增強(qiáng)點(diǎn)來提高速度。BALST算法是建立在嚴(yán)格的統(tǒng)計(jì)學(xué)的基礎(chǔ)之上的。它集中于發(fā)現(xiàn)具有較高的相似性的局部聯(lián)配，且局部聯(lián)配中不能含有空位。由于局部聯(lián)配的限制條件，在大多數(shù)情況下聯(lián)配會(huì)被分解為若干個(gè)明顯的HSP(High-score Sequence Pairs)。,首先確定一個(gè)終止

25、值S、步長參數(shù)w和一個(gè)閥值t ，S值通常是基于統(tǒng)計(jì)學(xué)的原理指明一個(gè)預(yù)期的終止E值，然后軟件會(huì)在考慮搜索背景性質(zhì)的基礎(chǔ)上計(jì)算出合適的S值。使要聯(lián)配的序列中包含一個(gè)分值不小于S的HSP。引入鄰近字串的思想：不需要字串確切地匹配，當(dāng)有一個(gè)字串的分值高于t時(shí)，BALST就宣稱找到了一個(gè)選中的字串。為了提高速度，允許較長的字串長度W。W值很少變化，這樣，t值就成為權(quán)衡速度和敏感度的參數(shù)。,,一個(gè)字串選中后，程序會(huì)進(jìn)行沒有空位的局部尋優(yōu)，聯(lián)配的最

26、低分值是S，當(dāng)聯(lián)配延伸時(shí)會(huì)遇到一些負(fù)的分值，使得聯(lián)配的分值下降，當(dāng)下降的分值小于S時(shí)，命中的延伸就會(huì)終止。這樣系統(tǒng)會(huì)減少消耗于毫無指望的選中延伸的時(shí)間，使系統(tǒng)的性能得以改進(jìn)。,,PSI-BLAST,Position Specific Iterated BLAST在1997年提出了對(duì)BLAST程序的改進(jìn)算法，在序列的聯(lián)配中允許出現(xiàn)空位，提高了搜索速度、敏感度和實(shí)用性。對(duì)一個(gè)選中字串長度標(biāo)準(zhǔn)的延伸 利用profile(表頭文件)的數(shù)據(jù)

27、結(jié)構(gòu)來進(jìn)行搜索,在利用動(dòng)態(tài)規(guī)劃矩陣進(jìn)行聯(lián)配時(shí)，以步長為2w來搜索每一條對(duì)角線，這樣可以找到一些長度為2w的選中的字。 改進(jìn)的算法中允許位于不同的對(duì)角線的兩個(gè)片段拼接在一起。在拼接的過程中要涉及到indel操作，與FASTA算法只產(chǎn)生一個(gè)最佳聯(lián)配不同，位于不同對(duì)角線的兩個(gè)片段拼接在一起的前提條件是：拼接后片段的分值不小于某一個(gè)終止值。,,執(zhí)行通常的BLAST算法，使用一種不同的記分方式，根據(jù)高度顯著聯(lián)配(HSPs)的最高分值建立一個(gè)最初

28、的profile。 根據(jù)該profile反復(fù)利用BLAST算法對(duì)數(shù)據(jù)庫進(jìn)行搜索，這一步實(shí)際上是根據(jù)表頭文件的統(tǒng)計(jì)結(jié)果擴(kuò)展局部聯(lián)配。這一過程是反復(fù)進(jìn)行的，直到再?zèng)]有發(fā)現(xiàn)新的有意義的匹配為止。由于在每一輪都會(huì)有新的片段加入，因此在操作過程中profile需要在每一個(gè)循環(huán)結(jié)束之后更新。,生物序列聯(lián)配中的并行算法,兩條序列聯(lián)配的并行算法,根據(jù)序列的相似性比較，找出兩者的最佳匹配 找出從一條序列轉(zhuǎn)化到另一條序列的最佳路徑核心：動(dòng)態(tài)規(guī)劃,L

29、ander的處理方法,如果已知第i-1和 i-2行反對(duì)角線上 的各項(xiàng)值，那么 第i行反對(duì)角線上 各項(xiàng)的值可以并 行計(jì)算。,S-W的并行處理,Row-Wavefront算法和Diagonal Wavefront算法,Row-Wavefront算法是讓每個(gè)處理器順序處理動(dòng)態(tài)規(guī)劃矩陣中的每一行，當(dāng)一個(gè)處理器檢測(cè)到它所需要的上一行矩陣的元素還沒有計(jì)算出來時(shí)，該處理器就阻塞自己 ，直到所需要的元素計(jì)算出來 。,S-

30、W的并行處理,Row-Wavefront算法和Diagonal Wavefront算法,Diagonal Wavefront算法中所有的處理器同時(shí)沿著矩陣的一條反對(duì)角線進(jìn)行計(jì)算，當(dāng)一條反對(duì)角線的元素全部計(jì)算完，才轉(zhuǎn)到下一條反對(duì)角線計(jì)算。當(dāng)一個(gè)處理器空閑的時(shí)候，它從當(dāng)前的反對(duì)角線上尋找一個(gè)還沒有計(jì)算的元素進(jìn)行計(jì)算。這種算法沒有嚴(yán)格的處理器計(jì)算順序,Huang的處理方法,采用了分而治之的策略對(duì)回溯算法進(jìn)行了改進(jìn) :按照中間的一條反對(duì)角線來分

31、割動(dòng)態(tài)規(guī)劃矩陣。,動(dòng)態(tài)規(guī)劃的并行計(jì)算,基于流水線的動(dòng)態(tài)規(guī)劃算法反對(duì)角線的動(dòng)態(tài)規(guī)劃算法 反對(duì)角線分塊的動(dòng)態(tài)規(guī)劃算法 粗粒度分塊策略 細(xì)粒度分塊策略H-V(Horizontal-Vertical )分塊策略,基于流水線的動(dòng)態(tài)規(guī)劃算法,,反對(duì)角線分塊的動(dòng)態(tài)規(guī)劃算法,粗粒度分塊策略,,細(xì)粒度分塊策略,,H-V(Horizontal-Vertical )分塊策略,利用前趨計(jì)算的并行算法,Prefix Computation,,多序列聯(lián)配

32、的并行算法,利用預(yù)測(cè)計(jì)算（Speculative computation）技術(shù)的并行算法 分而治之策略的并行算法,Berger_Munson算法,best_score=initial_score();while(stop criteria is not met){ 1 current_score=calculate(seq,gap_positions); 2 flag=decide(current_sco

33、re, best_score); 3 seq=modify(seq,flag,gap_positions);},第一步：首先輸入n條序列，任意分成兩組。然后開始計(jì)算這兩組序列之間的兩條序列的聯(lián)配分值，在這一步中新的空位的位置gap_positions被保存下來以備第三步調(diào)用 。,第二步：設(shè)置一個(gè)標(biāo)記flag，如果一個(gè)新的聯(lián)配的分值高于當(dāng)前的最佳聯(lián)配隊(duì)列的分值，則flag標(biāo)記為A( Accepted )，否則標(biāo)記為R( Re

34、jected )；,第三步：如果在第二步中flag標(biāo)記為A，則根據(jù)第一步中的gap_positions來修改當(dāng)前的聯(lián)配隊(duì)列，以作為下一次迭代的初始值，同時(shí)更新聯(lián)配隊(duì)列的分值，當(dāng)連續(xù)遇到q個(gè)R時(shí)，算法結(jié)束。,Berger_Munson算法的并行化,每一次迭代內(nèi)的三個(gè)階段是相互依賴；第i次迭代可能會(huì)用到第i-1次迭代的結(jié)果。預(yù)測(cè)計(jì)算的基本思想是：利用當(dāng)前狀態(tài)的輸入?yún)?shù)來推斷未來狀態(tài)的計(jì)算結(jié)果。如果有連續(xù)的迭代被標(biāo)記為R，那么這些迭代之間