序列分析與聯(lián)配課件

第三章序列分析與聯(lián)配第一節(jié)序列組成和單一序列分析第二節(jié)序列聯(lián)配第三節(jié)數(shù)據(jù)庫搜索引擎——BLAST和FASTA應(yīng)用第四節(jié)寡核苷酸設(shè)計序列分析是生物信息學(xué)最主要的研究內(nèi)容之一，它可以分為兩個主要局部：一是序列組成〔特別是涉及到基因組層次上〕分析，二是序列之間的比較分析。兩條序列或多條序列間的比對或聯(lián)配(alignment)的目的，是對它們的序列相似性進(jìn)行評估，找出這些序列中結(jié)構(gòu)或功能相似性區(qū)域等。通過聯(lián)配未知序列與序列(其功能或結(jié)構(gòu)等)的相似程度，我們可以判斷或推測未知序列的結(jié)構(gòu)與功能。第一節(jié)序列組成及單一序列分析一、堿基組成二、堿基相鄰頻率三、同向重復(fù)序列分析四、DNA序列的幾何學(xué)分析——Z曲線一、堿基組成DNA序列一個顯而易見的特征是四種堿基類型的分布。盡管四種堿基的頻率相等時對數(shù)學(xué)模型的建立可能是方便的，但幾乎所有的研究都證明堿基是以不同頻率分布的。表3.1包含了9條完整DNA分子序列的資料，表3.2的數(shù)據(jù)來自兩個胎兒球蛋白基因(Gr和Ar)，每個基因具有三個外顯子和兩個內(nèi)含子(shen等1981)。這兩個例子說明序列內(nèi)和序列間堿基具有不同的頻率。在基因每一側(cè)的500個任意堿基區(qū)域被稱為“側(cè)翼〞，基因間區(qū)域是指兩個基因間的其余序列。二、堿基相鄰頻率分析DNA序列的主要困難之一是堿基相鄰的頻率不是獨(dú)立的。堿基相鄰的頻率一般不等于單個堿基頻率的乘積：如果Pu是序列中堿基u的頻率，且Puv為兩個相鄰堿基u和v的頻率，那么Puv≠PuPvNussinov(1984)研究了兩堿基相鄰的頻率(表3.3)。數(shù)據(jù)來自166個脊椎動物的DNA序列，總長136731個堿基。表中的比值為16種二個堿基相鄰的頻率除以相應(yīng)的單個堿基頻率的乘積。作為一個特別的例子，圖3.1給出了雞血紅蛋白β鏈的mRNA編碼區(qū)的438個堿基。表3.4列出了4種堿基和16種兩堿基的數(shù)目。將該表看作4×4的表，計算行列獨(dú)立性的卡方統(tǒng)計量，得到x2=59.3〔x20.05,9=16.92〕說明行(第一堿基)列(第二堿基)之間存在明顯的關(guān)聯(lián)。在編碼區(qū)，存在某種約束來限制DNA序列編碼氨基酸。在密碼子水平上，這一約束與堿基相鄰頻率有關(guān)。表3.5列出了遺傳密碼和圖3.1序列中各密碼子數(shù)量。盡管數(shù)目很小，難以作出有力的統(tǒng)計結(jié)論，但編碼同一氨基酸的不同密碼子(同義密碼子)好似不是等同存在的。這種密碼子偏倚必定與兩堿基相鄰頻率水平有關(guān)。表3.5還清楚地說明，由于密碼子第3位置上堿基的改變常常不會改變氨基酸的類型，因而對第3位置上堿基的約束要比第2位堿基小得多。相鄰堿基之間的關(guān)聯(lián)將導(dǎo)致更遠(yuǎn)堿基之間的關(guān)聯(lián)，這些關(guān)聯(lián)延伸距離的估計可以從馬爾科夫鏈(Markovchain)理論得到(Javare和Giddings，1989)。在不援引任何生物學(xué)機(jī)制的情況下，第k階馬爾科夫鏈假定在序列中某一位置上堿基的存在只取決于前面k個位置上的堿基。一階鏈假定一個特定堿基存在于位置i的概率只取取決于在位置i-1的4種堿基概率。相互獨(dú)立的堿基所組成的序列將與0階馬爾科夫鏈相對應(yīng)。階可以通過似然法估計。同時，馬爾科夫鏈分析更適應(yīng)于基因組水平，而非單一序列(基因)。三、同向重復(fù)序列分析

除了分析整個序列堿基關(guān)聯(lián)程度的特征外，我們常對尋找同向重復(fù)序列(directrepeats)之類的問題感興趣。Karlin等(1983)給出了完成這一分析的有效算法。該法采用由特定的幾組堿基字母組成的不同亞序列或稱為字碼(word)。只需要對整個序列搜索一次。給一堿基賦以值α,例如A、C、G、T的值為0、1、2、3。由X1、X2、...、Xk共k個字母組成的每一種不同的字碼按計算字碼值。這些值的取值范圍為1到4k。例如，5字碼TGACC的值為1+3×44+2×43+0×42+1×41+1×40=459?？上葟牡蚹值的字碼開始搜索。記錄序列中每一個位置k字碼的字碼值。只有在發(fā)現(xiàn)k字碼長度重復(fù)的那些位置考慮進(jìn)行長度大于k的字碼搜索。表3.6列出了序列TGGAAATAAAACGTAAGTAG中所有堿基2字碼(k=2)的初始位置和字碼值。對于完全重復(fù)、長度大于2的同向重復(fù)或亞序列的搜索可只限于2字碼重復(fù)的初始位置。在本例中只有4個重復(fù)的2堿基重復(fù)序列。例如，在位置4、5、8、9、10和15均發(fā)現(xiàn)了字碼值為1的堿基重復(fù)序列。從有重復(fù)的第2個堿基為起點(diǎn)的3字碼值及位置列于表3.7，其中發(fā)現(xiàn)字碼值為1、45和49的序列有重復(fù)。以每一重復(fù)的3堿基為起點(diǎn)的4字碼搜索未能發(fā)現(xiàn)更長的重復(fù)序列。因此最長的同向重復(fù)為4、8、9位置上的AAA，13、17位置上的GTA以及7、14位置上的TAA。同樣對圖3.1雞β球蛋白DNA序列進(jìn)行同向重復(fù)序列搜索，一些最長同向重復(fù)序列列于表3.8。Karlin等(1983)提出了序列內(nèi)存在的最長同向重復(fù)序列的統(tǒng)計顯著性評價方法。在核苷酸的位置為獨(dú)立的假定下(相當(dāng)于階次為0的馬爾科夫鏈)，長度為n的序列中，最長同向重復(fù)Ln的期望長度和方差為：其中，P為序列中堿基頻率的平方和：用盡可能接近最大長度的期望均值的字碼(即R→μL)來開始同向重復(fù)序列的搜索計算可能節(jié)省計算量。方程3.1可以用一個近似方法來驗證以上統(tǒng)計假說。假定同向重復(fù)序列的長度呈正態(tài)分布。對于圖3.1雞β蛋白序列，A、C、G、T四個堿基的次數(shù)分別為87、144、118和89，因而P=0.2614，最長重復(fù)序列的期望長度為8.13且具有期望方差0.9138。根據(jù)95%的正態(tài)分布概率，理論上可以預(yù)期最長同向重復(fù)序列不超過10。四、DNA序列的幾何學(xué)分析—Z曲線DNA序列實際上是一種用4種字母表達(dá)的“語言〞，只是其“詞法〞和“語法〞規(guī)那么目前還沒有搞清楚。人類的語言有文字、聲音兩種根本表現(xiàn)形式，此外還有手語、旗語甚至圖畫語等特殊表達(dá)形式。同樣，DNA序列作為一種語言，其表達(dá)形式也不是唯一的。傳統(tǒng)上，DNA序列是用4種字母符號表達(dá)的一維序列。這是一種抽象形式，適合于存儲、印刷和代數(shù)算法的處理，包括比較、排列和查找特殊序列等。我國學(xué)者張春霆等開展了DNA序列三維空間曲線表示形式，即DNA序列幾何表示形式的研究。幾何形式雖然與符號形式完全等價，但顯示了DNA序列的新特征。兩種形式各有其特點(diǎn)，相互補(bǔ)充。這一新方法，為解讀DNA序列信息提供了嶄新的手段。他們的研究始于對4種堿基對稱性的觀察，提出了用正面體表示堿基對稱性。1994年，他們利用這種形式來表示任意長度的DNA序列。現(xiàn)將這種序列表示方法簡述如下。考察一個長為L的單股DNA序列，方向(5'→3'或3'→5')不限。從第一個堿基開始，依次考察此序列，每次只考察一個堿基。當(dāng)考察到第n個堿基時(n=1,2，.,L)，數(shù)一下從1到n這個子序列中四種堿基各自出現(xiàn)的次數(shù)。設(shè)4種堿基A、C、G、T出現(xiàn)的次數(shù)分別以An、Cn、Gn、Tn表示之，這里下標(biāo)“n〞是說明這些整數(shù)是從1到n這個子序列中數(shù)出來的，如圖3.2所示。顯然，它們都是正整數(shù)。根據(jù)正四面體的對稱性可以證明，在正面體內(nèi)存在唯一的一個點(diǎn)Pn與這四個正整數(shù)對應(yīng)。點(diǎn)Pn構(gòu)成了四個正整數(shù)的一一對應(yīng)映射。點(diǎn)Pn坐標(biāo)可用四正整數(shù)表達(dá)：方程3.2其中xn，yn和zn為點(diǎn)Pn的三個坐標(biāo)分量。當(dāng)n從1跑到L時，我們依次得到P1，P2，P3，...，PL共L個點(diǎn)。將相鄰兩點(diǎn)用適當(dāng)?shù)那€連接所得到的整條曲線，就稱為表示DNA序列的Z曲線?？梢宰C明，Z曲線與所表示的DNA序列是一一對應(yīng)的，即給定一DNA序列，存在唯一的一條Z曲線與之對應(yīng)；反之，給定一條Z曲線，可找到唯一的一個DNA序列與之對應(yīng)。換言之，Z曲線包含了DNA序列的全部信息。Z曲線是與符號DNA序列等價的另一種表示形式，一種幾何形式?？梢酝ㄟ^Z曲線對DNA序列進(jìn)行研究。Z曲線的三個分量(方程3.2)具有明確的生物學(xué)意義：xn表示嘌呤/嘧啶堿基沿序列的分布。當(dāng)從1到n的這個子序列中(圖3.2)嘌呤堿基多于嘧啶堿基時，xn>0，否那么，xn<0，當(dāng)兩者相等時xn=0。同樣，yn表示氨基/酮基堿基沿序列的分布。當(dāng)在子序列中氨基堿基多于酮基堿基時，yn>0，否那么，yn<0，當(dāng)兩者相等時yn=0。zn表示強(qiáng)/弱氫鍵堿基沿序列的分布。當(dāng)弱氫鍵堿基多于強(qiáng)氫鍵堿基時，zn>0，否那么zn<0，當(dāng)兩者相等時，zn=0。這三種分布是相互獨(dú)立的，表現(xiàn)在以下事實上：任何一種分布不能由其它兩種分布的線性疊加表示出來。給定的DNA序列唯一地決定了這三種分布；三種分布唯一地描述了DNA序列。對DNA序列的研究就是通過對這三種分布的研究來進(jìn)行。從方法學(xué)的角度來看，這是DNA序列的一種幾何學(xué)研究途徑。圖3.3給出了大腸桿菌ayoP基因族序列Z曲線的三個分量，即三種分布圖。該基因族包含了大腸桿菌5個基因aroP，A，aceFE，aceF和lpd，總長度為9501bp，分別編碼芳香族氨基酸運(yùn)輸?shù)鞍譨roP，蛋白質(zhì)A(功能不詳)和三種酶，即丙酮酸脫氫酶，二氫硫辛?；D(zhuǎn)移酶和二氫硫辛酰脫氫酶。它們位于此序列的0039-1406，1947-2654，2870-5527，5545-7434，7759-9183區(qū)間。在圖中X軸的下方的基因排列圖上已分別用陰影標(biāo)出相應(yīng)基因。在這些基因之間有三個啟動子區(qū)〔pm1、pm2和pm3〕，其中aceE和aceF基因?qū)儆赼ce操縱子，共用一個啟動子。三個啟動子區(qū)亦在圖中標(biāo)出。非常令人感興趣的是，在5個編碼區(qū)，Z曲線的z分量根本上都是單調(diào)下降的，而在三個啟動子區(qū)根本上都是單調(diào)上升的。x，y分量亦有變化，但不如z分量明顯。在上升、下降的交界處，Z曲線均發(fā)生了重大的轉(zhuǎn)折，據(jù)此有可能用Z曲線識別這些位置。由此圖可見，用Z曲線這種幾何方法顯示DNA序列不僅直觀，而且作為一種識別序列中的不同基因和功能區(qū)的新方法，展現(xiàn)了廣闊的應(yīng)用前景。第二節(jié)序列聯(lián)配一、Needleman-Wunsch算法二、Smith-Waterman算法三、序列相似性統(tǒng)計特征四、替換矩陣五、多序列聯(lián)配一、Needleman-Wunsch算法有2種經(jīng)典方法可以計算兩條序列間的最適聯(lián)配。Needleman-Wunsch算法是一種整體聯(lián)配(globalalignment)算法，最正確聯(lián)配中包括了全部的最短匹配序列。Smith-Wateman算法是在Needleman-Wunsch算法根底上開展而來的，它是一種局部聯(lián)配(Localalignment)算法。這二種算法均可以用于核酸和蛋白質(zhì)序列。在給定空位罰值和替換矩陣情況下，它們總是能給出具有最高〔優(yōu)〕聯(lián)配值的聯(lián)配。但是，這個聯(lián)配并不需要到達(dá)生物學(xué)意義上的顯著水平。GCG軟件包中，BESFIT和GAP程序，EMBOSS的needle等可用于該聯(lián)配。一些網(wǎng)站可以通過遞交序列進(jìn)行兩條序列的聯(lián)配分析。從整體上分析兩個序列的關(guān)系，即考慮序列總長的整體比較，用類似于使整體相似(globalsimilarity)最大化的方式，對序列進(jìn)行聯(lián)配。兩個不等長度序列的聯(lián)配分析必需考慮在一個序列中圈掉一些堿基或在另一序列作空位(gap)處理。Needleman和Wunsch(1970)的法那么為這些步驟提供了實例。這一算法是為氨基酸序列開展的，但也可以用于核苷酸序列。算法最初尋求的是使兩條序列間的距離最小。盡管這類距離的元素是以一種特定的方式定義的，但該算法的良好特性在于它確定了最短距離。這是一個動態(tài)規(guī)劃(dynamicprogramming)的方法。將兩條聯(lián)配的序列沿雙向表的軸放置，兩條序列的所有可能的聯(lián)配方式都將在它們所形成的方形圖中〔見以下圖〕。從任一堿基對，即表中的任一單元開始，聯(lián)配可延三種可能的方式延伸：如果堿基不匹配，那么每一序列加上一個堿基，并給其增加一個規(guī)定的距離權(quán)重；或在一個序列中增加一個堿基而在另一序列中增加一個空位或反之亦然。引入一個空位時也將增加一個規(guī)定的距離權(quán)重。因此，表中的一個單元可以從(至多)三個相鄰的單元到達(dá)。我們把達(dá)左上角單元距離最小的方向看作相似序列延伸的方向。等距離時意味著存在兩種可能的方向。將這些方向記錄下來，并在研究了所有的單元之后，沿著記錄的方向就有一條路徑可從右下角(兩個序列的末端)追蹤到左上角(兩個序列的起點(diǎn))。由此所產(chǎn)生的路徑將給出具有最短距離的序列聯(lián)配。以兩個短序列CTGTATC和CTATAATCCC為例，將上述過程說明于圖3.4。設(shè)堿基錯配時距離權(quán)重為1，引入一個空位時距離權(quán)重為3。該圖邊緣的行和列作為起始條件增加到表中。在單元5行3列，即相應(yīng)較短序列(第二序列)的第2個T堿基和較長序列(第一序列)的第1個T堿基位置，有三種可能的距離增量。設(shè)在各序列中增加堿基T時(從4行2列移動)對距離的奉獻(xiàn)為0。從5行2列的位置作水平移動(等價于增加第二序列的堿基T而在第一序列引入一個空位)，在本例中增加一個罰值3。從3列4行向該單元作垂直移動，使第一序列增加堿基T而第二序列引入一個空位，結(jié)果也得到一個罰值3。因此從該單元(5行3列)所得到的最小距離的延伸方向是沿對角線和水平方向。在表中這兩個方向用箭頭表示。這兩種最短方向都使從左上角到該單元的距離為6。在上述6種聯(lián)配中，距離均為10，即在較短序列中有6個匹配堿基、1個錯配堿基和3個空位。沿箭頭所指方向在表中從右下角向左上角追蹤，得到6種可能的聯(lián)配：該算法可以用代數(shù)形式來描述。設(shè)具有堿基ai和bj的兩個序列a和b，這兩個序列間距離為d(a,b)。通過評價序列a中前i個位置和序列b前j位置的距離，遞歸地得到距離d(ai,bj)。如果a和b的長度為m和n，那么其期望距離為d(am,bn)。上表中引入的第1行1列單元的距離為0(相當(dāng)于空序列)，在單元(i,j)內(nèi)，使到達(dá)該單元距離增加的三種可能事件為：1.從單元(i-1,j)向(i,j)的垂直移動，相當(dāng)于在b序列中插入一個空位使相似序列延伸。換言之，b序列由a序列中ai的缺失所產(chǎn)生，這一事件的權(quán)重記作w_(ai)。2.從單元(i-1,j-1)向(i,j)的對角線移動，相當(dāng)于增加堿基ai和bj使相似序列延伸。換言之，b序列由a序列中的ai被bj取代所產(chǎn)生，這一事件的權(quán)重記為w_(ai,bj)。3.從單元(i,j-1)向(i,j)的水平移動，相當(dāng)于在序列b中插入一個空位使相似序列延伸。換言之，b序列由bj插入a序列所產(chǎn)生，這一事件的權(quán)重記為w+(bj)。因此，單元(i,j)的距離可看成三個相鄰單元的距離加上相應(yīng)權(quán)重后的最小者，即方程3.3且初始條件為在圖3.4的實例中當(dāng)兩個序列被聯(lián)配時，通過計算其重排序列(shuffedversion)的聯(lián)配距離，可以得到這兩個序列間的最小距離估計。如果實際得到的聯(lián)配距離小于重排序列距離的95%，那么說明實際的聯(lián)配距離到達(dá)了5%的顯著水平，是不可能由機(jī)誤造成的。二、Smith-Waterman算法由于親緣關(guān)系較遠(yuǎn)的蛋白質(zhì)序列可能只有一些相互獨(dú)立的相同片段，所以進(jìn)行局部相似性分析有時可能比整體相似性分析更合理。Smith和Waterman描述了一種查找具有最高相似性片段的算法。對于序列A=(a1,a2,…,am)和B=(b1,b2,…,bn)，Hij被定義為以ai和bj堿基對結(jié)束的片段(亞序列)的相似性值。與Needle-Wunsch算法一樣，Smith-Waterman算法也要利用遞推關(guān)系來確定H值，H的初始值為：相似性計算中包括2個統(tǒng)計量：堿基對(序列因子)ai,bj的相似性值S(ai,bj)和空位權(quán)重wk=v+uk(k為空位長度)。Smith-Waterman算法可以給出2條序列的最大相似性值。以ai,bj堿基對結(jié)束的片段可以由以ai-1和bj-1結(jié)束片段增加堿基(因子)來獲得，或者ai可以刪除k長度的堿基片段，bj可刪除l長度堿基片段。具體算法如下：方程3.4那么方程3.5其中該算法可以確保具有最大Hij值的序列片段是相似性最好的。從(ai,bj)為起點(diǎn)，向后追蹤矩陣，直到到達(dá)某一負(fù)值。對于具有最大相似性片段以外部分的差異性不會影響到該片段的H值。舉例說明了這一算法。我們同樣以上節(jié)Needleman-Wunsch算法中的兩條短序列為例。兩條序列(CTGTATC和CTATAATCCC)排于表3.9的兩側(cè)，相應(yīng)的和值分別列入表中。本例的權(quán)重等根據(jù)Smith和Waterman(1981)以前的例子設(shè)定為：方程3.6對于4個堿基具有相同頻率的隨機(jī)長序列，S(ai,bj)值的平均值為零。wk值應(yīng)至少不小于匹配與不匹配權(quán)重的差值。表3.9的最大Hij為4.33(8行與7列相交處)，星號(*)表示出具有最大相似性的片段匹配方式：三、序列相似性統(tǒng)計特征到目前為止，對局部聯(lián)配的統(tǒng)計學(xué)問題已根本搞清楚，特別是那些不含有空位(gap)的局部聯(lián)配更是如此。我們不妨首先考慮不含有空位的局部聯(lián)配問題，BLAST最初的搜索程序便是以此為根底的。無空位局部聯(lián)配涉及的是等長度的一對序列