生物序列聯(lián)配中的算法-中科院計算所生物信息學實驗室

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

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

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

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

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

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

7、－定義,定義1：兩個任意的字符 x和y，?(x,y)表示表x和y比較時的分值。 ?(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),動態(tài)規(guī)劃DP(Dynamic Programming)Smith-Waterman 算法

9、計算出兩個序列的相似分值，存于一個矩陣中。(相似度矩陣、DP矩陣)根據(jù)此矩陣，按照動態(tài)規(guī)劃的方法尋找最優(yōu)的聯(lián)配序列。,全局聯(lián)配(4),前提條件遞歸關(guān)系,,,全局聯(lián)配(5),在得到相似度矩陣后，通過動態(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)具有很高的相似度。在生物學中局部

11、聯(lián)配比全局聯(lián)配更具有實際的意義。兩條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),對全局聯(lián)配策略稍作修改可得到局部最優(yōu)聯(lián)配算法。聯(lián)配的路徑不需要到達搜索圖的盡頭

12、 ,如果某種聯(lián)配的分值不會因為增加聯(lián)配的數(shù)量而增加時，這種聯(lián)配就是最佳的。依賴于記分系統(tǒng)的性質(zhì)：因為某種路徑的記分會在不匹配的序列段減少 ,當分值降為零時，路徑的延展將會終止，一個新的路徑就會產(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ù)的盡可能長的空格。 空位的引入是為了補償那些插入或缺失，但是在序列的聯(lián)配中引入的空位不能太多，否則會使序列的排列變得面目全非。每

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

32、的并行算法,利用預測計算（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條序列，任意分成兩組。然后開始計算這兩組序列之間的兩條序列的聯(lián)配分值，在這一步中新的空位的位置gap_positions被保存下來以備第三步調(diào)用 。,第二步：設(shè)置一個標記flag，如果一個新的聯(lián)配的分值高于當前的最佳聯(lián)配隊列的分值，則flag標記為A( Accepted )，否則標記為R( Re

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