三周期性優(yōu)選下的蛋白質(zhì)表達(dá)區(qū)域預(yù)測算法研究?

夏盛瑜 辛月振 陸志卿

（中國石油大學(xué)（華東）計(jì)算機(jī)與通信工程學(xué)院 青島 266000）

1 引言

DNA序列中帶有遺傳效應(yīng)的片段稱為基因，是遺傳的結(jié)構(gòu)和功能單位。真核生物的DNA序列中，如果把基因劃分成許多分開的片段，那么把能夠編碼蛋白質(zhì)的序列即基因的編碼區(qū)，稱作外顯子；把不能編碼蛋白質(zhì)的基因序列稱為內(nèi)含子。序列要經(jīng)過剪接、轉(zhuǎn)錄、復(fù)制過程才能合成蛋白質(zhì)。而在合成蛋白質(zhì)的過程中，外顯子攜帶的遺傳信息卻得保存，準(zhǔn)確地將遺傳信息轉(zhuǎn)遞到蛋白質(zhì)中，并參與指示各項(xiàng)生命活動(dòng)。對編碼區(qū)和非編碼區(qū)的預(yù)測是基因預(yù)測的關(guān)鍵［1］。

弄清整個(gè)DNA序列中每一段基因所攜帶的遺傳訊息對后續(xù)研究工作的進(jìn)行起著關(guān)鍵作用，同時(shí)這也是人類研究的起初目的所在。由于真核生物的外顯子和內(nèi)含子之間缺乏明顯的特征，基因組的蛋白質(zhì)編碼區(qū)域通常不連續(xù)［2］。

而原核生物，其蛋白質(zhì)編碼區(qū)域相當(dāng)于真核生物的外顯子區(qū)域。原核生物的基因預(yù)測已經(jīng)有較多方法。總體可分為:1）基于同源性搜索的外部算法，稱為間接遺傳識別方法。現(xiàn)在最流行的序列對齊工具是FASTA和BLAST［3］。2）基于內(nèi)在統(tǒng)計(jì)預(yù)測的算法。這種搜尋基因的方法與模式識別非常相似，分為基于隱式馬爾科夫模型的經(jīng)典的原核識別方法和基于神經(jīng)網(wǎng)絡(luò)的算法和基于曲線辨識分析的算法。3）混合法。它是內(nèi)在算法和外部算法的組合。

根據(jù)3堿基周期理論，大多數(shù)外顯子序列具有3堿基周期性，而內(nèi)含子序列沒有這個(gè)獨(dú)特的特征。傳統(tǒng)的三期分析方法將基因序列映射成數(shù)字序列，而蛋白質(zhì)編碼區(qū)域往往顯示出較為明顯的峰值。3周期性分析結(jié)合先驗(yàn)生物學(xué)知識，對預(yù)測結(jié)果有很好的解釋。此外，功率譜分析方法具有良好的性能和良好的穩(wěn)定性，因?yàn)樾枰绕渌椒ǜ俚膮?shù)，可以滿足個(gè)人需要［5］。王其強(qiáng)利用三周期性對抑癌基因方面做了建模和算法分析［9］。最近，有學(xué)者利用韋爾奇函數(shù)法進(jìn)一步完善了現(xiàn)有的功率譜分析方法，并提出了一種更有效的三次循環(huán)判斷模型，并利用K-means聚類［17］和相關(guān)分析與3堿基周期結(jié)。因此信號處理方法也可用于基因組DNA研究，因?yàn)樗鼈兛梢宰R別常規(guī)統(tǒng)計(jì)學(xué)方法無法檢測的隱藏特征。將DNA序列轉(zhuǎn)換為數(shù)字信號后，可以應(yīng)用一些信號處理策略。作為生物蛋白序列的典型特征，V.R.Chechetkin于1995年已經(jīng)對傅里葉諧波的大小和3堿基周期性的依賴關(guān)系做了討論［4］。

在文中提出了一種可擴(kuò)展的3堿基周期性預(yù)測方法來預(yù)測基因組蛋白編碼區(qū)的優(yōu)化方法。通過傅立葉分析，DNA序列轉(zhuǎn)化為信號量。根據(jù)NCBI庫中的實(shí)驗(yàn)驗(yàn)證蛋白質(zhì)數(shù)，利用滑動(dòng)窗口預(yù)測初始滑動(dòng)框大小和閾值。對結(jié)果進(jìn)行比較和篩選以縮小預(yù)測區(qū)域，然后使用可擴(kuò)展滑動(dòng)窗口策略來提高預(yù)測精度。

2 相關(guān)知識

2.1 三周期性

基因在編碼蛋白質(zhì)時(shí)是通過64個(gè)三聯(lián)碼（其中包含3個(gè)終止子）來確定20種氨基酸，蛋白質(zhì)編碼區(qū)具有周期3行為正是基于三聯(lián)碼的緣故［15］。即使許多外顯子不具有周期三行為，但將它們連在一起編碼蛋白質(zhì)時(shí)大都表現(xiàn)出周期三行為。三周期行為可以表述為:將基因序列映射成數(shù)字序列，蛋白質(zhì)編碼區(qū)域往往顯示出較為明顯的峰值，這個(gè)峰值往往出現(xiàn)在序列的1/3和2/3處，而非編碼區(qū)和隨機(jī)序列則表現(xiàn)為無規(guī)則序列［16］。

由于直接處理DNA序列很困難，需要將DNA序列轉(zhuǎn)化為數(shù)字信號，通過數(shù)值方法更方便處理DNA序列。這里使用的是Voss方法［7］。令S為基于長度為N的符號集I的DNA序列，其中I=｛A，T，C，G｝。 DNA 序列 S可以由 S：S［0］，S［1］，…，S［N-1］表示，其中S［i］∈I。對于任何確定b∈I。

二進(jìn)制序列uA（n），uT（n），uC（n），uG（n）的A，T，C和G的不論存在或不存在都可以在序列的第n位表示出來。

離散傅立葉變換（DFT）中，可以將信號序列轉(zhuǎn)換為頻域中的一組新值［6］。對于長度為N的信號序列，其DFT定義如下。

在k處的信號的DFT功率譜被定義為

DNA序列的DFT功率譜是其四個(gè)二進(jìn)制指示符序列的功率譜之和。

其中 PA（k），PT（k），PC（k）和 PG（k）是四個(gè)指示符序列 uA（n），uT（n），uC（n），uG（n）的傅里葉功率譜［12］。

2.2 滑動(dòng)窗口

滑動(dòng)窗口是基于DFT的一種優(yōu)化，為了便于區(qū)分和識別功率譜。對于DNA序列S，其指示序列為｛ub［n］｝。長度M設(shè)置（設(shè)置為3的整數(shù)倍）為固定窗口的長度。對于以n（對于任何n（0≤n≤N-1）為中心的長度為M（n-（M-1）/2，n+（M-1）/2）的序列片段，當(dāng)n接近開始和結(jié)束的DNA序列時(shí)，窗口實(shí)際有效長度可能小于M，那么可以有四個(gè)指令序列的 DFT［10］。

那么在M/3處的頻譜值為p（n，M/3）。

最后，使用歸一化處理光譜值的結(jié)果。將結(jié)果除以頻譜的最大值［11］。

3 相基于三周期性的蛋白質(zhì)區(qū)域預(yù)測

3.1 基本外顯子區(qū)域獲取

頻譜圖的基礎(chǔ)為三周期全序列功率圖。在此基礎(chǔ)上的頻譜圖是簡化處理的一種實(shí)驗(yàn)手段。

圖1和圖2分別為實(shí)驗(yàn)中的編碼區(qū)和非編碼區(qū)的3種周期性功率譜的總和曲線。

圖1 編碼區(qū)三周期圖

圖2 非編碼區(qū)三周期圖

通過滑動(dòng)窗口的算法在頻譜中將不斷出現(xiàn)較為明顯的波動(dòng)值，對預(yù)測有很好的幫助，對獲取外顯子區(qū)域是算法的基礎(chǔ)。對于通過滑動(dòng)窗口計(jì)算的長度為N的DNA序列，使用閾值A(chǔ)（用于區(qū)分外顯子和其他編碼區(qū)域），閾值B（用于子區(qū)域的控制長度）和頻譜序列K的子區(qū)域發(fā)現(xiàn)外顯子區(qū)域。為了提高傅里葉變換的效率，本文使用中的快速傅里葉變換［8］（FFT）。

外顯子獲取的基本過程如下:

1）使用滑動(dòng)窗口M計(jì)算DNA序列，得到P（n）的光譜序列；

2）設(shè)置K=0；

3）從起始點(diǎn)開始，P（n）與閾值A(chǔ)比較:如果P（n）＞ A，光譜序列可能區(qū)域K+1；

4）否則:如果光譜序列的子區(qū)域的長度大于閾值B，輸出核苷酸序列；

5）否則:重新計(jì)數(shù)K；

6）增加n，重復(fù)步驟3）到步驟4）。

選取過程中，對各別異常點(diǎn)需要做包容處理，因?yàn)樯飳?shí)驗(yàn)中不可避免的實(shí)驗(yàn)誤差以及生物基因突變和實(shí)驗(yàn)環(huán)境的影響，外顯子區(qū)域的獲取盡量在一定容忍誤差內(nèi)。

在步驟4）中進(jìn)行優(yōu)化:

若P（n）＜A，不再立刻重置K，而是暫時(shí)保存K值，判斷之后β個(gè)位點(diǎn)值，若連續(xù)出現(xiàn)P（n）＜A情況則重置K，若在β個(gè)位點(diǎn)內(nèi)，僅出現(xiàn)δ個(gè)不滿足情況，則此β段也為滿足序列，即當(dāng)前滿足序列長度為K+β。

3.2 可擴(kuò)展的滑動(dòng)窗口

由于基因長度的不確定性，滑動(dòng)窗口的大小應(yīng)該是不相同的，以防止預(yù)測信息的遺漏。對于一個(gè)DNA序列，應(yīng)該使用不同的滑動(dòng)窗口來獲得光譜。通過計(jì)算NCBI中的已知序列獲得滑動(dòng)窗口的大小。不同滑動(dòng)窗口下的序列將使用外顯子預(yù)測來獲得大致的外顯子區(qū)域。小滑動(dòng)窗口的結(jié)果可能包括在大滑動(dòng)窗口的結(jié)果中，因此確定該區(qū)域中外顯子的可能尺寸。下面提出了一種新的優(yōu)化策略:

1）重疊優(yōu)先

重疊率越高意味著外顯子存在的可能性就越大。所以優(yōu)先選擇高度重疊的區(qū)域，然后重疊區(qū)域的順序?qū)⒈缓喜ⅰＪ褂迷撔蛄凶鳛閮蓚€(gè)末端的堿基序列和插入片段序列。使用等距滑動(dòng)窗口操作來縮小外顯子區(qū)域。重疊區(qū)域有很多可能性，但在實(shí)際情況中，重疊率并不是很高。當(dāng)有三個(gè)或更多的重疊以滿足條件時(shí)，進(jìn)行進(jìn)一步處理。

2）滑動(dòng)框優(yōu)化方法

在獲得高重疊區(qū)域后也就獲得該區(qū)域的滑動(dòng)框范圍，根據(jù)滑動(dòng)框的起始大小，選擇固定的步長對此區(qū)域進(jìn)行不同大小的滑動(dòng)框運(yùn)算［13］。在不同的滑動(dòng)框下獲得的區(qū)域會(huì)逐漸穩(wěn)定，即當(dāng)?shù)竭_(dá)頭尾端點(diǎn)值誤差在某個(gè)閾值范圍內(nèi)時(shí)。在此，若滑動(dòng)框在高重疊區(qū)域的滑動(dòng)框長度取值為［A，B］，步長選取x，則需要進(jìn)行（B-A）/x次的的滑動(dòng)框運(yùn)算，當(dāng)高重疊區(qū)域很大時(shí)，次數(shù)會(huì)大大增加，影響預(yù)測效果。對此采用二分法思想:每次步長取當(dāng)前滑動(dòng)框長度差值的1/4，1/2，3/4，根據(jù)不同步長下的頭尾端點(diǎn)閾值變化情況選取，可減少滑動(dòng)框的運(yùn)算次數(shù)。

4 實(shí)驗(yàn)結(jié)果

桑黃實(shí)驗(yàn)原始數(shù)據(jù)由中國石油大學(xué)（華東）生物中心提供，測試數(shù)據(jù)為Fasta格式純文本文件。這些文件中的序列大小為初選數(shù)據(jù)為5.4G，測序公司統(tǒng)計(jì)優(yōu)選文本45M，通過本文方法獲取的可能序列與NCBI庫進(jìn)行了比對，其中與Translation elongation factor 1 alpha和ATP synthase subunit a蛋白質(zhì)有較好的匹配度，最高分為356，E值為3e-99，Ident為85%，表明算法有較好的可靠性。整體算法首先預(yù)處理文本數(shù)據(jù)；然后選擇閾值（閾值包括滑動(dòng)窗口的大小和外顯子區(qū)分值）；通過區(qū)域映射和滑動(dòng)窗口操作DNA序列獲取原始頻譜圖；接著對外顯子區(qū)域進(jìn)行預(yù)測；最后使用可擴(kuò)展的滑動(dòng)框策略縮小預(yù)測區(qū)域直至滿足閾值條件如表1所示（單位是堿基個(gè)數(shù)）。

表1 閾值選擇

滑動(dòng)框大小的選取結(jié)合一些先驗(yàn)知識，并結(jié)合NCBI庫中鞘氨醇單胞菌的公開數(shù)據(jù)中已知基因長度進(jìn)行的數(shù)理統(tǒng)計(jì)結(jié)果。在NCBI庫中，通過固定滑動(dòng)窗口分析不同長度的參考序列，并在序列的兩端添加背景序列的固定長度。計(jì)算序列的頭部和尾部的光譜值。光譜如圖3所示。對頭尾光譜的計(jì)算統(tǒng)計(jì)獲得符合實(shí)驗(yàn)要求的閾值。

圖3 頻譜圖中的光譜值

由于NCBI和Uniport等公開數(shù)據(jù)庫中的桑黃基因較短，所以測試桑黃二級產(chǎn)物的Phosphoenolpyruvate carboxylase蛋白質(zhì)ppc基因，888個(gè)氨基酸序列，預(yù)測蛋白質(zhì)編碼區(qū)的閾值為0.2。在大多數(shù)情況下測試結(jié)果表現(xiàn)良好，特別是當(dāng)序列較長時(shí)。在本文中，不使用傳統(tǒng)的經(jīng)驗(yàn)閾值0.4，因?yàn)椴煌纳锞哂胁煌倪z傳結(jié)構(gòu)特征。不同的窗口尺寸，不同的功率譜密度計(jì)算方法也將導(dǎo)致功率譜預(yù)測曲線幅度的變化。因此，經(jīng)驗(yàn)閾值將帶來顯著的預(yù)測偏差（在下圖中背景序列下的起始位點(diǎn)為1000）。

圖4

為了測試提出的方法，本文添加了仿真實(shí)驗(yàn)。從Uniprot數(shù)據(jù)庫中隨機(jī)選擇5種桑黃及其近緣的蛋白質(zhì)，并將它們經(jīng)過指定次數(shù)的變異后插入到長的背景序列中。經(jīng)過粗略的外顯子預(yù)測后，不同滑動(dòng)窗口的可能區(qū)域的數(shù)量如下表所示。然后使用端點(diǎn)對齊和合并重疊序列，并選擇了兩個(gè)新的序列，對于新序列1，使用300到100個(gè)堿基長度的等距滑動(dòng)窗口來縮小區(qū)域并刪除了非表達(dá)區(qū)域。對于新序列2，使用等距滑動(dòng)窗從500到200個(gè)堿基長度。根據(jù)頭尾位點(diǎn)反饋信息，選取新的滑動(dòng)框。

表2 背景序列下的測試結(jié)果

最后得到兩個(gè)可能的區(qū)域，預(yù)測結(jié)果與其中一個(gè)插入序列有良好的匹配度。而不經(jīng)過優(yōu)化的算法獲得的初步區(qū)域數(shù)量較多，雖然包含了目標(biāo)序列，但是缺乏實(shí)際應(yīng)用的意義。特別是滑動(dòng)框?yàn)?00和200的序列，其數(shù)量較多，但是只是滑動(dòng)框?yàn)?00的初選結(jié)果的一部分，重疊區(qū)域過多，影響精確度。

5 結(jié)語

本文針對生物的蛋白表達(dá)區(qū)域問題，提出了一種基于3堿基周期性蛋白表達(dá)區(qū)域的優(yōu)化方法。實(shí)驗(yàn)中，功率譜分析方法的算法性能良好，方法所需的各種參數(shù)設(shè)置比其它方法少，穩(wěn)定性好［14］。應(yīng)用通過已知的DNA序列和實(shí)際實(shí)驗(yàn)數(shù)據(jù)的預(yù)測結(jié)果經(jīng)過Blast后，部分結(jié)果數(shù)據(jù)和NCBI中的已知數(shù)據(jù)相符，有些將具有進(jìn)一步驗(yàn)證的生物實(shí)驗(yàn)。通過這些結(jié)果比對，驗(yàn)證此方法能良好地預(yù)測蛋白質(zhì)編碼區(qū)。其中可擴(kuò)展的滑動(dòng)窗口是一個(gè)重要的部分，使用區(qū)域選擇優(yōu)化策略降低數(shù)據(jù)復(fù)雜度，縮小可能區(qū)域。在下一步的研究中，將通過添加自適應(yīng)滑動(dòng)窗口來減少人為的主觀干預(yù)，使該方法的預(yù)測精度進(jìn)一步地提高。時(shí)間復(fù)雜度和誤差檢驗(yàn)［18］也是未來優(yōu)化的方向。