趙 陽

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

0 引言

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

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

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

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

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

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

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

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

表1 各堿基的EIIP值

1.2 基因序列的去直流化

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

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

1.3 自回歸模型階次確定

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

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

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

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

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

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

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

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

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

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

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

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

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

2 實(shí)驗(yàn)及結(jié)果分析

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

2.1 實(shí)驗(yàn)過程

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

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

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

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

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

圖2 PSE識(shí)別法得到的

圖3 序列Y-27H39的譜估計(jì)

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

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

2.2 實(shí)驗(yàn)結(jié)果分析

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

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

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

表2 定位精度比較

2 結(jié)論

基于參數(shù)譜估計(jì)對(duì)已有的PSE識(shí)別方法進(jìn)行了改進(jìn)，采用堿基的EIIP作為序列數(shù)字映射的依據(jù)，大大減小了譜估計(jì)的計(jì)算量;并對(duì)參數(shù)估計(jì)方法進(jìn)行了改進(jìn)，采用改進(jìn)的協(xié)方差方法作為參數(shù)估計(jì)方法，有效避免了可能會(huì)出現(xiàn)的譜峰分裂現(xiàn)象。今后將對(duì)識(shí)別的精度作進(jì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在時(shí)間序列分析中的應(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)代信號(hào)處理教程［M］.北京:清華大學(xué)出版社，2006:52-61.

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