基于最長公共子序列的DNA序列相似性分析

吳東根， 周小安

(深圳大學 信息工程學院， 廣東 深圳 518000)

引言

人類從未停止過對生物本質的好奇和探索，隨著過去幾十年來生物分子學以及基因技術的飛速發展，人類得到的生物信息呈現爆炸性的增長趨勢。生物技術尤其是分子生物技術的快速發展，已經將人們帶入一個新的時代——基因組時代。在人類基因組計劃[1](Human Genome Project，簡稱HGP)的實施過程中，科學家們獲得了海量的人體基因信息，這極大促進了生物學尤其是基因學的發展，讓人類步入了后基因時代。

目前關于分析DNA相似性的方法已有很多，早期的有：Gibbs和Mclntre在1970年提出的點陣圖分析方法[2]、傅里葉變換分析法[3]、Lipman和Pearson在1988年提出了FASTA算法[4]、Gatlin提出的條件熵方法[5]；近期的有：三角函數法、譜動態法等等[6-11]。傳統研究思路基本分為兩個步驟：序列轉換、轉換后相似分析。在圖形表示法中，第一步將DNA序列轉化為圖形模式，如轉化為二維、三維坐標[12-13]中的曲線，第二步利用圖形構造矩陣，再用矩陣的不變量進行DNA序列相似性分析；在數值表示法中，第一步將DNA序列轉化為數值模式，如轉化為二進制數或整數[14-15]，第二步通過計算序列之間向量的歐氏距離或夾角余弦值[16]來分析DNA序列的相似性。

傳統的方法存在一定的缺陷，都需要進行序列轉換，這可能會在轉換過程中丟失原始DNA序列的信息，并增加工作量。基于此，本文研究旨在尋找并研究一種不需要序列轉換并且降低工作量的新方法，同時能達到相似性分析的要求。

1 基本原理

1.1 最長公共子序列

兩個DNA序列之間的最長公共子序列越長，說明這兩個DNA序列之間具有更多的相似性[1]。最長公共子序列(Longest Common Subsequence)是一個在一個序列集合中(通常為兩個序列)用來查找所有序列中最長子序列的問題。這與查找最長公共子串的問題不同的地方是：子序列不需要在原序列中占用連續的位置 。最長公共子序列問題是一個經典的計算機科學問題，也是數據比較程序，比如Diff工具，和生物信息學應用的基礎。最長公共子序列的定義如下：

一個序列S，如果分別是兩個或多個已知序列的子序列，且是所有匹配此條件序列中最長的，則S稱為已知序列的最長公共子序列[17]。

例如，序列1A2C3D4B56和B1D23CA45B6A的最長公共子序列為123456、12C4B6、12C456和1234B6(不一定只有一條)，最長公共子序列長度為6。實驗中只選取其中一條最長公共子序列即可。

求解序列的最長公共子序列問題首先需要求解動態規劃表。設有兩個序列：A=a[0]a[1]…a[m-1],A中任意一個元素記為a[i]，i∈[0,m-1]；B=b[0]b[1]…b[n-1]，B中任意一個元素記為b[j]，j∈[0,n-1]。A的長度為m，B的長度為n，生成大小為m×n的矩陣記為tr，tr的行數為m、列數為n。記tr[i][j]的含義為序列a[0]…a[i]與序列b[0]…b[j]的最長公共子序列的長度。從左到右，再從上到下計算矩陣tr中的各個元素。

(1)tr[i][0]的含義是a[0]…a[i]與b[0]的最長公共子序列長度。b[0]只有一個字符，所以tr[i][0]最大為1。如果a[i]=b[0]，令tr[i][0]=1，一旦tr[i][0]被設置為1，之后的tr[i+1][0]、…、tr[m-1][0]也都為1。

(2)矩陣tr第一行與步驟1同理，如果a[0]=b[j],則令tr[0][j]=1，一旦tr[0][j]被設置為1，之后的tr[0][j+1]…tr[0][n-1]也都為1。

(3)對其它位置(i,j)，tr[i][j]的值可能來自以下3種情況：

①可能是tr[i-1][j]，代表a[0]…a[i-1]與b[0]…b[j]的最長公共子序列長度。

②可能是tr[i][j-1]，代表a[0]…a[i]與b[0]…b[j-1]的最長公共子序列長度。

③如果a[i]=b[j]，還可能是tr[i-1][j-1]+1。

在這3個可能值中，選取最大的作為tr[i][j]的值即可。

tr矩陣中右下角的元素值代表A和B的最長公共子序列的長度。通過整個tr矩陣的狀態，可以得到最長公共子序列。具體方法如下：

步驟1從矩陣的右下角開始，有3種移動方式：向上、向左、向左上。假設移動的過程中，用p表示此時的行數，q表示此時的列數，用一個變量C(其中C=c[0]c[1]…)來表示最長公共子序列。

步驟2如果tr[p][q]大于tr[p-1][q]和tr[p][q-1]，說明之前計算tr[p][q]的時候，一定是選擇了決策tr[p-1][q-1]+1，可以確定a[p]等于b[q]，并且這個字符一定屬于最長公共子序列，把這個字符放進C，然后向左上方移動。

步驟3如果tr[p][q]等于tr[p-1][q]，說明之前在計算tr[p][q]的時候,tr[p-1][q-1]+1這個決策不是必須選擇的決策，向上方移動即可。

步驟4如果tr[p][q]等于tr[p][q-1]，與步驟3同理，向左方移動。

步驟5如果tr[p][q]同時等于tr[p-1][q]和tr[p][q-1]，向上移動或向左移動都可以。

通過上面的方法，就可以用C++程序得到兩個序列其中一條最長公共子序列，并得到這條最長公共子序列的長度。

1.2 相似系數

本文將一個DNA序列視為以{A、G、T、C}為符號集的信號序列，記x、y為物種的名稱，如mouse、human等等，其中物種名稱x與物種名稱y可以相同。記d(x)為x物種的DNA序列的長度，可以定義兩物種DNA序列中長度更大的那條序列的長度D(x,y)：

D(x,y) = max{d(x),d(y)}

(1)

記B(x,y)為兩個物種DNA序列的最長公共子序列長度，其中B(x,y)≤D(x,y)。可以定義兩個物種的相似系數H(x,y)，用來評判兩物種的相似程度：

(2)

其中0≤H(x,y)≤1，H(x,y)越接近1，表示x、y兩物種相似度越高，當H(x,y)等于1時，表示兩個物種是同一物種；H(x,y)越接近0，表示x、y兩物種相似度越低，當H(x,y)等于0時，表示兩個物種是不同物種。

2 分析對象

本文從NCBI數據庫(http://www.nubi.nlm.nih.gov)中找到并下載了7種病毒的DNA數據，表1中列出了本文中所需要研究的7種DNA序列片段的信息。表2引用文獻[1]的10個物種的β-globin基因第一外顯子序列信息，其中序列的長度為86到105不等。

表2 10種β-globin基因第一外顯子序列

3 實驗結果

3.1 兩物種DNA序列的最長公共子序列

根據前文所描述求解最長公共子序列的方法，可以求出兩個物種的最長公共子序列，稱為物種對的最長公共子序列。表3只列出了一部分物種對的最長公共子序列,每個物種對的最長公共子序列可能不止一條，僅選取其中任意一條即可。其它物種對的最長公共子序列的求法一樣，在此不一一列出。

表3 兩物種DNA序列的最長公共子序列

3.2 7種病毒DNA序列實驗情況

利用最長公共子序列算法，求出7種病毒DNA序列兩兩間的最長公共子序列，利用定義的相似系數，求出7種病毒DNA序列的相似系數，所得相似系數見表4。

觀察表4中的數據， H5N1(1)與H5N1(2)的相似系數是0.975，非常接近1，說明H5N1(1)與H5N1(2)相似性很高。這與實際結論一致，相似系數分析結果與事實符合，說明基于最長公共子序列的方法對判斷DNA序列的相似性有效。

3.3 10種β-globin基因第一外顯子序列實驗情況

對文獻[1]中10種不同β-globin基因第一外顯子序列進行了實驗，計算出各物種對的相似系數，所得數據見表5。

觀察表5中的數據，mouse-rat物種對、human-gorilla對和goat-bovine對的相似系數比較高，都十分接近1。在實際中，mouse和rat為鼠類，其DNA序列應該具備最大的相似性，實驗結果與實際情況相符；human和gorilla的相似系數最為接近1，即相似性最大，與實際情況相符。其次發現物種gallus和物種opossum分別與其他9個物種的相似系數比較小，這也是與實際相符的。證明本方法對相似性分析是有效的。另一方面，human-mouse對、mouse-gorilla對和gorilla-bovine對的相似系數接近0.9，比起其它物種對的相似系數，這些物種對更靠近1，說明這些物種對有很大的相似性，這與實際情況不符。從生物學的角度來看gorilla-bovine、gorilla屬于靈長目，bovine屬于偶蹄目，不是同源物種，不應該有很大的相似性，因此實驗顯示的數據是有偏差的。也許是由于這些物種在某些功能或某些結構中存在相似性。在其它的文獻里也曾出現過類似的情況，見文獻[18]。

表4 7種病毒DNA序列相似系數矩陣

表5 10種β-globin基因第一外顯子序列相似系數矩陣

將表4與表5數據對比，可以發現，當處理的DNA序列長度較大時，相似系數會更加準確，不同物種的相似系數會更加趨向0，同源物種的相似系數會更趨向1，這是本方法的優勢；相比之下，處理的DNA序列較短時，優勢就沒那么明顯，會出現不同物種的相似系數也十分趨近1的情況。另外，本文的方法既可用作相同長度DNA序列的相似性分析；也可以用作長度不同的DNA序列的相似性分析，這也是該方法的優勢。

4 結束語

本文提出一種基于最長公共子序列的DNA序列相似性分析方法，并通過實驗證明了此方法的有效性，同大多數傳統的方法相比，本方法無需進行DNA序列的圖形轉換或數值轉換，直接從原始DNA序列出發，大大降低了工作量。且方法既可以對比長度相等的DNA序列相似性，也可以對比長度不相等的DNA序列相似性。處理的序列長度越大，評判相似性的程度越好。同時也為DNA序列相似性分析增添了一種新的研究方法。本文只使用了其中一條最長公共子序列，是否能通過求得所有最長公共子序列進而加深相似性值得進一步研究。