歡迎關(guān)注：oddxix

本章主要講序列分析與聯(lián)配

序列分析是生物信息學(xué)最主要的研究內(nèi)容之一，它可以分為兩個(gè)主要部分：

一、是序列組成（特別是涉及到基因組層次上）分析。

二、是序列之間的比較分析。

兩條序列或多條序列間的比對(duì)或聯(lián)配(alignment)的目的，是對(duì)它們的序列相似性進(jìn)行評(píng)估，找出這些序列中結(jié)構(gòu)或功能相似性區(qū)域等。通過聯(lián)配未知序列與已知序列(其功能或結(jié)構(gòu)等已知)的相似程度，我們可以判斷或推測(cè)未知序列的結(jié)構(gòu)與功能。

序列組成及單一序列分析

在基因每一側(cè)的500 個(gè)任意堿基區(qū)域被稱為“側(cè)翼”，基因間區(qū)域是指兩個(gè)基因間的其余序列，DNA序列內(nèi)和序列間堿基具有不同的頻率。

分析 DNA 序列的主要困難之一是堿基相鄰的頻率不是獨(dú)立的。堿基相鄰的頻率一般不等于單個(gè)堿基頻率的乘積：如果 P u 是序列中堿基 u 的頻率，且 P uv 為兩個(gè)相鄰堿基 u 和 v 的頻率，則

Puv≠Pu*Pv

在編碼區(qū)，存在某種約束來限制 DNA 序列編碼氨基酸。在密碼子水平上，這一約束與堿基相鄰頻率有關(guān)。編碼同一氨基酸的不同密碼子(同義密碼子)好像不是等同存在的。這種密碼子偏倚必定與兩堿基相鄰頻率水平有關(guān)。

相鄰堿基之間的關(guān)聯(lián)將導(dǎo)致更遠(yuǎn)堿基之間的關(guān)聯(lián)，這些關(guān)聯(lián)延伸距離的估計(jì)可以從馬爾科夫鏈(Markovchain)理論得到。在不援引任何生物學(xué)機(jī)制的情況下，第 k 階馬爾科夫鏈假定在序列中某一位置上堿基的存在只取決于前面 k 個(gè)位置上的堿基。一階鏈假定一個(gè)特定堿基存在于位置 i的概率只取取決于在位置 i-1 的 4 種堿基概率。相互獨(dú)立的堿基所組成的序列將與 0 階馬爾科夫鏈相對(duì)應(yīng)。階可以通過似然法估計(jì)。同時(shí)，馬爾科夫鏈分析更適應(yīng)于基因組水平，而非單一序列(基因)。

Needleman-Wunsch 算法是一種整體聯(lián)配(global alignment)算法，最佳聯(lián)配中包括了全部的最短匹配序列。Smith-Wateman算法是在 Needleman-Wunsch 算法基礎(chǔ)上發(fā)展而來的，它是一種局部聯(lián)配(Local alignment)算法。這二種算法均可以用于核酸和蛋白質(zhì)序列。在給定空位罰值和替換矩陣情況下，它們總是能給出具有最高（優(yōu)）聯(lián)配值的聯(lián)配。

從整體上分析兩個(gè)序列的關(guān)系，即考慮序列總長的整體比較，用類似于使整體相似(global similarity)最大化的方式，對(duì)序列進(jìn)行聯(lián)配。兩個(gè)不等長度序列的聯(lián)配分析必需考慮在一個(gè)序列中圈掉一些堿基或在另一序列作空位(gap)處理。

Needleman-Wunsch 算法

Needleman-Wunsch 算法是為氨基酸序列發(fā)展的，但也可以用于核苷酸序列。算法最初尋求的是使兩條序列間的距離最小。盡管這類距離的元素是以一種特定的方式定義的，但該算法的良好特性在于它確定了最短距離。這是一個(gè)動(dòng)態(tài)規(guī)劃(dynamic programming)的方法。

將兩條聯(lián)配的序列沿雙向表的軸放置，兩條序列的所有可能的聯(lián)配方式都將在它們所形成的方形圖中。從任一堿基對(duì)，即表中的任一單元開始，聯(lián)配可延三種可能的方式延伸：如果堿基不匹配，則每一序列加上一個(gè)堿基，并給其增加一個(gè)規(guī)定的距離權(quán)重；或在一個(gè)序列中增加一個(gè)堿基而在另一序列中增加一個(gè)空位或反之亦然。引入一個(gè)空位時(shí)也將增加一個(gè)規(guī)定的距離權(quán)重。因此，表中的一個(gè)單元可以從(至多)三個(gè)相鄰的單元達(dá)到。我們把達(dá)左上角單元距離最小的方向看作相似序列延伸的方向。等距離時(shí)意味著存在兩種可能的方向。將這些方向記錄下來，并在研究了所有的單元之后，沿著記錄的方向就有一條路徑可從右下角(兩個(gè)序列的末端)追蹤到左上角 (兩個(gè)序列的起點(diǎn))。由此所產(chǎn)生的路徑將給出具有最短距離的序列聯(lián)配。

Smith-Waterman 算法

由于親緣關(guān)系較遠(yuǎn)的蛋白質(zhì)序列可能只有一些相互獨(dú)立的相同片段，所以進(jìn)行局部相似性分析有時(shí)可能比整體相似性分析更合理。對(duì) 于序 列 A=(a 1 ,a 2 , … ,a m ) 和B=(b 1 ,b 2 , … ,b n )，H ij 被定義為以a i 和b j 堿基對(duì)結(jié)束的片段(亞序列)的相似性值。與Needle-Wunsch算法一樣，Smith-Waterman算法也要利用遞推關(guān)系來確定H值。

相似性計(jì)算中包括 2 個(gè)統(tǒng)計(jì)量：堿基對(duì)(序列因子) 的相似性值和空位權(quán)重 (k 為空位長度)。Smith-Waterman算法可以給出2 條序列的最大相似性值。以 堿基對(duì)結(jié)束的片段可以由以 和 結(jié)束片段增加堿基(因子)來獲得，或者 可以刪除 k 長度的堿基片段，可刪除 l 長度堿基片段。

該算法可以確保具有最大H ij 值的序列片段是相似性最好的。從 為起點(diǎn)，向后追蹤矩陣，直到到達(dá)某一負(fù)值。對(duì)于具有最大相似性片段以外部分的差異性不會(huì)影響到該片段的H值。

序列相似性的統(tǒng)計(jì)特性

先考慮不含有空位的局部聯(lián)配問題

無空位局部聯(lián)配涉及的是等長度的一對(duì)序列片段，兩個(gè)片段的各部分彼此比較。一種Smith-Waterman 或 Sellers 算法的改進(jìn)算法可以找到所有高比值片段對(duì)(high-scoring segment pairs,HSPs)，即這些片段對(duì)的比較分值不會(huì)因片段的延伸而進(jìn)一步升高。

為了分析上述分值隨機(jī)性產(chǎn)生的幾率大小，需要建立一個(gè)隨機(jī)序列模型。對(duì)于蛋白質(zhì)而言，最簡(jiǎn)單的序列模型可通過從一條序列中隨機(jī)地選取氨基酸殘基，當(dāng)然這一條序列中各種殘基的頻率必需一定。另外，一對(duì)隨機(jī)氨基酸的聯(lián)配期望值必需為負(fù)值，否則不論聯(lián)配片段是否相關(guān)的，都會(huì)得到高比值。

在一定的序列長度 m 和 n 限定下，HSP 的統(tǒng)計(jì)值可由 2 個(gè)參數(shù)(k 和λ)確定。最簡(jiǎn)單的形式，即不小于比較值為 S的 HSP 個(gè)數(shù)，可由下列公式算得其期望值：

在給定比值的情況下，將比較序列長度加倍，則 HSP數(shù)(即 E 值)也將加倍，同樣，S 值為 2X 的某個(gè) HSP 長度必是 S 值為 X 的兩倍，所以 E 值將隨著 s 值的增大急劇減少。參數(shù) K 和λ可分別被簡(jiǎn)單地視為搜索步長(search spacesize)和計(jì)分系統(tǒng)(scoring system)的特征數(shù)。

P 值(P-Value)(概率值)

具有大于或等于某一比值 S 的隨機(jī)HSP 數(shù)可由泊松分布(Poisson distribution)確定。由此可以計(jì)算出搜索到某一比值大于或等于 S 的 HSP 的機(jī)率為

空位聯(lián)配(gapped alignment)的統(tǒng)計(jì)問題

對(duì)于非空位聯(lián)配，可用基于替換矩陣和比較序列的殘基頻率的辦法估計(jì)統(tǒng)計(jì)參數(shù)；對(duì)于空位聯(lián)配，參數(shù)的估計(jì)則必須根據(jù)“隨機(jī)”序列的大尺度比較。

空位罰值(gappenalties)

聯(lián)配中另一個(gè)重要問題是空位問題?？瘴惶幚硎轻槍?duì)序列進(jìn)化過程中可能發(fā)生的插入和缺失而設(shè)計(jì)的。插入和缺失可能只涉及 1 個(gè)或 2 個(gè)殘基，也可能是整個(gè)功能域(domain)，所以，在進(jìn)行空位罰值設(shè)計(jì)時(shí)必須反映這些情況。有 2 個(gè)參數(shù)應(yīng)用于空位罰值設(shè)定，一個(gè)與空位設(shè)置(gap opening)有關(guān)，另一個(gè)與空位擴(kuò)展(gap extension)有關(guān)。任一空位的出現(xiàn)均處以空位設(shè)置罰值，而任一空位的擴(kuò)大必須處于空位擴(kuò)展罰值。

對(duì)于一個(gè)空位長度為k的罰值W K 可用下式表示：

其中 a 是空位設(shè)置罰值，b 為空位擴(kuò)展罰值。這兩個(gè)參數(shù)值設(shè)置的變化對(duì)聯(lián)配產(chǎn)生影響

替換矩陣的一般原理

我們并不能直接計(jì)算出兩條序列的最佳聯(lián)配，我們需要找到一個(gè)可以估計(jì)任何聯(lián)配的某一統(tǒng)計(jì)數(shù)，使生物學(xué)關(guān)系匹配最顯著的聯(lián)配統(tǒng)計(jì)數(shù)最大。替換矩陣(substitution matrices)包括了在聯(lián)配中各種匹配方式如何賦分的信息，故替換矩陣又常被稱為計(jì)分矩陣(scoring matrices)。

用于蛋白質(zhì)聯(lián)配的替換矩陣要復(fù)雜一些，因?yàn)闆]有一個(gè)矩陣可以適用各種情況。構(gòu)建矩陣時(shí)應(yīng)考慮不同的蛋白質(zhì)家族在進(jìn)化過程中，一種氨基酸突變成另一種氨基酸概率的差異，根據(jù)不同的蛋白質(zhì)家族和預(yù)期的相似程度構(gòu)建不同的替換矩陣。2 個(gè)最有名的蛋白質(zhì)替換矩陣是 PAM和 BLOSUM

PAM 氨基酸替換矩陣

在進(jìn)行蛋白質(zhì)序列聯(lián)配時(shí)，必須通過一定的方法給聯(lián)配的殘基對(duì)賦予一定的分值，替換矩陣便是其中最重要的方法。已故 Dayhoff 是蛋白質(zhì)列序比較的先驅(qū)，她和她的同事們通過對(duì)蛋白質(zhì)進(jìn)化模式的研究，建立了一組被廣泛應(yīng)用的替換矩陣，這些矩陣常被稱為Dayhoff，MDM(Mutation Data Matrix)或 PAM(Percent Accepted Mutation)矩陣。

由于蛋白質(zhì)最有可能是自然選擇的目標(biāo)，可以認(rèn)為蛋白質(zhì)序列的分析比DNA分析更具有生物學(xué)意義。蛋白質(zhì)分析完全避免了幾個(gè)三聯(lián)體可能編碼同一氨基酸的遺傳密碼簡(jiǎn)并問題。有必要進(jìn)一步分析各種氨基酸間的同源性程度，以及在進(jìn)化過程中一種氨基酸被另一種氨基酸替換的概率大小。也許把氨基酸按一定特性分成若干組更便于以上分析，例如氨基酸可分成中性疏水(G、A、V、L、I、F、P、M)、中性親水(S、T、Y、W、N、E、C)、堿性(K、R、H)和酸性(D、E)氨基酸等。在比較許多具有相似性蛋白質(zhì)序列的基礎(chǔ)上，Dayhoff等于 1979 年構(gòu)建了一個(gè)突變概率矩陣M(mutation probability matrix)。

BLOSUM 氨基酸替換矩陣

另外一種構(gòu)建矩陣的方法是由 Henikoff 等于 1992 年提出的，建成的矩陣為BLOSUM(Blocks SubstitutionMatrices)。他們直接利用多序列聯(lián)配(multiple alignment)分析親緣關(guān)系較遠(yuǎn)的蛋白質(zhì)，而不是用相近的序列。這方法的優(yōu)點(diǎn)是符合實(shí)際觀測(cè)結(jié)果，不足之處是它不能和進(jìn)化掛起鉤來。大量的試驗(yàn)表明，BLOSUM矩陣總體比 PAM 矩陣更適合于生物學(xué)關(guān)系的分析和局部相似性搜索。

DNA 替換矩陣

以上有關(guān)替換矩陣的討論僅僅提及蛋白質(zhì)序列的比較，但是，相關(guān)的原則同樣適用于 DNA 序列的比較。在進(jìn)行比較時(shí)應(yīng)該意識(shí)到，用翻譯而來的蛋白質(zhì)序列總是好于直接用 DNA 序列。這是因?yàn)?DNA 序列的進(jìn)化變化很少，在使用簡(jiǎn)單的DNA 替換矩陣比較時(shí)，獲得的同源性信息遠(yuǎn)少于蛋白質(zhì)序列。DNA 替換矩陣非常簡(jiǎn)單，所有 4 個(gè)堿基的匹配與不匹配的數(shù)值均設(shè)為相同，不同的只有匹配與否(0.9 和-0.1)。一個(gè)較復(fù)雜的模型是把轉(zhuǎn)換(transition,兩種嘧啶或兩種嘌呤間的突變)頻率設(shè)為高于顛換(tranversion,嘧啶與嘌呤間的突變)頻率。

多序列聯(lián)配

通過以上的兩條序列算法，總是可以返回一個(gè)最佳匹配的聯(lián)配結(jié)果。但是，當(dāng)我們將兩條以上的序列放在一起聯(lián)配時(shí)，情況就就不一樣了。現(xiàn)有實(shí)用的多序列聯(lián)配方法還不能保證一定給出最優(yōu)聯(lián)配結(jié)果，只能給出一個(gè)近似值——往往人為的修正可以使聯(lián)配結(jié)果更佳。

同源序列的多序列聯(lián)配是生物信息學(xué)一個(gè)重要課題。通過多序列聯(lián)配結(jié)果，允許你觀察殘基可以改變到什么程度而蛋白質(zhì)仍保持功能；它也可以使你得到圍繞某一殘基的三級(jí)結(jié)構(gòu)信息。有不少多序列聯(lián)配程序可通過匿名 ftp 等服務(wù)獲得，例如：ClustalW等。三條或三條以上序列的聯(lián)配方法可分為幾類，如用于兩條序列聯(lián)配的Needleman-Munsch 等算法的改進(jìn)算法、等級(jí)法(hierachical method)、片段法(segment method)、一致或區(qū)段法(consensus or regions'method)等。這些方法中，等級(jí)法是目前應(yīng)用最為廣泛的方法。

等級(jí)法又稱為樹法(tree method)，是由 Feng 和 Doolittge(1987)等人發(fā)展的（ClustalW 程序）。由于兩條序列的聯(lián)配結(jié)果可以很容易地獲得，多序列聯(lián)配便可以在連續(xù)使用兩條序列聯(lián)配算法(如 Needleman-Wunsch 算法)基礎(chǔ)上，通過先建“樹”的思路來進(jìn)行多序列聯(lián)配。這一方法同樣是一種動(dòng)態(tài)規(guī)劃方法。具體步驟如下：

• ①對(duì)所有序列進(jìn)行兩兩聯(lián)配分析，N 條序列應(yīng)有 N×(N-1)/2 對(duì)；
• ②對(duì)兩兩聯(lián)配的數(shù)據(jù)進(jìn)行聚類分析，產(chǎn)生聯(lián)配等級(jí)。該等級(jí)可用分叉樹(binary tree)形式或簡(jiǎn)單的排序來表示；
• ③根據(jù)以上聯(lián)配結(jié)果，首先從所有聯(lián)配中相似性最好的兩條序列開始，然后是剩余聯(lián)配中相似性最好的兩條序列??依次類推，直至多序列聯(lián)配結(jié)束。一旦兩條序列的聯(lián)配被列入，則序列的位置就被固定下來。例如，對(duì)于序列 A、B、C、D，如果 A 與 C、B 與 D 分別是兩兩聯(lián)配的最佳聯(lián)配結(jié)果，則 A、B、C、D 四條序列的聯(lián)配則通過比對(duì) A-C 和 B-D 兩個(gè)聯(lián)配(每個(gè)聯(lián)配位置取平均值)來確定。這一組合方法對(duì)大量序列的多序列聯(lián)配提供了實(shí)用的空位聯(lián)配手段，除了最初的兩序列間的聯(lián)配過程，整個(gè)多序列聯(lián)配過程是很快的。

數(shù)據(jù)庫搜索——BLAST 和 FASTA 應(yīng)用

比較和確定某一數(shù)據(jù)庫中的序列與某一給定序列的相似性是生物信息學(xué)中最頻繁使用和最有價(jià)值的操作。本質(zhì)上這與兩條序列的比較沒有什么兩樣，只是要重復(fù)成千上萬次。但是要嚴(yán)格地進(jìn)行一次比較必定需要一定的耗時(shí)，所以必需考慮在一個(gè)合理的時(shí)間內(nèi)完成搜索比較操作。目前有二個(gè)最為常用的程序服務(wù)于未知序列的數(shù)據(jù)庫相似性搜索，即 BLAST 和 FASTA。FASTA 使用的是Wilbur-Lipman 算法的改進(jìn)算法，進(jìn)行整體聯(lián)配，重點(diǎn)查找那些可能達(dá)到匹配顯著的聯(lián)配。雖然 FASTA 不會(huì)錯(cuò)過那些匹配極好的序列，但有時(shí)會(huì)漏過一些匹配程度不高但達(dá)顯著水平的序列。BLAST(Basic Local Alignment Search Tool，基本局部聯(lián)配搜索工具)是基于匹配短序列片段，用一種強(qiáng)有力的統(tǒng)計(jì)模型來確定未知序列與數(shù)據(jù)庫序列的最佳局部聯(lián)配。

大多數(shù)研究目前都通過國際互聯(lián)網(wǎng) Internet 應(yīng)用 NCBI 研制的 BLAST 程序(Basic Local Alignment Search Tool)來進(jìn)行 DNA 和蛋白質(zhì)序列相似性搜索。用一組 BLAST 程序聯(lián)配可以快速進(jìn)行核酸和蛋白質(zhì)序列庫的相似性檢索。
采用BLAST 的基本算法編成了若干各不同的程序，分別使用特定的序列庫和用于特定類型的輸入序列。BLASTN 是在核苷酸序列庫搜索核苷酸序列。BLASTP 是在蛋白質(zhì)序列庫中搜索氨基酸序列。TBLASTN 則可以在核酸序列庫中搜索氨基酸序列，此時(shí)序列庫在搜索之前要按所有 6 種讀框即時(shí)翻譯。與此相反的一項(xiàng)分析則由BLASTX 來完成，它要將所輸入的核酸序列按所有 6 種讀框翻譯，然后再以之搜索蛋白質(zhì)序列庫。PSI-BLAST 可以對(duì)數(shù)據(jù)庫進(jìn)行多輪循環(huán)檢索，每一輪的檢索速度都大約是 BLAST 的兩倍，但每一輪都能提高檢索的敏感性。它是目前 BLAST 程序家族中敏感性性最高的成員。

如果目的序列中有蛋白質(zhì)編碼區(qū)，則用翻譯的蛋白質(zhì)序列來搜索蛋白質(zhì)序列庫要比用 DNA 序列搜索核酸序列庫更有價(jià)值。由于蛋白質(zhì)序列的進(jìn)化要比 DNA序列慢一些，在蛋白質(zhì)序列水平上的遠(yuǎn)緣關(guān)系在 DNA 水平上可能被錯(cuò)過。如果無法確定編碼區(qū)，則可利用 BLASTX 按所有 6 種讀框來翻譯 DNA 序列，然后用它搜索蛋白質(zhì)序列庫。由于蛋白質(zhì)序列庫僅包含已鑒定的蛋白質(zhì)，所以必須采用TBLASTN 程序在現(xiàn)有的 GenBank、EMBL 或 DDBJ DNA 序列庫中檢索新確定的氨基酸或翻譯過的 DNA 序列。這種檢索有時(shí)可以找到一些顯著相似的 DNA 序列，而原本并不知道這些序列可編碼蛋白質(zhì)。

另一個(gè)常用的核酸和蛋白質(zhì)序列庫搜索程序是 FASTA，即 FASTN 和 FASTP 程序的新版本。FASTA 首先在序列庫中進(jìn)行快速的初檢，找出與待檢序列高度相似的序列。這一快速檢索局限于待檢序列和序列庫序列之間較短的完全相同序列區(qū)段上。

無論采用 FASTA 或 BLAST，推斷相似性是否具有生物學(xué)意義都取決于研究者。要作出決斷，必須充分考慮蛋白質(zhì)已知的或推斷的功能，與已知活性位點(diǎn)或模序的相似程度等等。因?yàn)锽LAST和FASTA采用不同的算法，同時(shí)用這兩種搜索引擎重新檢索某一特定序列往往是可取的。如果用其中一種找不到顯著相似序列，不妨試一試另一程序。如果BLAST和FASTA均找不到顯著匹配的序列，還可以選擇第 3 條比較費(fèi)時(shí)的搜索策略。一些網(wǎng)站允許用戶使用基于Smith-Waterman算法的搜索程序，如BLITZ。BLITZ( www.ebi.ac.uk/searchs/blitz.html) 被設(shè)計(jì)在大型并行計(jì)算機(jī)上運(yùn)行，因此使檢索更靈敏。雖然運(yùn)行這樣的程序比較費(fèi)時(shí)，但它們有時(shí)會(huì)發(fā)現(xiàn)一些被BLAST和FASTA錯(cuò)過的勉強(qiáng)達(dá)到顯著的聯(lián)配。