基于數字信號處理的串聯重復序列識別方法

趙 陽

（太原理工大學，山西 太原 030024）

0 引言

隨著許多生物基因組測序工作的完成，大量的重復序列被發現。重復基因序列在生物進化過程中起著非常重要的作用，這些重復序列在病毒和原核生物中很少出現，在真核生物中則大量存在。目前科學證實，人類基因組中大約含有50%以上的重復基因序列［1］。這些串聯重復序列可能會與一些轉錄因子結合位點相互作用改變染色體的結構或者作為蛋白質結合位點與其他蛋白質相結合。由于串聯重復序列的多態性，個體間串聯重復序列中重復拷貝的個數可能不同，因此一些串聯重復序列可以用來研究基因標識、基因圖譜、個體識別等［2］。生物學家將重復序列作為研究非編碼區的突破口，掌握基因組非編碼區規律的重要途徑。

2009年HongXia Zhou等人［3］提出了一種基于參數譜估計的串聯重復序列識別方法，PSE（Parametric Spectral Estimation）識別法，作者采用了現代頻譜估計中的自回歸模型（AR Auto-Regressive）作為功率譜估計的模型。通過對該方法的進一步分析，發現該方法在求解基因串聯重復序列時還存在不足之處。首先，基因序列的頻譜圖會出現譜峰分裂現象，不利于觀察串聯重復序列的重復周期。其次，識別速度有待提高，PSE識別法中采用二進制表示法將基因序列映射成數字序列，導致計算量增大，識別速度不高。最后，模型階次的確定準則有待改進，階次的準確與否直接關系到譜估計的分辨率。PSE識別法中根據每一種定階準則分別求出基因序列的階次，并根據實驗估算出適合該模型的階次，沒有明確指出具體應該如何為模型確定階次，導致容易出現估計誤差。

通過對參數譜估計識別法進行深入研究并結合串聯重復序列識別的理論知識，針對該方法存在的不足，本文提出了基于自回歸模型的串聯重復序列識別方法，對以上不足進行了改進。

1 基于自回歸模型的串聯重復序列識別

基于參數估計的功率譜估計是現代功率譜估計的重要內容，其目的就是為了改善功率譜估計的頻率分辨率。它主要包括自回歸模型（AR Auto-Regressive）、滑動平均模型（MA Moving Average）以及自回歸滑動平均模型（ARMA Auto-Regressive Moving Average）［4］。

由于AR模型的參數估計只需要解一組線性方程，而ARMA以及MA模型的參數估計通常需要解一組非線性方程，因此求解AR模型參數估計的過程會相對容易一些?；贏R模型的功率譜估計是最常用的一種參數估計頻譜分析方法，本文將采用AR模型識別基因序列中的串聯重復序列。其功率譜可表示為:

1.1 基因序列的數字映射

基因序列是由四種堿基（A、T、G、C）組成的字符串序列。為了在DNA序列分析中使用數字信號處理的方法，首先需要將DNA序列中的四個堿基分別映射成數字。本文根據各個堿基的EIIP［5］值將DNA序列映射成為一條數字序列。堿基的EIIP是堿基的物理屬性，表示其價電子的平均能量，可唯一表示一種堿基。該映射法可將一條基因序列惟一的映射為一條數字序列，相比二進制映射法，其計算量減少了3/4。各堿基的EIIP值如表1所示。

表1 各堿基的EIIP值

1.2 基因序列的去直流化

在求數字基因序列的譜估計時，為了避免基因信號中直

其中Ne［n］為原數字序列，N［n］為去直流化的基因數字序列，并最終采用N［n］進行階次估計及參數估計。

1.3 自回歸模型階次確定

模型階數p是一個非常重要的參數，當階數選擇過小時，會導致相應的AR模型的極值點不夠精確，表現為譜峰較少，頻譜較為平滑，頻譜分辨率下降;而當階數選擇過大時，譜估計值的分辨率雖然會有所提高，但是會產生虛假的譜峰，導致其統計特性的不穩定。因此我們需要為模型確定一個合適的階次，既要保證譜估計的分辨率較高，又不至于出現虛假譜峰。下面是幾種階數估計準則。

（1）最終預測誤差準則（FPE）流信號的干擾，這里需要對基因數字信號進行處理，將每一個數字信號分別減去整個數字信號序列的平均值得到一個新的數字序列:

（2）Akaike信息準則（AIC）

上述三式中其中L為數據長度，p為模型階數，ρp表示p階AR模型的預測誤差功率估計值。

根據已有理論可知，各準則函數取得最小值時的階次p即為AR模型階次。每一種定階準則都有其優缺點，并且對于同一輸入采用不同的準則得到的階數可能不同。為了減小因為階次選擇不準確導致的譜估計誤差，本文將三種準則所確定階次的平均值作為模型的階次，經過反復試驗發現，通過該方法確定的階次較為準確，得到的譜估計結果能夠滿足串聯重復序列識別所需要的分辨率。

1.4 AR模型的參數估計

由（1）式可知，只要得到該模型的參數估計{a（1），a（2），LL，a（p）}及噪聲方差σ2就可以得到該模型的功率譜密度。常用的模型參數估計方法主要有尤克-沃勒（Yule-Walker）算法，Burg算法，協方差算法以及改進的協方差算法。本文將采用改進的協方差算法作為AR模型參數估計的方法。改進的協方差方法克服了Burg算法的缺點，采用該方法能夠有效地避免譜線分裂現象的出現［6］。

在得到的功率譜密度圖中橫坐標表示頻率，縱坐標表示功率譜密度值。當一個基因序列中包含周期為n的串聯重復序列拷貝時，在功率譜密度圖中對應頻率為ω=2π/n處將會出現一個顯著的波峰。在得到基因數字序列的譜估計之后根據信噪比的設置從功率譜密度圖中選擇出對研究有意義的波峰，此時便可得到串聯重復序列拷貝在基因序列中出現的頻率，進而得到串聯重復序列拷貝的周期即拷貝

（3）貝葉斯信息準則（BIC）長度。

在得到的功率譜密度圖中能夠很好地顯示基因序列中的頻率信息，即圖中波峰出現的地方對應的頻率就是串聯重復序列拷貝在基因序列中出現的頻率。但是由于在功率譜密度圖中可能會出現若干個波峰，并且有些波峰并不能被認為是對識別串聯重復序列有意義的信息，因此需要有針對性地選擇一些對識別串聯重復序列有用的波峰。本文采用信噪比來確定每一個波峰的重要性，即當S（k）/Sm＞S/N成立時便認為功率譜密度圖中的波峰是有意義的。研究表明，當信噪比設置為4時，能更好地識別串聯重復序列的位置信息，因此這里將信噪比設置為4。根據信噪比確定出在功率譜密度圖中有意義的波峰，則波峰對應的頻率即為串聯重復序列拷貝在基因序列中出現的頻率，求其倒數便可得到串聯重復序列拷貝的周期。

1.5 串聯重復序列定位

在分析過基因序列的功率譜密度圖之后，得到了序列中存在的串聯重復序列拷貝的長度，這里將采用短時傅里葉變換對序列進行分析，便可得到串聯重復序列在基因序列中出現的位置［7］。序列Ne［n］的短時傅里葉變換:

其中k=0，1，L，M-1。在求基因序列的短時傅里葉變換時需要根據信號的特點選擇適合的窗函數以及窗口大小。適合的窗函數以及窗口大小能夠使得頻譜更加精確，分辨率更高。

2 實驗及結果分析

為了驗證本文方法的正確性及有效性，本文從美國國家生物技術中心（NCBI）維護的基因數據庫GenBank［8］中提取了若干條基因序列進行分析。本文將以序列Y-27H39為例進行詳細實驗過程分析。

2.1 實驗過程

首先將序列Y-27H39根據各堿基的EIIP值映射為數字序列。然后采用基于AR模型的譜估計法對數字序列進行譜估計。觀察序列中串聯重復序列拷貝出現的頻率。

在求數字序列的譜估計時，需要先確定該數字序列的階次。如圖1所示分別為采用FPE、AIC、BIC三種定階準則估計的階次結果。從圖中可以觀察到由三種定階準則確定的階次分別為:22，22，9。根據上文中關于模型階次的分析，這里將采用的階次為18。

其次，結合已確定的階次利用參數譜估計方法進行基因序列的譜估計。為了對比方便，我們分別采用已有的PSE識別法以及新提出的方法分別求出了序列Y-27H39的譜估計，如圖2及圖3所示。從二者的對比可以看出，本文提出的方法避免了PSE識別法中存在的譜峰分裂現象。

最后通過對該序列求短時傅里葉變換便可得到該序列中串聯重復序列出現的位置，如圖4所示，矩形中標注的部分便是序列Y-27H39中串聯重復序列出現的位置。

圖1 采用FPE、AIC以及BIC準則分別對序列Y-27H39定階

圖2 PSE識別法得到的

圖3 序列Y-27H39的譜估計

圖4 序列Y-27H39的短時傅里葉變換

但是PSE識別方法采用了二進制表示法進行基因序列的數字映射，在進行譜估計時將一條基因序列映射成為四條數字序列，即針對每類堿基各得到一條數字序列，需要分別對得到的四條序列求四次譜估計才可以得到最終的譜估計，計算量較大。本文提出的方法可唯一地將一條基因序列映射成一條數字序列，并針對該序列進行譜估計，計算量相對前者減少了75%，由于基因序列數據量較大，因此計算量是進行基因序列分析必須考慮的問題之一。另外，從圖2及圖3的對比可以看出，本文提出的方法解決了PSE識別法中存在的譜峰分裂現象，使得得出的串聯重復序列出現頻率更加精確，便于進一步進行串聯重復序列位置分析。

2.2 實驗結果分析

從圖2中的波峰處可以看到譜峰出現了分裂現象，但是仍能模糊地判斷出此處波峰對應的頻率大約為F=0.25 Hz，這是因為該序列是一個短基因序列，只有194個堿基。當序列長度較長或者序列中包含的串聯重復拷貝較多時就很難準確地斷定波峰處對應的頻率究竟應該是多少，這樣就對進一步識別串聯重復序列出現的位置帶來了困難。出現譜峰分裂現象的主要原因是PSE識別法中采用的模型參數估計方法不當。PSE識別法中采用了Burg算法作為參數估計的方法，不能保證對所有的基因序列都能得到其準確的功率譜估計。

觀察圖3可以發現并沒有出現譜峰分裂現象，從圖中波峰位置可以清楚地判斷該序列中串聯重復序列出現的頻率為0.25 Hz，進而可以判斷該序列中串聯重復序列出現的周期為4，即串聯重復拷貝長度為4 bp。

表2中列出了采用本文提出的識別方法對序列Y-27H39進行分析，得到的串聯重復序列位置，并與PSE識別法以及GenBank中標注的位置進行了對比。從表2中可以看出，兩種方法均能較準確地定位出串聯重復序列的位置。

表2 定位精度比較

2 結論

基于參數譜估計對已有的PSE識別方法進行了改進，采用堿基的EIIP作為序列數字映射的依據，大大減小了譜估計的計算量;并對參數估計方法進行了改進，采用改進的協方差方法作為參數估計方法，有效避免了可能會出現的譜峰分裂現象。今后將對識別的精度作進一步的研究。

［1］Lander E S，Linton L M，Birren B，et al.Initial Sequencing and Analysis of the Human Genome［J］.Nature，2001，409:860-921.

［2］Naruse K，Tanaka M，Mita K，et al.A Medaka Gene Map:the Trace of Ancestral Vertebrate Proto-chromosomes Revealed by Comparative Gene Mapping［J］.Genome Research，2004，14（5）:820-828.

［3］Zhou H X，Du L P，Yan H.Detection of Tandem Repeats in DNA Sequences Based on Parametric Spectral Estimation［J］.IEEE Transactions on Information Technology in Biomedicine，2009，13（5）:747-755.

［4］張善文，雷英杰，馮有前.Matlab在時間序列分析中的應用［M］.西安:西安電子科技大學出版社，2007:130-139.

［5］Irena Cosic.Macromolecular Bioactivity:Is It Resonant Interaction Between Macromolecules–Theory and Applications［J］.IEEE Transactions on Biomedical Engineering，1994，41（12）:1101-1114.

［6］Akhtar M，Ambikairajah E，Epps J.Comprehensive Autoregressive Modeling for Classification of Genomic Sequence［C］.Proceedings of the IEEE 6th International Conference on Information，Communications＆ Signal Processing，Singapore，2007:1-5.

［7］胡廣書.現代信號處理教程［M］.北京:清華大學出版社，2006:52-61.

［8］GenBank［OL］http://www.ncbi.nlm.nih.gov/genbank/2011.