一種提高外顯子預(yù)測(cè)的改進(jìn)周期3消噪策略

杜竹青

(江蘇科技大學(xué) 計(jì)算機(jī)科學(xué)與工程學(xué)院,江蘇 鎮(zhèn)江 212003)

DNA是生物遺傳信息的載體,其化學(xué)名稱(chēng)為脫氧核糖核酸.DNA序列由腺嘌呤(A),鳥(niǎo)嘌呤(G),胞嘧啶(C),胸腺嘧啶(T)這4種核苷酸符號(hào)按一定的順序連接而成.其中帶有遺傳訊息的DNA片段稱(chēng)為基因.在真核基因中,基因通常被劃分為許多間隔的片段,其中編碼蛋白質(zhì)的部分稱(chēng)為外顯子,不編碼的部分稱(chēng)為內(nèi)含子.外顯子在DNA序列剪接后仍然會(huì)被保存下來(lái),并可在蛋白質(zhì)合成過(guò)程中被轉(zhuǎn)錄、復(fù)制合成為蛋白質(zhì).DNA序列通過(guò)遺傳編碼來(lái)儲(chǔ)存信息,指導(dǎo)蛋白質(zhì)的合成,把遺傳信息準(zhǔn)確無(wú)誤地傳遞到蛋白質(zhì)上并實(shí)現(xiàn)各種生命功能[1].

隨著世界人類(lèi)基因組工程計(jì)劃的順利完成,通過(guò)物理或數(shù)學(xué)的方法從大量的DNA序列中獲取豐富的生物信息,對(duì)生物學(xué)、醫(yī)學(xué)、藥學(xué)等諸多方面都具有重要的理論意義和實(shí)際價(jià)值,也是目前生物信息學(xué)領(lǐng)域的一個(gè)研究熱點(diǎn)[2].文中使用一個(gè)算法來(lái)提高現(xiàn)有模型的檢測(cè)性能,同時(shí)不影響模型的獨(dú)立性.該算法簡(jiǎn)稱(chēng)為后處理算法(post-processing algorthms,PPA),主要用來(lái)處理變換的周期3系數(shù),并實(shí)現(xiàn)周期3去噪.通過(guò)實(shí)驗(yàn)分析,驗(yàn)證了PPA算法能夠有效地提高外顯子的預(yù)測(cè)精度,同時(shí)基于小波變換的策略確保了模型的獨(dú)立性.

1 DNA數(shù)字映射及周期3頻譜

在DNA序列研究中,首先需要把A,T,G,C這4種核苷酸的符號(hào)序列,根據(jù)一定的規(guī)則映射成相應(yīng)的數(shù)值序列,以便于對(duì)其作數(shù)字處理[3-4].值+1和-1被分配給A-T和C-G堿基對(duì),形成數(shù)字序列x(n).

對(duì)數(shù)字序列x(n)作離散Fourier變換(DFT)[5-6].

k=0,1,…,L-1

(1)

式中：w(n)是窗口函數(shù)[7].窗口函數(shù)是由一個(gè)堿基對(duì)沿長(zhǎng)度為N的基因組序列從0 到N-1滑動(dòng)而成.由堿基的周期3性[8]知，外顯子序列的功率譜曲線在頻率k=N/3處具有較大的頻譜峰值,而內(nèi)含子則沒(méi)有類(lèi)似的峰值.利用此特性獲得正向頻譜成分[2]：

FSC=|XF[k]|2

(2)

結(jié)合文獻(xiàn)[9]中的零填充方案,將滑動(dòng)窗函數(shù)的位置加入式(2)中得：

FSC[n]=|XF[n,k]|2

n=0,1,2,…,N-1

(3)

式中：N是DNA序列的長(zhǎng)度.

為了提高FSC算法的預(yù)測(cè)準(zhǔn)確性,Akhtar等人提出了成對(duì)頻譜分析算法(PSC)[9].在PSC算法中,分別從DNA序列的正反兩個(gè)方向分析滑動(dòng)窗口的DFT算法.表達(dá)式(4)給出了PSC算法:

PSC[n]=|XF[n,k]|2+|XR[n,k]|2

(4)

式(4)中，XF[n,k]和XR[n,k]是指示序列x(n) 在正反兩個(gè)方向的傅里葉變換(k=L/3,n為序列長(zhǎng)度).

基于窗口的DFT方法[10-11]，依賴(lài)于窗口形狀和大小,易導(dǎo)致預(yù)測(cè)精度降低.為了克服DFT算法的缺點(diǎn),文獻(xiàn)[12]中提出的基于改進(jìn)的Gabor變換(MGWT)方案.此方案使用不同的尺度系數(shù),算法定義為：

(5)

式(5)中：n為DNA序列的位置變量,a為尺度系數(shù),w0取序列N/3處的頻率.MGWT的信號(hào)u(t)可表示為：

(6)

式(6)中信號(hào)u(t)的平方功率譜為:

M(n,a)=|U(n,a)|2

(7)

對(duì)于長(zhǎng)度為N的DNA序列,MGWT算法中的周期3分量即為投影頻譜,可通過(guò)如下關(guān)系式獲得:

(8)

該頻譜平均值涵蓋了不同大小的窗口,將其投射到位置區(qū)域,投影的局部最大值即為外顯子區(qū)域.

2 后處理算法(PPA)及應(yīng)用

PPA算法的流程如圖1.數(shù)字x(n)作為帶阻濾波器的輸入,濾波器的輸出為xf(n).S1(n)和S2(n)分別代表濾波前后x(n)的周期3幅度值.在外顯子的區(qū)域中存在堿基的周期3性,而內(nèi)含子區(qū)域沒(méi)有.因而S1(n)和S2(n)存在明顯差異.若S1(n)的和S2(n)之間的差異小于選定的閾值電平DTL,則假定此處序列為內(nèi)含子區(qū)域.

圖1 預(yù)處理算法Fig.1 Preprocessing algorithm

內(nèi)含子區(qū)域確定后,令差分信號(hào)SD(n)=S1(n)-S2(n)為內(nèi)含子區(qū)域的信號(hào)值,外顯子區(qū)域的信號(hào)電平保持不變.SD(n)的替換相對(duì)提升了外顯子區(qū)域的信號(hào)電平.

為獲得最佳的閾值電平,從數(shù)據(jù)集中選取部分小序列作為PPA算法的輸入,求得PPA算法的輸出.利用PPA的輸出繪制不同的ROC曲線,ROC曲線下的最大值即為DTL的值.DTL值確定之后,PPA算法將被應(yīng)用到提高FSC,PSC,MGWT算法的預(yù)測(cè)精度.

2.1 FSC

S1(n)和S2(n)的計(jì)算式:

0≤n≤L-1

(9)

0≤n≤L-1

(10)

根據(jù)S1(n),S2(n)和DTL,求得信號(hào)SR(n).

2.2 PSC

窗口首先向正方向滑動(dòng).S1(n)和S2(n)分別由式(9,10)計(jì)算,PPA的輸出作為正方向的頻譜值,即|XF[N,F]|2.在序列末端,窗口反方向滑動(dòng).再次使用式(9,10),計(jì)算S1(n)和S2(n).此時(shí),求得反向頻譜值|XR[N,F]|2.PSC算法和PPA算法通過(guò)式(4)結(jié)合形成PSC-PPA方案.

2.3 MGWT

對(duì)于MGWT,|U(n,a)|是序列周期3分量在尺度系數(shù)為a時(shí)的幅值.對(duì)給定的n,S1(n)是小波系數(shù)在所有尺度下的求和,即

(11)

U(n,a)由式(6)計(jì)算而得.同理,求得S2(n).將S1(n)和S2(n)應(yīng)用到PPA算法中形成MGWT-PPA方案.

3 數(shù)據(jù)集和性能評(píng)估方案

為評(píng)定PPA算法的性能,文中用到了3個(gè)數(shù)據(jù)集,即HMR195[12-13],Burset/Guigo(BG)570,Asp67.HMR195數(shù)據(jù)集包含了人、小鼠和大鼠共195條基因.其中人、小鼠、大鼠的比例為103∶82∶10.BG570數(shù)據(jù)集包含570條脊椎動(dòng)物的多外顯子基因序列.Asp67數(shù)據(jù)集包含多煙曲霉的67條混合基因序列.實(shí)驗(yàn)中,將核苷酸水平上的預(yù)測(cè)指標(biāo)分為4類(lèi)(圖2):其中，TP代表實(shí)際編碼區(qū)序列被預(yù)測(cè)為編碼區(qū);FP代表實(shí)際非編碼區(qū)序列被預(yù)測(cè)為編碼區(qū)；FN代表實(shí)際編碼區(qū)序列被預(yù)測(cè)為非編碼區(qū)；TN代表實(shí)際非編碼區(qū)序列被預(yù)測(cè)為非編碼區(qū).基因結(jié)構(gòu)預(yù)測(cè)方案的評(píng)價(jià)使用了特異性、靈敏度曲線和ROC曲線為標(biāo)準(zhǔn).關(guān)于外顯子-內(nèi)含子的分界問(wèn)題,ROC曲線探討了分類(lèi)閾值的調(diào)整對(duì)真陽(yáng)性和假陽(yáng)性比率的影響.

圖2 核苷酸水平上的4類(lèi)預(yù)測(cè)指標(biāo)Fig.2 Four categories of predictors in nucleotide level

4 實(shí)驗(yàn)結(jié)果討論和比較

文中PPA分別被應(yīng)用到FSC,PSC和MGWT3種方案中檢測(cè)外顯子.該算法的性能已在數(shù)據(jù)集BG570,HMR195和Asp67上進(jìn)行了評(píng)估.計(jì)算式(9)和(10)中用到的矩形窗口的大小為351個(gè)堿基對(duì).對(duì)FSC和PSC方案而言,需計(jì)算滑動(dòng)窗口的DFT.方案MGWT中的R度系數(shù)[12]分別為0.2和0.7.文中使用二階全通IIR帶限濾波器,濾波器的傳遞函數(shù)在Re±jθ處的極點(diǎn)由式(12)給出.根據(jù)合成濾波器的性質(zhì),濾波器的中心頻率為2π/3 Hz,R值取為0.992.

(12)

為獲取FSC-PPA,PSC-PPA,MGWT-PPA 3種算法的DTL的值.從HMR195數(shù)據(jù)集中任選24條小序列作為PPA算法的輸入，并以PPA算法的輸出SD(n)為參考DTL值繪制ROC曲線.對(duì)于FSC-PPA,PSC-PPA算法，由于二者在求取DTL時(shí)使用的方法相同，因而可得到相同的ROC曲線，如圖3所示.當(dāng)DTL為20時(shí)可得最大的AUC，因而FSC-PPA和PSC-PPA算法中的DTL值為20.對(duì)于MGWT-PPA算法，其ROC曲線如圖4所示.由圖4知，DTL值為5時(shí)，最大的AUC值為0.753.因此，MGWT-PPA算法中的DTL值為5.

圖3 FSC-PPA和PSC-PPA算法的ROC曲線Fig.3 ROC curves of FSC-PPA and PSC-PPA algorithm

圖4 MGWT算法的ROC曲線Fig.4 ROC curve of MGWT algorithm

為了進(jìn)一步分析改進(jìn)算法的性能,文中給出了功率譜對(duì)比圖.首先計(jì)算FSC,PSC,MGWT算法中序列的功率譜，然后計(jì)算加入后處理算法后序列的功率譜,即FSC-PPA,PSC-PPA,MGWT-PPA算法的功率譜.兩類(lèi)功率譜的比較分別如圖5～7.從圖中可以看出，PPA算法在內(nèi)含子區(qū)域的噪聲抑制作用.

圖5 FSC和FSC-PPA功率譜圖Fig.5 Power spectrum of FSC and FSC-PPA

圖6 PSC和PSC-PPA功率譜圖Fig.6 Power spectrum of PSC and PSC-PPA

圖7 MGWT和MGWT-PPA功率譜圖Fig.7 Power spectrum of MGWT and MGWT-PPA

方案的ROC曲線以及PPA算法處理后的功率譜對(duì)比圖都已給出.根據(jù)圖中曲線可求得不同RFP下外顯子檢測(cè)百分比,如表1.

最后,文中就方案的靈敏度和特異性?xún)煞矫孢M(jìn)行比較.靈敏度SN,表示預(yù)測(cè)成功的核酸序列中編碼區(qū)的比例;特異性Sp表示預(yù)測(cè)為編碼區(qū)的序列中，實(shí)際為編碼區(qū)的比例.3個(gè)數(shù)據(jù)集的FSC,PSC和MGWT方案以及加入PPA后的FSC-PPA,PSC-PPA和MGWT-PPA方案,他們的特異性SP與靈敏性SN關(guān)系曲線分別如圖8～10.由圖知,引入后處理后方案的檢測(cè)性能得到了提高.

表1 3個(gè)數(shù)據(jù)集上FSC,PSC,MGWT與其處理后的結(jié)果比較Table 1 Comparison of FSC,PSC and MGWT results with their post-processed values for the three datasets

圖8 數(shù)據(jù)集HMR195上SP與SN關(guān)系Fig.8 SP of SN on HMR195 Dataset

圖9 數(shù)據(jù)集BG570上SP與SN關(guān)系Fig.9 SP of SN on BG570 Dataset

圖10 數(shù)據(jù)集Asp67上SP與SN關(guān)系Fig.10 SP of SN on Asp67 Dataset

5 結(jié)論

文中使用了一種改進(jìn)的PPA算法,主要處理周期3系數(shù).該算法抑制了由偽周期得到的無(wú)關(guān)峰,從而提高了外顯子的預(yù)測(cè)精度.在噪聲抑制過(guò)程中,該算法并不需要大量的訓(xùn)練集獲得經(jīng)驗(yàn)參數(shù).該算法中唯一需憑經(jīng)驗(yàn)確定的參數(shù)為DTL值.DTL值可以從HMR195序列中任意選取一段小序列而獲得.將此DTL值應(yīng)用到數(shù)據(jù)集BG570和Asp67上,同樣得到了令人滿意的結(jié)果.

[1] 邵建峰,嚴(yán)曉華,邵偉,等.DNA序列信號(hào)3-周期特性[J].南京工業(yè)大學(xué)學(xué)報(bào):自然科學(xué)版,2012,34(4):133-137.

Shao Jianfeng,Yan Xiaohua,Shao Wei,et al.Signal 3-cyclical nature of the DNA sequence[J].JournalofNanjingUniversityofTechnology:NaturalScienceEdition,2012,34(4):133-137.(in Chinese)

[2] 周立前.基因識(shí)別算法研究與基因組進(jìn)化分析[D].湖南湘潭:湘潭大學(xué),2008:27-30.

[3] 饒妮妮,邱麗君.DNA序列數(shù)值映射方法的研究[J].生物醫(yī)學(xué)工程學(xué)雜志,2005,22(4):681-685.

Rao Nini,Qiu Lijun.The study of numerical mapping method of DNA sequences[J].JournalofBiomedicalEngineering,2005,22(4):681-685.(in Chinese)

[4] 王宏漫,歐宗瑛.一種新的DNA序列映射規(guī)則及其分析應(yīng)用[J].信號(hào)處理,2002,18(2):133-136.

Wang Hongman,Ou Zongying.A new mapping rules of DNA sequence and their analysis application[J].SignalProcessing,2002,18(2):133-136.(in Chinese)

[5] 楊晉生,熊積慧.引入閾值和判決指導(dǎo)的DFT信道估計(jì)算法[J].系統(tǒng)工程與電子技術(shù),2012,34(2):375-378.

Yang Jinsheng,Xiong Jihui.Introduction of thresholds and decision directed channel estimation algorithm of DFT [J].SystemsEngineeringandElectronics,2012,34(2):375-378.(in Chinese)

[6] 巴特爾,江彬,范曉俊.基于對(duì)稱(chēng)擴(kuò)展DFT變換的OFDM系統(tǒng)信道估計(jì)方法[J].東南大學(xué)學(xué)報(bào):自然科學(xué)版,2012,42(1):1-6.

Ba Te′er,Jiang Bin,Fan Xiaojun.Symmetric extension DFT-based channel estimation for OFDM systems [J].JournalofSoutheastUniversity:NaturalScienceEdition,2012,42(1):1-6.(in Chinese)

[7] Rapuano s,Harris F J.An introduction to FFT and time domain windows[J].IEEEInstrumentation&MeasurementMagazine,2007,10 (6):32-44.

[8] 田元新,陳超,鄒小勇.外顯子周期三行為特征的研究[J].化學(xué)學(xué)報(bào),2005,63(13):1215-1219.

Tian Yuanxin,Chen Chao,Zhou Xiaoyong.The study of Period-3 behavior of exon [J].ActaChimicaSinica,2005,63(13): 1215-1219.(in Chinese)

[9] Akhtar M， Epps J,Ambikairajah E.Signal processing in sequence analysis:advances in eukaryotic gene prediction[J].JSelTopSignalProcess,2008,2 (3):310-321.

[10] Akhtar M,Ambikairajah E,Epps J.Optimizing period-3 methods for eukaryotic gene prediction[C]∥IEEEInternationalConferenceonAconstics,SpeechandSignalProcessing,Phoenix,Arizona.USA:[s.n.],2008: 621-624.

[11] Gunawan T S.On the optimal window shape for genomic signal processing[C]∥InternationalConferenceonComputerandCommunicationEngineering.[S.l.]:IEEE,2008: 252-255.

[12] Mena-Chalco J P,Carrer H,Zana Y,et al.Identification of protein coding regions using the modified Gabor-wavelet transform[J].IEEE/ACMTransComputBiolBioinform,2008,5 (2):198-207.

[13] Rogic S,Mackworth A K,Ouellette B F.Evaluation of gene finding program onmammalian sequence[J].GenomicRes,2001,11 (5): 817-832.

[14] 徐尚蕾,饒妮妮.Bootstrap算法在基因預(yù)測(cè)中的閾值選取[D].成都:電子科技大學(xué),2011:19-21.