串聯重復序列識別方法研究

摘要：非編碼區重復序列分析在基因組研究中起著重要作用，其基礎就是在非編碼DNA序列中識別和定位所有的重復結構。重復序列識別問題在計算機科學中主要體現為字符串匹配問題。在分析了后綴樹和后綴數組字符串匹配算法的基礎上，詳細闡述了基于后綴數組的精確串聯重復序列識別方法。實驗表明，該方法適合用于非編碼DNA序列分析。

關鍵詞：串聯重復序列；后綴樹；后綴數組；最大串聯重復序列

中圖分類號：TP274文獻標識碼：A 文章編號：1009-3044(2008)31-0930-02

The Research of Methods for Identifying Tandem Repeat

CHEN Chang-ping， LIU Zi-wei， ZHOU Wen-juan， PENG Chun-yan

(Southwest University of Science and Technology， College of Computer Science and Technology， Mianyang 621000， China)

Abstract: The repeat sequences analysis of non-coding area plays an important role in the research of genomes， its foundation is to identify and locate the periodic patterns. It addresses the method of identifying the accurate tandem repeat in detail after analyzing suffix tree and suffix array algorithms of string matching. The experiment indicates that the method adapts to non-coding DNA sequence analysis.

Key words: tandem repeat; suffix tree; suffix array; longest tandem repeat

1 引言

研究表明，在高等生物和人的基因組中非編碼區都占到基因組序列的絕大部分。這說明，非編碼序列的存在與真核生物基因表達調控密切相關。絕大部分非編碼序列以高度重復序列的形式存在，如衛星、小衛星、微衛星、長散布元件、短散布元件等。從物種進化的觀點來看，發生染色體數目改變的最有可能的機制就是在特定的點發生斷裂。產生這種特定的點的方式極有可能是非編碼區DNA重復序列的加倍，而這種加倍的產生必定是由具有一定催化功能的非編碼DNA功能元件調控的。根據這個機制，可以先識別出非編碼DNA序列中的重復序列，然后再搜索其前的功能元件。因此，識別非編碼區的重復序列是分析非編碼區功能結構特征的基礎。

串聯重復序列(tandem repeat)是指以一定的堿基數作為重復單元，首尾相連排列在一起形成聚集區的重復序列。串聯重復序列識別問題在計算機科學中主要體現為字符串匹配問題，在分析了后綴樹和后綴數組字符串匹配算法的基礎上，詳細闡述了基于后綴數組的精確串聯重復序列識別方法。

2 字符串匹配算法

在計算機科學領域，字符串匹配算法一直都是研究焦點之一。比較經典的字符串匹配算法有BF[1](Brute-Force)算法，KMP[1]算法，后綴樹和后綴數組算法等。前兩種算法廣泛用于單模式匹配中，后兩種算法用于多模式匹配、最長回文串查找等方面。

2.1 后綴樹算法

后綴樹將目標字符串的所有后綴按字符的字典順序構造成“字典樹”。此后綴樹的每一個節點(除葉子結點外)，都有可能有|Σ|個分支。因此，從后綴樹的根結點出發，按照對應的分支向下查找，最多需要O(m)的比較，就可以確定模式串是否在目標字符串中出現。但是構造此后綴樹的時空花費為O(n2)。為了減少構造后綴樹的時空花費，E.M.Mclreight[2]提出一種采用路徑壓縮方法構造后綴樹的方法，即每條邊不只是代表一個字符，還有可能是一個字符串，它的時間復雜度為O(n)，但它是非在線的。于是E.Ukkonen[3]又提出了在輸入的同時構造出后綴樹的方法，通過分析可以知道基于后綴樹的字符串匹配算法的時間復雜度為O(m+n)。但是當它應用到實際的基因組序列中時，內存大小的限制成為主要瓶頸。因此，研究者開始尋找一個比后綴樹更節省空間的結構－后綴數組。

2.2 后綴數組算法

后綴數組SA(Suffix array)是一個一維數組，它保存的是字符串S的n個后綴從小到大進行排序之后的后綴的開頭位置。后綴數組與后綴樹的思想基本相同，但它比較容易理解，也易于編程實現，更重要的是它占用的空間僅是后綴樹的1/3~1/5。因此，它適用于基因組非編碼區重復序列的識別。目前構造后綴數組普遍使用的方法就是倍增算法，它充分利用了各個后綴之間的聯系，將構造后綴數組的最壞時間復雜度成功降至O(nlogn)。

3 倍增算法思想

3.1 基本概念

定義3.1 一個字符表Σ是一個建立了全序關系的集合，字符串S是將n個字符順次排列形成的數組，n稱為S的長度。S的第i個字符表示為S[i-1]，且S[n-1]='$'。'$'小于Σ中的任何一個字符。除了S[n-1]之外，S中的其他字符都屬于Σ。字符串S從i開始的后綴表示為Suffix(i)，即Suffix(i)=S[i…|S|-1]。

定義3.2 名次數組：Rank=SA-1，也就是說若SA[i]=j，則Rank[j]=i。

定義3.3 對于一個字符串u，定義u的k-前綴，k-前綴的比較關系

性質3.1對k≥n-1，Suffix(i)

性質3.2 Suffix(i)=2k Suffix(j)等價于Suffix(i)=k Suffix(j)且Suffix(i+k)=k Suffix(j+k)

性質3.3 Suffix(i)≤2k Suffix(j)等價于Suffix(i)

定義3.4 k-后綴數組：SAk保存0...n-1的某個排列SAk[0]，SAk[1]，…，SAk[n-1]使得Suffix(SAk[i])≤k Suffix(SAk[i+1])，0≤i

定義3.5 k-名次數組Rankk：Rankk[i]代表Suffix(i)在k-前綴關系下從小到大的“名次”，

3.2 構造方法

假設已經求出了SAk和Rankk，根據性質3.2和3.3，2k-前綴比較關系可以由常數個k-前綴比較關系組合起來等價地表達，因此可以很方便地求出SA2k和Rank2k，而Rankk數組實際上給出了在常數時間內進行

Suffix(i)

Suffix(i)=k Suffix(j)當且僅當Rankk[i]=k Rankk[j]

因此，對所有的后綴在≤k關系下進行排序就相當于每個Suffix(i)有一個主關鍵字Rankk[i]和一個次關鍵字Rankk[i+k]。采用快速排序，則從SAk和Rankk構造出SA2k的時間復雜度為O(nlogn)，然后可以在O(n)時間內根據SA2k構造出Rank2k。

由定義3.3可知，當k＝1時，實際上就是把每個后綴按照它的第一個字符進行排序就可以求出SA1。這樣，依次求出SA2和Rank2，SA4和Rank4，SA8和Rank8，…直到SAm和Rankm，其中m=2k且m≥n-1。而根據性質3.1，SAm和SA是等價的。這樣一共需要進行logn次O(n)的過程，因此，可以在O(nlogn)的時間內計算出后綴數組SA和名次數組Rank。

3.3 最長公共前綴

最長公共前綴是后綴數組的一個輔助工具。

定義3.6 lcp(i，j)=(Suffix(i)，Suffix(j))，表示后綴Suffix(i)和Suffix(j)的最長公共前綴的長度。數組LCP表示后綴數組中每個后綴和它的祖先間的最長公共前綴的長度，即LCP[i]=lcp(SA[i-1]，SA[i])，1≤i≤n-1，約定LCP[0]=0。

性質3.4對于任意的0≤i

設h[i]=LCP[Rank[i]]，即LCP[i]=h[SA[i]]。h數組滿足下面的性質：

性質3.5 對于i>0且Rank[i]>1，一定有h[i]≥h[i-1]-1。

根據性質3.5，可以令i從0循環到n-1按照如下方法依次算出h[i]：

若Rank[i]=0，則h[i]=0。字符比較次數為0。

若i=1或者h[i-1]≤1，則直接將Suffix(i)和Suffix(Rank[i]-1)從第一個字符開始依次比較直到有字符不相同，由此計算出h[i]。

否則，說明i>1，Rank[i]>1，h[i-1]>1，根據性質3.5，Suffix(i)和Suffix(Rank[i]-1)至少有前h[i-1]-1個字符是相同的，于是字符比較可以從h[i-1]開始，直到某個字符不相同，由此計算出h[i]。

整個算法的復雜度為O(n)，求出了h數組，根據關系式LCP[i]=h[SA[i]]可以在O(n)時間內求出LCP數組。

從空間效率來看，計算后綴數組SA和名詞數組Rank總的空間占用為4n個整數，加上LCP數組共需要5n個整數空間。因此，從空間復雜度考慮，后綴數組比后綴樹更適合用于非編碼DNA序列分析中串聯重復序列的識別。

4 精確串聯重復序列識別

精確串聯重復序列是指核心序列以不間斷的重復方式首尾相連而成，字符串之間完全匹配。

4.1 問題描述

設字符串S=s0s1...sn-1，字母表Σ={A，C，G，T}或Σ={a，c，g，t}，S的長度為|S|=n。對于i≤j，子串S[i…j]表示從位置i開始到位置j結束的字符串，若S[i…j]滿足前三個定義，則稱它為串聯重復序列。

定義4.1 重復周期：如果?坌i，0≤i

定義4.2 重復單元：若r是S的最小周期，則子串sisi+1...si+r-1稱為重復序列的重復單元。例如S=ATATATAT，此重復序列的最小周期是2，重復單元為AT。

定義4.3 重復次數：e=|S|/r，當e≥2時，稱e為S的重復次數。

定義4.4 最大串聯重復序列：令S[i…j]是最小周期為r的串聯重復序列，若S[i-1]≠S[i-1+r](i≥1)且S[j+1]≠S[j+1-r](j≤n-2)，則稱S[i…j]是最大串聯重復序列。把滿足第一個條件的S[i…j]稱為左最大重復序列，滿足第二個條件的S[i…j]稱為右最大重復序列。

下面介紹利用后綴數組進行精確串聯重復序列識別，查找出滿足用戶指定條件的所有最大串聯重復序列。

4.2 理論基礎

根據后綴數組和LCP的定義，發現一個串聯重復序列中的所有拷貝是排列在一起的。因此，通過分析后綴數組和LCP數組元素之間的關系可定位所有的串聯重復序列。先介紹下面幾個引理[3]。

引理4.1 設d=|SA[i-1]-SA[i]|，如果d≤LCP[i]，則必定存在一個從位置i=min{SA[i-1]，SA[i]}開始的右最大重復序列。這個重復序列的長度為length=d+LCP[i]。其中d是這個重復序列的周期，但不一定是最小周期。

根據性質3.4，將引理4.1進行擴展，得到引理4.2。

引理4.2 設d=|SA[i]-SA[j]|(i

利用引理4.1和引理4.2，可以找到所有從任意位置開始的右最大重復序列。然而，如果試圖找到開始于每個位置的最大重復序列，則空間消耗非常大。

引理4.3 設d=|SA[i]-SA[j]|(i

由引理4.3可知找出所有長度大于ml的最大串聯重復序列，不必從每個位置開始搜索，而是先定位LCP數組中LCP[i]<(ml+1)/2的元素，則剩下的就是LCP[i]≥ml/2的元素。

4.3 具體算法

步驟1讀取參數minperiod、mintime、minlength和DNA序列S

步驟2 計算S的后綴數組SA[i]和LCP數組LCP[i](0≤i≤n-1)

步驟3 當LCP[i]<(ml+1)/2時，標記i，則未標記的LCP[j]滿足LCP[j]≥ml/2。

步驟4 若兩個標記過的LCP數組元素之間有k個未標記的元素LCP[m]，LCP[m+1]，…LCP[m+k-1]，則對任意的j(m-1

步驟5 判斷MR是否滿足最小長度和最小周期的設置，若滿足，則根據MR輸出相應的周期、重復次數和重復單元。

步驟6 重復步驟4和步驟5，直到處理完所有未標記的元素。

算法結束

4.4 實驗及其結果分析

參數minperiod、mintime、minlength分別表示串聯重復序列的最小周期、最小循環次數、最小長度。上述算法用VC++6.0編程實現，并應用于小鼠X染色體1720條非編碼DNA序列的串聯重復序列查找。當參數設置為minperiod=5，mintime=4，minlength=20時，共找到4387條重復序列，其分布如表1所示：

表1 小鼠X染色體非編碼區重復序列分布

■

實驗結果表明，在非編碼區大量存在五、六堿基重復序列類型，這與Tóth G等人[5]的研究結果一致。同時，重復序列的數量隨著重復單元堿基數目的增加而減少，這是由于串聯重復序列在生物機體內很不穩定，經常發生變異，即重復單位的擴展或縮小。

5 結束語

本文將基于后綴數組的字符串匹配算法用來識別非編碼DNA序列中的串聯重復序列，該方法的時空復雜度都能滿足實際應用的需要。下一步工作就是根據串聯重復序列的位置特征，尋找其前具有一定催化功能的序列，進而分析非編碼區功能結構特征。

參考文獻：

[1] 朱戰立.數據結構—使用C++語言[M].西安:西安電子科技大學出版社.2001:129-134.

[2] Mclreight E M. A Space-Economical Suffix Tree Construction Algorithm[J]. Jrul of Algorithms， 1976，23(2):262-272.

[3] Ukkonen E. On-line Construction of Suffix Trees[J]. Algorithmica， 1995，14(3):249-260.

[4] 王曉敏，王正志.基于后綴列的基因序列最大串聯重復查找技術(英文)[J].生物信息學，2006，(2):72-75.

[5] Tóth G，Gáspári Z，Jurka J.Microsatellites in different eukaryotic genomes:Survey and analysis[J].Genet Res，2000，(76):323-326.