基因識別問題及其算法實現(xiàn)

1、基因識別問題及其算法實現(xiàn)4、背景介紹DNA是生物遺傳信息的載體，其化學(xué)名稱為脫氧核糖核酸（Deoxyribo nucleic acid，縮寫為DNA ）°DNA分子是一種長鏈聚合物，DNA序列由腺嘌呤（Ade nine, A），鳥嘌呤（Gua nine,G），胞嘧啶（Cytosine, C），胸腺嘧啶（Thymine, T）這四種核苷酸（nucleotide）符號按一定 的順序連接而成。其中帶有遺傳訊息的DNA片段稱為基因（Gene）（見圖1第一行）。其他的DNA序列片段，有些直接以自身構(gòu)造發(fā)揮作用，有些則參與調(diào)控遺傳訊息的表現(xiàn)。在真核生物的DNA序列中，基因通常被劃分為許多間隔的片

2、段 （見圖1第二行），其中編 碼蛋白質(zhì)的部分，即編碼序列（ Coding Sequenee）片段，稱為外顯子（Exon）,不編碼的部 分稱為內(nèi)含子（Intron ）。外顯子在DNA序列剪接（Splicing）后仍然會被保存下來，并可在圖1真核生物DNA序列（基因序列）結(jié)構(gòu)示意圖蛋白質(zhì)合成過程中被轉(zhuǎn)錄（transcription ）、復(fù)制（replication ）而合成為蛋白質(zhì)（見圖2）。DNA序列通過遺傳編碼來儲存信息，指導(dǎo)蛋白質(zhì)的合成，把遺傳信息準(zhǔn)確無誤地傳遞到蛋 白質(zhì)（protein）上去并實現(xiàn)各種生命功能。DNA序列蛋白質(zhì)序列圖2蛋白質(zhì)結(jié)構(gòu)示意圖對大量、復(fù)雜的基因序列的分析，傳統(tǒng)生物學(xué)

3、解決問題的方式是基于分子實驗的方法，其代價高昂。諾貝爾獎獲得者 W.吉爾伯特（Walter Gilbert，1932 ;【美】，第一個制備 出混合脫氧核糖核酸的科學(xué)家） 1991年曾經(jīng)指出： 現(xiàn)在，基于全部基因序列都將知曉，并以電子可操作的方式駐留在數(shù)據(jù)庫中，新的生物學(xué)研究模式的出發(fā)點應(yīng)是理論的。一個科學(xué)家將從理論推測出發(fā)， 然后再回到實驗中去， 追蹤或驗證這些理論假設(shè)。”隨著世界人類基 因組工程計劃的順利完成，通過物理或數(shù)學(xué)的方法從大量的DNA序列中獲取豐富的生物信息，對生物學(xué)、醫(yī)學(xué)、藥學(xué)等諸多方面都具有重要的理論意義和實際價值，也是目前生物信 息學(xué)領(lǐng)域的一個研究熱點 。1、數(shù)字序列映射與頻

4、譜 3周期性:對給定的DNA序列，怎么去識別出其中的編碼序列（即外顯子），也稱為基因預(yù)測， 是一個尚未完全解決的問題，也是當(dāng)前生物信息學(xué)的一個最基礎(chǔ)、最首要的問題?；蝾A(yù)測問題的一類方法是基于統(tǒng)計學(xué)的1。很多國際生物數(shù)據(jù)網(wǎng)站上也有“基因識別”的算法。比如知名的數(shù)據(jù)網(wǎng)站提供的基因識別軟件GENSCAN由斯坦福大學(xué)研究人員研發(fā)的、可免費使用的基因預(yù)測軟件），主要就是基于隱馬爾科夫鏈（HMM方法。但是，它預(yù)測人的基因組中有45000個基因，相當(dāng)于現(xiàn)在普遍認(rèn)可數(shù)目的兩倍。另外，統(tǒng)計預(yù)測方法通常需要將編碼序列信息已知的DNA序列作為訓(xùn)練數(shù)據(jù)集來確定模型中的參數(shù)，從而提高模型的預(yù)測水平。但在對基因信息了解

5、不多的情況下，基因識別的準(zhǔn)確率會明顯下降。因此在目前基因預(yù)測研究中，采用信號處理與分析方法來發(fā)現(xiàn)基因編碼序列也受到廣泛 重視。1.數(shù)字序列映射在DNA序列研究中，首先需要把 A、T、G、C四種核苷酸的符號序列，根據(jù)一定的規(guī) 則映射成相應(yīng)的數(shù)值序列，以便于對其作數(shù)字處理。令I(lǐng) =A,T,G,C，長度（即核苷酸符號個數(shù)，又稱堿基對（ Base Pair ）長度，單位 記為bp ）為N的任意DNA序列，可表達(dá)為S = Sn | Sn I, n =0,1,2,|IN -1即A、T、G、C的符號序列S : S0, S1,H|,SN -1?，F(xiàn)對于任意確定的 b I，令1Ub滬。,Sn二 bSn 嚴(yán)bn =

6、0,1,2,111 N -1稱之為Voss映射5，于是生成相應(yīng)的0-1 序列（即二進制序列）ubn : ub0, ub1J|l,，ubN -1 (b l)。例如，假設(shè)給定的一段DNA序列片段為S= ATCGTACTG，則所生成的四個0-1序列分別為：uAn : 1,0,0,0,0,1,0,0,0 ；Ugn : 0,0,0,1,0,0,0,0,1 ；ucn : 0,0,1,0,0,0,1,0,0 ；比門 : 0,1,0,0,1,0,0,1,0。這樣產(chǎn)生的四個數(shù)字序列又稱為DNA序列的指示序列(indicator Sequenee)。2.頻譜3-周期性為研究DNA編碼序列(外顯子)的特性，對指示序

7、列分別做離散Fourier變換(DFT)N -1.2-nkUbk = E %門, k = 0,1|,N-1(1)n=0以此可得到四個長度均為n的復(fù)數(shù)序列ubk，b I。計算每個復(fù)序列Ubk的平方功率譜，并相加則得到整個DNA序列S的功率譜序列Pk:2 2 2 2Pk = UA【k +UT【k +UG【k +Uck , k=0,1,川N 1(2)對于同一段DNA序列，其外顯子與內(nèi)含子序列片段的功率譜通常表現(xiàn)出不同的特性10000500050060000100200300400k10000500000100200300400500600k圖3編號為BK006948.2的酵母基因DNA序列的功率譜(

8、因為對稱性，實際這里只給出了功率譜圖 的一半)。(a)上圖是基因上一段外顯子(區(qū)間為81787，82920，長1134bp)對應(yīng)的指示序列映射的功率 譜，它具有3-周期性；(b)下圖是基因上一段內(nèi)含子(區(qū)間為96361,97551，長1191bp)的指示序列的功率譜，它不具有3-周期性?？梢钥吹剑和怙@子序列的功率譜曲線在頻率k=N處，具有較大的頻譜峰值(Peak3Value)，而內(nèi)含子則沒有類似的峰值。這種統(tǒng)計現(xiàn)象被稱為堿基的3-周期(3-base Periodicity)23記DNA序列S的總功率譜的平均值為' Pkk _0E =N而將DNA序列在特定位置，即 k處的功率譜值，與整個

9、序列S的總功率譜的平均值的3比率稱為DNA序列的“信噪比” (Sig nal Noise Ratio，SNR),即DNA序列的信噪比值的大小，既表示頻譜峰值(Peak Value)的相對高度，也反映編碼或非編碼序列3-周期性的強弱。信噪比R大于某個適當(dāng)選定的閾值 Rq (比如Ro = 2 )，是DNA序列上編碼序列片段(外顯子)通常滿足的特性，而內(nèi)含子則一般不具有該性質(zhì)6。在DNA序列Sn, n =0,1,2川IN -1中，若N為3的倍數(shù)，將核苷酸符號 b I =A,T,G,C出現(xiàn)在該序列的0,3,6,N 3與1,4,7,N 2以及2,5,8,N 1等位置上的頻數(shù)分別記為x.,yb和，則N處的

10、總功率譜值即為363P3N2N 4Nj 2 兀'a2NJJ2 叫=2u b=E遲 Ubn eN=z遲 Ubn e 3b旺3b社nb令2八 (xb yb Z： - Xbyb - XbZb - ybZb) b三122 兀j2_nyb e 3Zb e 3b.I易見，當(dāng)四種核苷酸符號b ( b I )在序列的上述第一、第二、第三個子序列上出現(xiàn)的頻數(shù)xb, yb, zb越接近相等時，處的譜值也就越接近于零。所以，基因外顯子序列的功率譜3曲線，在頻率處具有較大的頻譜峰值 (Peak Value),反映了在基因外顯子片段上，四種核3苷酸符號在序列的三個子序列上分布的“非均衡性”。通常認(rèn)為這種現(xiàn)象源于

11、編碼基因序列“密碼子” (code n)使用的偏向性(bias)。雖然目前對此現(xiàn)象產(chǎn)生的“機理”還不是十分地清楚，但是頻譜的3-周期性被普遍認(rèn)為是可用于識別基因編碼序列(外顯子)的一個重要的特征信息。3.基因識別頻譜峰值特征的發(fā)現(xiàn)，或者頻譜與信噪比概念的引入，其最終目的是要探測、預(yù)報一個尚未被注釋的完整的 DNA序列的所有基因編碼序列（外顯子）片段。外顯子 判別分類預(yù)測結(jié)果已經(jīng)有一些研究者提出了識別基因的算法（如參見及其后面的文獻）。目前利用信噪比的基因識別算法通常有兩種：一是固定長度窗口滑動法2 3;另一是移動信噪比曲線識別法6?；诠潭ㄩL度滑動窗口上頻譜曲線的基因識別方法對一個DNA序列S

12、和它的指示序列 ub n , I，n =0,1,2,川N-1。取長度M （通 常取為3的倍數(shù)，例如 M=99, 129, 255, 513等）作為固定窗口長度。對任意（0 n乞N -1 ），在以n為中心的長度為 M的序列片段n 吐7 n 也2 2上（當(dāng)n接近序列的兩端時，窗口實際有效長度可能會小于M），作四個指示序列的離散Fourier 變換（DFT）-4i 才 T_.2二ikUbk二 'ubieF,k = 0,1川,M-1M -4并求出它在 處總頻譜p（n；M），即33M2M2M2M2 A咗+7+Ug弓+UC弓=P（n；M3）把這樣得到的頻譜值號）川N,，經(jīng)過標(biāo)準(zhǔn)化處理（即除以最大頻

13、譜值0mmaxp（n;Mh），并畫出其頻譜曲線030000.135004000450050005500nucleotide position n6000650070009 87 6 5 o o O4 3 2 o o O.-cpmd 莒圖5固定長度滑動窗口的頻譜 p = p（n； §）曲線（人類線粒體基因，NC_012920_1.fasta）圖中紅色水平細(xì)線條是DNA序列實際的基因外顯子的區(qū)間。滑動窗口頻譜 p（n；M）曲線的3峰與基因外顯子區(qū)間具有“對應(yīng)”關(guān)系?；贒NA序列上“移動序列”信噪比曲線的基因識別方法:設(shè)已知DNA序列S和它的指示序列ubn，I，n =0,1,2,l|lN

14、 -1。對任意n（0 ： n込N -1），通常n取3的倍數(shù)并逐漸增大。 在n的左邊一個長度為 n的序列片段0， n-1上，相應(yīng)的子序列 Son稱為DNA序列S的“移動子序列”，作該移動子序列對應(yīng)的四個指示序列的離散 Fourier變換（DFT）j2 二ikUbk=E Ubiew,k = 0,1,川,n 1i =0并求出移動子序列 S02，n =0,1,Hl,N-1上的信噪比RngQUc?0 ： n < N -1Enn-1其中En為移動子序列 S0nj的功率譜的平均值_ 送 P kEn=。在坐標(biāo)系中畫出移動序n列Son4的信噪比曲線 Rn（稱為信噪比移動曲線（SNR walk curve）

15、，見圖6）n ucleotide positi on n圖6 DNA移動序列其指示序列的信噪比曲線。（人類線粒體基因，NC_012920_1.fasta）圖中紅色水平細(xì)線條是 DNA序列實際的基因外顯子的區(qū)間。DNA序列的信噪比移動曲線的峰、谷與基因外顯子區(qū)間的端點也具有較“明顯的”的對應(yīng)關(guān)系。三、請研究的幾個問題:1. 功率譜與信噪比的快速算法對于很長的DNA序列，在計算其功率譜或信噪比時， 離散Fourier變換（DFT）的總體計 算量仍然很大，會影響到所設(shè)計的基因識別算法的效率。大家能否對Voss映射，探求功率譜與信噪比的某種快速計算方法？在基因識別研究中，為了通過引入更好的數(shù)值映射而獲

16、取DNA序列更多的信息，除了上面介紹的 Voss映射外，實際上人們還研究過許多不同的數(shù)值映射方法。例如，著名的 Z-curve映射（參見5或者附件1 ）。試探討Z-curve映射的頻譜與信噪比和 Voss映射下的頻 譜與信噪比之間的關(guān)系；此外，能否對實數(shù)映射，如：A > 0,C > 1,G > 2,T >3，也給出功率譜與信噪比的快速計算公式？2. 對不同物種類型基因的閾值確定對特定的基因類型的 DNA序列，將其信噪比R的判別閾值取為 =2，帶有一定的主 觀性、經(jīng)驗性。對不同的基因類型， 所選取的判別閾值也許應(yīng)該是不同的。附件中給出了來自于著名的生物數(shù)據(jù)網(wǎng)站： http

17、:/www. ncbi. nl m /guide/ 的幾個基因序列數(shù)據(jù)，另外也 給出了帶有編碼外顯子信息的 100個人和鼠類的，以及200個哺乳動物類的基因序列的樣本 數(shù)據(jù)集合。大家還可以從生物數(shù)據(jù)庫下載更多的數(shù)據(jù)，找你們認(rèn)為具有代表性的基因序列， 并對每類基因研究其閾值確定方法和閾值結(jié)果。 此外， 對按照頻譜或信噪比特征將編碼與非 編碼區(qū)間分類的有效性，以及分類識別時所產(chǎn)生的分類錯誤作適當(dāng)分析。3. 基因識別算法的實現(xiàn)我們的目的是要探測、預(yù)報尚未被注釋的、完整的 DNA 序列的所有基因編碼序列（外 顯子）。目前基因識別方面的多數(shù)算法結(jié)果還不是很充分。例如前面所列舉的某些基因識別

18、 算法，由于 DNA 序列隨機噪聲的影響等原因，還很難“精確地”確定基因外顯子區(qū)間的兩 個端點。對此，你的建模團隊有沒有更好的解決方法？請對你們所設(shè)計的基因識別算法的準(zhǔn)確率做出適當(dāng)評估，并將算法用于對附件中給出的6個未被注釋的 DNA序列（gene6）的編碼區(qū)域的預(yù)測。4. 延展性研究 在基因識別研究中，還有很多問題有待深入探討。比如（ 1）采用頻譜或信噪比這樣單一的判別特征，也許是影響、限制基因識別正確率的一 個重要原因。 人們發(fā)現(xiàn)， 對某些 DNA 序列而言， 其部分編碼序列 （外顯子），尤其是短的 （長 度小于100bp）的編碼序列，就可能不具有頻譜或者信噪比顯著性。你們團隊能否總結(jié)，甚

19、 至獨自提出一些識別基因編碼序列的其它特征指數(shù)，并對此做相關(guān)的分析？（ 2）“基因突變 ”是生物醫(yī)學(xué)等方面的一個關(guān)注熱點。基因突變包括 DNA 序列中單個核 苷酸的替換， 刪除或者插入等。 那么， 能否利用頻譜或信噪比方法去發(fā)現(xiàn)基因編碼序列可能 存在的突變呢？上面提出的基于頻譜 3-周期性的基因預(yù)測四個方面問題中，“快速算法” 與“閾值確定”是為設(shè)計基因預(yù)測算法做準(zhǔn)備的。 此外， 在最后的延展性研究中， 各隊也可以對你們自己認(rèn) 為有價值的其它相關(guān)問題展開探討。參考文獻：【1 】Burge, C., Karlin, S., 1997. Prediction of complete gene structures in human genomic DNA.J. Mol. Biol. 268, 78 -94.【2】Anastassiou, D., 2000. Frequency-domain analysis of biomolecular sequences.Bioi nformatics 16, 1073 081.【3】Kotlar, D.,