基于數(shù)字信號處理的串聯(lián)重復(fù)序列識別方法

趙 陽

（太原理工大學(xué)，山西 太原 030024）

0 引言

隨著許多生物基因組測序工作的完成，大量的重復(fù)序列被發(fā)現(xiàn)。重復(fù)基因序列在生物進化過程中起著非常重要的作用，這些重復(fù)序列在病毒和原核生物中很少出現(xiàn)，在真核生物中則大量存在。目前科學(xué)證實，人類基因組中大約含有50%以上的重復(fù)基因序列［1］。這些串聯(lián)重復(fù)序列可能會與一些轉(zhuǎn)錄因子結(jié)合位點相互作用改變?nèi)旧w的結(jié)構(gòu)或者作為蛋白質(zhì)結(jié)合位點與其他蛋白質(zhì)相結(jié)合。由于串聯(lián)重復(fù)序列的多態(tài)性，個體間串聯(lián)重復(fù)序列中重復(fù)拷貝的個數(shù)可能不同，因此一些串聯(lián)重復(fù)序列可以用來研究基因標(biāo)識、基因圖譜、個體識別等［2］。生物學(xué)家將重復(fù)序列作為研究非編碼區(qū)的突破口，掌握基因組非編碼區(qū)規(guī)律的重要途徑。

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

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

1 基于自回歸模型的串聯(lián)重復(fù)序列識別

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

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

1.1 基因序列的數(shù)字映射

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

表1 各堿基的EIIP值

1.2 基因序列的去直流化

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

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

1.3 自回歸模型階次確定

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

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

（2）Akaike信息準(zhǔn)則（AIC）

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

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

1.4 AR模型的參數(shù)估計

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

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

（3）貝葉斯信息準(zhǔn)則（BIC）長度。

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

1.5 串聯(lián)重復(fù)序列定位

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

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

2 實驗及結(jié)果分析

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

2.1 實驗過程

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

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

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

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

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

圖2 PSE識別法得到的

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

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

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

2.2 實驗結(jié)果分析

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

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

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

表2 定位精度比較

2 結(jié)論

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

［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在時間序列分析中的應(yīng)用［M］.西安:西安電子科技大學(xué)出版社，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］胡廣書.現(xiàn)代信號處理教程［M］.北京:清華大學(xué)出版社，2006:52-61.

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