基于置換距離度量的蛋白質多序列比對算法性能評估*

高 峰， 李防震， 王 珺，董騮煥

(1.山東經濟學院計算機科學與技術學院∥山東省數字媒體技術重點實驗室，山東 濟南 250014；2.中國科學院-馬普學會計算生物學伙伴研究所，上海 200031)

蛋白質多序列比對是蛋白質組學研究中的一項重要而具有基礎性的研究課題，在整個生物信息學和系統生物學研究中亦是至關重要的一步。比對算法的好壞直接影響到其它生物信息學研究結果的準確性，低質量的比對算法甚至會給相關研究帶來導致錯誤結論的危險。因此，研究對蛋白質多序列比對算法的評估方法就顯得尤為重要。

目前，蛋白質多序列比對算法的打分方法有很多[1-2]，比如：WSP(weighted sum-of-pairs)，最大似然，最小熵等等。現在最流行的打分方法是WSP，但是在使用CluatalW，Tcoffee等方法的比對測試中，63%的比對都比正確的比對有較高的WSP值，也就是說WSP分數與比對質量之間的對應是比較差的。后來Thompson等[3]引入兩種計分來比較多序列比對：列分數和SPS (sum-of-pairs score)。

列分數是計算參考比對與待檢驗比對序列中相同列的個數。這種衡量反映了程序正確比對所有序列的能力，所以任何一條序列的比對錯誤既可導致比對得分為零。SPS計算待檢驗比對與參考比對間相同的殘基對的比例。后來，Karplus等[4]進一步修改了這種打分方法，給它賦予了權值。Lassmann等則對SPS方法做了另外的修改，提出了交疊分數(overlap score)法。該方法把參考比對中的所有殘基對及待檢驗比對中的所有殘基對分別存入兩個集合中，然后計算這兩個集合的交與并的比值。所有這些評估算法都基于序列比對結果的局部性質，因而都存在著易受數據噪聲影響的缺陷。

通常情況下，研究者得到的DNA和氨基酸序列數據總是含有噪聲的。噪聲主要來源于測序過程，目前的DNA和氨基酸測序技術還難以達到100%的準確度；另外，測序后的序列數據在存儲和傳輸過程中也有可能引入誤差，例如，由于操作人員的粗心導致的字母輸入錯誤。因此，對序列數據的處理程序，需要有一定的容錯能力和抗噪聲能力。

本文提出了一種蛋白質多序列比對算法的性能評估新方法，即置換距離法。由于置換距離法僅關注于不同蛋白質-蛋白質進化距離之間的相對次序，而不考慮這些進化距離之間的細微差異，因而具有更強的魯棒性。它能夠克服氨基酸序列數據中噪聲的影響，對多種蛋白質序列比對算法做出客觀公正的評價。

另外用最長公共子序法度量置換距離，能夠比較準確的反映出置換距離之間的差異性。基于該算法，我們對Dialign， Tcoffee， ClustalW和Muscle多序列比對算法進行了比較評估。

1 材料與方法

1.1 數據源

1.1.1 BALIBASE數據庫數據 BALIBASE(Benchmark Alignment Database)是針對多序列比對問題而設計的數據資源，它提供的所有多序列比對都是來源于三維結構的疊加，是目前用作多序列比對結果比較的標準平臺，我們應用BALIBASE作為我們的數據源。BALIBASE包含了八個等級的參考系列，由于其中的RV6-RV8還不成熟，所以我們只應用了RV1-RV5共142組多序列進行測試。

1.1.2 ROSE軟件生成的數據 ROSE軟件可以生成蛋白質序列進化的概率模型。通過進化樹，可以經過插入、刪除和替換，從同一祖先生成具有一定相似性的序列。在人工的進化過程中，我們可以指定序列長度，進化距離，進化祖先等多項參數，來產生具有一定生物意義的序列。ROSE生成的序列適合于多序列比對，也廣泛應用于生物進化關系的預測。我們將其作為另一個數據源。

1.2 蛋白質多序列比對算法介紹

本文選取Dialign，Tcoffee，ClustalW和Muscle這四種目前比較流行的蛋白質多序列比對算法進行性能評估，從而驗證我們算法的有效性。

Dialign是最近幾年才發展起來的一個多序列比對算法[5]。它注重于尋找多條序列中相似的區域，是一種基于段的漸進比對算法。它對于非全局相關，只是局部相似的序列非常有效。在序列比對軟件當中，它的準確性是非常高的，但是運行時間很長。而且輸出格式不夠靈活，使得比對比較耗時，并需要進一步處理比對后的結果。

Tcoffee是一個多序列比對包[6]。給定一組蛋白質或者DNA序列，它便可以生成多序列比對。Tcoffee可以連接多序列或者一對序列，也可以連接全局比對序列或者局部比對序列到一條序列中。不論數據源來自哪里，它都能通過剩余的部分來判斷得出每個位置與新序列聯系的緊密性。這種緊密型通常是序列比對準確性的指示器。另外，它也提供多種格式的輸出結果。

Clustal W可以進行蛋白質與核酸的多序列比較[7]，分析不同序列之間的相似性關系，還可以繪制進化樹。由于其靈活的輸入輸出格式、方便的參數設定和選擇以及良好的可移植性，使得ClustalW在蛋白質與核酸的序列分析中得到了廣泛應用。它可以用來發現特征序列，進行蛋白分類，證明序列間的同源性，幫助預測新序列二級結構與三級結構，確定PCR引物，以及用于分子進化分析。

Muscle在比對過程中采取了兩次迭代過程來提高比對的精度[8]。對于一對序列Muscle提供了兩種距離測量方法：一種針對未比對序列的kmer距離和一種針對已比對序列的Kimura距離。Muscle算法會首先進行一次漸進式比對，然后利用迭代精細法對多序列比對結果進行進一步的優化。Muscle采用的是基于進化樹的分組迭代法。

1.3 置換距離法介紹

設Ω為n個自然數構成的集合，即Ω={1,2,3,…,n}。集合Ω到其自身上的一一映射叫做Ω上的一個置換[9]。集合Ω上的所有置換構成置換群，共有n!個元素。

例1： 集合Ω={1,2,3}共有3!=6個置換，分別為{1,2,3}，{1,3,2}，{2,1,3}，{2,3,1}，{3,1,2}，{3,2,1}。 其中{1,2,3}稱為恒等置換。

本文中置換的生成過程如下：從數據庫中提取一組序列數據，設含有m個氨基酸序列。由于數據庫自身已給出各組序列的參考比對，我們可基于這m個序列的參考比對生成任意一對序列之間的進化距離。顯然，m條序列共有n=m(m-1)/2個兩兩組合(稱為二元組)，構成一個n維的距離向量。我們把這n個距離由小到大進行排序，并分別給予編號1～n，得到n個自然數(構成集合Ω)的恒等置換π0={1,2,...,n}，同時得到這m條序列的n維二元組集合到自然數集合Ω上的一一映射。

選取一個待考察的比對算法，如Dialign, 對這組序列重新進行比對。然后按與前面相同的過程，基于此比對結果生成n維的進化距離向量，并把這n個距離由小到大進行排序。需要注意的是，這里得到的n個二元組的排列順序與前面是不一樣的。我們把這n個二元組按與前面相同的映射關系映射到集合Ω中，即得到集合Ω的對應于該比對算法的置換π1={π(1),π(2),...,π(n)}。

因此，對于每個比對算法，我們都可得到n個自然數集合Ω關于這個算法的置換。

例2：假設從數據庫中提取一組氨基酸序列，含有3條序列a，b和c。3條序列共有3×(3-1)/2=3個兩兩組合，其進化距離分別記為D(a,b),D(b,c),D(c,a)。

按照數據庫給出的參考比對得到的進化距離排序設為：D(a,b)

接著使用某個比對算法得到的距離排序設為：D(c,a)

這樣，通過進一步測量不同置換之間距離的大小，就提供了一種評價各種多序列比對算法性能的有效方法。采用該方法，數據的細節差異性被排序后的數字所隱藏了，我們關心的已不再是原始的數據，而是置換后的數字串。因而，置換距離法僅關注于不同蛋白質-蛋白質進化距離之間的相對次序，而不考慮這些進化距離之間的細微差異，因而具有更強的魯棒性。

1.4 置換距離的最長公共子序度量

關于置換的距離問題，人們已發展了多種度量方法，我們采用最長公共子序(Longest Common Subsequence, LCS)法度量置換之間的距離。子序即子序列，是與原字符串具有相同出現順序的字符的子集。本文所處理的字符串是由自然數組成的數字串。

例3：[3,5,7]是[1,2,3,4,5,6,7]的子序。

兩條串所共有的最長子序的長度可用來衡量兩條串的相似性[10]。設原串長度都為n，則當兩條串完全相同時，其最長公共子序長度達到最大值n；當兩條串字符順序完全相反時，其最長公共子序長度達到最小值1。我們定義最長公共子序距離為n減去最長公共子序的長度，用dlcs表示，那么dlcs的范圍就是：[0,n-1]。

規范化的最長公共子序距離為最長公共子序與n-1的比值，即

(1)

最長公共子序可以通過一種動態規劃方法實現[11]，算法復雜度為O(n2)。首先介紹LCS相似度[12]：字符串a和b的LCS相似度是a和b間的最大相同子序的長度，記為LCS(a,b)。顯然LCS(a,b)越大，a，b越相似。LCS(a,b)的動態規劃計算公式為

(2)

其中i和j分別為字符串a和b中字符位置的標記，LCS(i,j)表示兩條子串[a1,...,ai]和[b1,...bj]的最長公共子序，則LCS(a,b)就是字符串a和b的最長公共子序。然后，基于LCS(a,b)可計算得兩條串a，b間的最長公共子序距離dlcs，為

dlcs(a,b)=n-LCS(a,b)

(3)

使用最長公共子序的方法來比較兩個置換，比較靈敏準確。由于它面對的處理數據是數字串，而傳統的比較方法面對的處理數據是比對過的氨基酸序列。相對來說，最長公共子序所對應的數據量比較小，因而可以達到比較快的速度。

1.5 處理過程

1.5.1 序列比對 BALIBASE數據庫數據和ROSE生成的數據均是fasta格式的序列，這是蛋白質序列信息的通用存儲格式。從BALIBASE數據庫的物種類別中，我們得到了142組數據。使用ROSE軟件，我們又生成了3 721組序列，指定序列長度范圍：[50,650],步長為10;進化距離范圍：[0,300]，步長為5。所有序列數據都已給出了參考比對。對于兩部分數據，我們分別使用 Dialign，Tcoffee，ClustalW和Muscle 這4種多序列比對算法比對各組序列。從而每組序列數據都可以生成5種比對結果(包括參考比對)。比對后的序列，可以顯示出一些相似的區域以及差異較大的區域，反映了序列在進化上的保守區和非保守區。

1.5.2 計算距離矩陣 使用phylip 3.69程序包中的protdist.exe程序可以基于比對后的每組序列分別生成進化距離矩陣[13]。protdist.exe程序提供了多達5種基于比對序列計算進化距離的模型，本文采用其缺省設置Jones-Taylor-Thornton (JTT)模型[14]。對于含m條序列的分組(各組的m值可能是不同的)，得到的進化距離矩陣是m×m維的，每個元素對應一個序列對。顯然，這個矩陣是對稱的，且對角線元素為0，因而只有m(m-1)/2個獨立的非零元素，對應m條序列的m(m-1)/2個兩兩組合。這樣，每組序列都可以利用protdist.exe程序產生5個距離矩陣，分別為參考比對距離矩陣、Dialign距離矩陣、Tcoffee距離矩陣、ClustalW距離矩陣和Muscle距離矩陣。

1.5.3 排序并生成置換 每個m×m維的進化距離矩陣可提取出一個n=m(m-1)/2維的二元組向量，其元素為序列的兩兩組合。按前面所介紹的過程對每個二元組向量進行排序和映射，即得到每組序列數據在自然數集合Ω={1,2,...,n}的5種置換，分別記為：

參考比對產生的置換π0={1,2,...,n}，為恒等置換；

Dialign產生的置換πd={πd(1),πd(2),...,πd(n)}；

Tcoffee產生的置換πt={πt(1),πt(2),...,πt(n)}；

ClustalW產生的置換πc={πc(1),πc(2), ...,πc(n)}；

Muscle產生的置換πm={πm(1),πm(2), ...,πm(n)}。

2 試驗結果

2.1 基于BALIBASE數據庫的評價結果

圖1表示在BALIBASE數據庫的六類數據中，各比對算法表現出最優的比例。可以看出在所有6種類別的數據中，Tcoffee方法的表現都是最優的，Dialign次之，且這兩個算法的獲勝比例遠遠超過ClustalW和Muscle。在該圖中， Dialign和Tcoffee占了絕對優勢，Muscle和ClustalW的表現只局限于很小的數據范圍。

圖1 BALIBASE數據庫數據處理的最優比例圖

圖2為4種算法在各類別數據上的置換距離數值統計圖，注意圖中較小的數值對應著算法較好的性能表現。可以看出，Tcoffee和Dialign的數值相近，前者略優于后者，Muscle和ClustalW的表現差不多，Tcoffee和Dialign的數值顯著小于Muscle和ClustalW，這些結論與圖1也都是一致的。這些結果與其它評價方法所得結果是一致的[8]。

圖2 BALIBASE數據庫數據處理的各算法置換距離圖

2.2 基于ROSE數據的評價結果

我們使用ROSE軟件生成的數據進行實驗，結果見圖3和圖4。圖3為獲勝算法的分布示意圖，其中橫坐標為ROSE軟件的序列長度參數，縱坐標為進化距離參數。由圖3可以看出，在所有區間Dialign方法的最優比例都是最高的，Muscle，ClustalW和Tcoffee的比例都比較小。

圖3 基于ROSE數據的各算法獲勝分布示意圖

圖4為基于ROSE數據得到的置換距離統計圖。可以看出Dialign算法表現最好，顯著超過其它幾種算法。ClustalW算法、Muscle算法與Tcoffee算法表現差不多，它們的置換距離數值都處于較高的水平，這些結論與圖3也是一致的。

圖4 ROSE生成數據的各算法置換距離統計圖

3 討 論

在蛋白質多序列比對算法的性能評估方法的研究中，我們提出了置換距離法。置換距離法能夠克服氨基酸序列數據中噪聲的影響，對多種蛋白質多序列比對算法進行了客觀公正的評價；用最長公共子序法度量置換距離比較準確的反映出了不同置換之間的差異性。計算結果與其它序列比對評價算法所得結果是一致的[8]，另外，計算結果表明各算法的性能表現與處理的實際數據有關，不同算法具有不同的數據偏好。

我們分別用2個數據庫評價4種算法，得到的結論不盡一致。綜合來看，在BALIBASE數據庫上Dialign和Tcoffee表現較好，但Dialign又稍遜于Tcoffee；在ROSE數據上Dialign占絕對優勢，其它3種算法表現差不多都處于較低的數據范圍內。但需要指出的是，由圖2和圖4可以看出，實驗結果具有顯著的標準差。考慮這個因素之后，我們認為，總的來看，在這4個多序列比對算法中，Dialign表現出最佳的性能。

分析Dialign在BALIBASE上的表現稍差的原因，可能與Dialign屬于概率模型有關。BALIBASE數據來自真實的生物學實驗，而ROSE數據是根據概率模型人工產生的。因而，同樣基于概率模型的Dialign在ROSE數據中會表現更好一些。

本文采用了最長公共子序法來度量置換之間的距離，事實上，目前存在著多種置換距離度量方法。采用其他的置換距離度量方法，如編輯距離，或許也可以達到較好的效果。目前人們發展了大量蛋白質序列比對算法，本文僅選取了其中4種進行評估，計算結果初步證明了我們所采用的置換距離法的可行性。我們希望將來用置換距離法對目前存在著的各種序列比對算法做一個全面系統的性能評價研究，以考察不同比對算法所適用的數據范圍，從而對其他研究者在選擇序列比對算法時提供參考。另外，本文僅討論了蛋白質多序列比對算法的評估問題，但置換距離法顯然也適用于評估各種DNA多序列比對算法，表明該方法具有廣泛的適用性。

參考文獻：

[1]LASSMANN T, SONNHAMMER E L L. Quality assessment of multiple alignment programs[J]. FEBS, 2002,529:126-130.

[2]OSAMU G. Multiple sequence alignment: algorithm and applications[J]. Adv Biophys, 1999,36: 159-206.

[3]THOMPSON J D, PLEWNIAK F, POCH O. A comprehensive comparison of multiple sequence alignment programs[J]. Nucleic Acids Res, 1999,27: 2682-2690.

[4]KARPLUS K, HU B. Evaluation of protein multiple alignments by SAM-T99 using the BALIBASE multiple alignment test set[J]. Bioinformatics, 2001,17:713-720.

[5]MORGENSTERN B. DIALIGN: multiple DNA and protein sequence alignment at BiBiServ[J]. Nucl Acids Res, 2004,32: 33-36.

[6]CEDRIC N, DESMOND G H, JAAP H. T-Coffee: a novel method for fast and accurate multiple sequence alignment[J]. J Mol Biol, 2000, 302:205-217.

[7]YANG Jing, LI Chengyun, WANG Yunyue, et al. Computational analysis of signal peptide-dependent secreted proteins in saccharomyces cerevisiae[J]. Agricultural Sciences in China, 2006, 5 :221-227.

[8]谷俊峰,王希誠,趙金城.多序列漸進式比對算法研究與比較[J]. 生物信息學, 2005, 2:73-76.

[9]尚驊.代數運算與自然數[EB/OL].(2002-05-23)[2010-04-12].http://media.open.edu.cn/media_file/rm/ip2/2002_5_23/gdds/gdds2/htm/gdds05.htm.

[10]MARC S, KENNETH S. Permutation distance measures for memetic algorithms with population management[C]. The Sixth Metaheuristics International Conference.Vienna, Austria, 2005.

[11]HIRSCHBERG D S. A linear space algorithm for computing maximal common subsequences [J]. Communications of the ACM, 1975,18: 341-343.

[12]張陽, 李建良,胡正國. News Grouper: 一個自動抽取重要新聞的軟件工具[J]. 計算機工程, 2002,28 :83-84.

[13]Phylip軟件包[EB/OL].[2010-04-12].http://evolution.gs.washington.edu/phylip.html.

[14]JONES D T, TAYLOR W R, THORNTON J M. The rapid generation of mutation data matrices from protein sequences[J]. Comput Appl Biosci, 1992,8 :275-282.