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

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

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