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

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

2、許多間隔的片段（見圖 1第二行），其中編碼蛋白質(zhì)的部分，即編碼序列（Coding Sequence）片段，稱為外顯子（Exon），不編碼的部分稱為內(nèi)含子（Intron ）。外顯子在 DNA 序列剪接（ Splicing）后仍然會(huì)被保存下來，并可在基因（ Gene）DNA 序列外顯子 (Exon)內(nèi)含子 (Intron)圖 1 真核生物 DNA 序列（基因序列）結(jié)構(gòu)示意圖蛋白質(zhì)合成過程中被轉(zhuǎn)錄（transcription ）、復(fù)制（ replication ）而合成為蛋白質(zhì)（見圖2）。DNA 序列通過遺傳編碼來儲(chǔ)存信息，指導(dǎo)蛋白質(zhì)的合成，把遺傳信息準(zhǔn)確無誤地傳遞到蛋白質(zhì)（ protein ）上去

3、并實(shí)現(xiàn)各種生命功能?；颍?Gene）DNA 序列剪接、轉(zhuǎn)錄、復(fù)制蛋白質(zhì)序列圖 2 蛋白質(zhì)結(jié)構(gòu)示意圖對(duì)大量、復(fù)雜的基因序列的分析，傳統(tǒng)生物學(xué)解決問題的方式是基于分子實(shí)驗(yàn)的方法，其代價(jià)高昂。諾貝爾獎(jiǎng)獲得者W. 吉爾伯特（ Walter Gilbert ， 1932 ；【美】，第一個(gè)制備出混合脫氧核糖核酸的科學(xué)家）1991 年曾經(jīng)指出： “現(xiàn)在，基于全部基因序列都將知曉，并以電子可操作的方式駐留在數(shù)據(jù)庫中，新的生物學(xué)研究模式的出發(fā)點(diǎn)應(yīng)是理論的。一個(gè)科學(xué)家將從理論推測(cè)出發(fā)，然后再回到實(shí)驗(yàn)中去，追蹤或驗(yàn)證這些理論假設(shè)。” 隨著世界人類基因組工程計(jì)劃的順利完成，通過物理或數(shù)學(xué)的方法從大量的DNA 序列中

4、獲取豐富的生物信1息，對(duì)生物學(xué)、 醫(yī)學(xué)、 藥學(xué)等諸多方面都具有重要的理論意義和實(shí)際價(jià)值，也是目前生物信息學(xué)領(lǐng)域的一個(gè)研究熱點(diǎn)。二、數(shù)字序列映射與頻譜3- 周期性：對(duì)給定的 DNA 序列， 怎么去識(shí)別出其中的編碼序列（即外顯子），也稱為基因預(yù)測(cè)，是一個(gè)尚未完全解決的問題，也是當(dāng)前生物信息學(xué)的一個(gè)最基礎(chǔ)、最首要的問題。基因預(yù)測(cè)問題的一類方法是基于統(tǒng)計(jì)學(xué)的1 。很多國際生物數(shù)據(jù)網(wǎng)站上也有“基因識(shí)別”的算法。比如知名的數(shù)據(jù)網(wǎng)站提供的基因識(shí)別軟件GENSCAN（由斯坦福大學(xué)研究人員研發(fā)的、可免費(fèi)使用的基因預(yù)測(cè)軟件）, 主要就是基于隱馬爾科夫鏈（ HMM）方法。但是，它預(yù)測(cè)人的基因組中有45000 個(gè)基因

5、，相當(dāng)于現(xiàn)在普遍認(rèn)可數(shù)目的兩倍。另外，統(tǒng)計(jì)預(yù)測(cè)方法通常需要將編碼序列信息已知的DNA 序列作為訓(xùn)練數(shù)據(jù)集來確定模型中的參數(shù)，從而提高模型的預(yù)測(cè)水平。但在對(duì)基因信息了解不多的情況下，基因識(shí)別的準(zhǔn)確率會(huì)明顯下降。因此在目前基因預(yù)測(cè)研究中，采用信號(hào)處理與分析方法來發(fā)現(xiàn)基因編碼序列也受到廣泛重視 4 。1. 數(shù)字序列映射在 DNA 序列研究中，首先需要把 A、 T、 G、C 四種核苷酸的符號(hào)序列，根據(jù)一定的規(guī)則映射成相應(yīng)的數(shù)值序列，以便于對(duì)其作數(shù)字處理。令 I A,T, G,C ，長度（即核苷酸符號(hào)個(gè)數(shù)，又稱堿基對(duì)（Base Pair ）長度，單位記為 bp）為 N 的任意 DNA 序列，可表達(dá)為S

6、Sn| SnI , n0,1,2, N1即 A、 T、G、C 的符號(hào)序列 S ： S0, S1,SN1 ?，F(xiàn)對(duì)于任意確定的 bI ，令1,S nbN 1ub nS n， n 0,1,2,0,b稱之為 Voss映射 5 ，于是生成相應(yīng)的0-1 序列（即二進(jìn)制序列 ） ub n ： ub 0,ub1, , ，ub N1 ( b I )。例如，假設(shè)給定的一段 DNA 序列片段為 S = ATCGTACTG，則所生成的四個(gè)0-1 序列分別為： uA n ： 1,0,0,0,0,1,0,0,0 ；uG n ： 0,0,0,1,0,0,0,0,1 ；2 uC n ： 0,0,1,0,0,0,1,0,0；

7、uT n ： 0,1,0,0,1,0,0,1,0。這樣產(chǎn)生的四個(gè)數(shù)字序列又稱為DNA 序列的指示序列(indicator Sequence)。2. 頻譜 3-周期性為研究 DNA 編碼序列（外顯子）的特性，對(duì)指示序列分別做離散Fourier 變換 (DFT )N 1j2nkU ku neN , k0,1, ,N 1(1)bbn 0以此可得到四個(gè)長度均為N 的復(fù)數(shù)序列 U b k ，bI 。計(jì)算每個(gè)復(fù)序列 U b k 的平方功率譜，并相加則得到整個(gè)DNA 序列 S 的功率譜序列 P k ：2222P k U A k U T k U G k U C k, k 0,1, N 1(2)對(duì)于同一段DNA

8、 序列，其外顯子與內(nèi)含子序列片段的功率譜通常表現(xiàn)出不同的特性)k(P10000500000100200300400500600k10000)k(P500000100200300400500600k圖 3 編號(hào)為 BK006948.2 的酵母基因 DNA 序列的功率譜（ 因?yàn)閷?duì)稱性，實(shí)際這里只給出了功率譜圖的一半 ）。 ( a) 上圖是基因上一段外顯子 (區(qū)間為 81787 ，82920，長 1134bp) 對(duì)應(yīng)的指示序列映射的功率譜，它具有 3- 周期性； (b) 下圖是基因上一段內(nèi)含子（區(qū)間為96361,97551 ，長 1191bp）的指示序列的功率譜，它不具有 3- 周期性。可以看到：外

9、顯子序列的功率譜曲線在頻率kN處，具有較大的頻譜峰值(Peak3Value)，而內(nèi)含子則沒有類似的峰值。這種統(tǒng)計(jì)現(xiàn)象被稱為堿基的3-周期 ( 3-base Periodicity )23 。記 DNA 序列 S 的總功率譜的平均值為N 1P k Ek 0(3)N3而將 DNA 序列在特定位置，即 kNS 的總功率譜的平均值的處的功率譜值，與整個(gè)序列3比率稱為 DNA 序列的“信噪比” （ Signal Noise Ratio ，SNR），即P NR3(4)EDNA 序列的信噪比值的大小，既表示頻譜峰值(Peak Value)的相對(duì)高度，也反映編碼或非編碼序列 3-周期性的強(qiáng)弱。信噪比 R 大于

10、某個(gè)適當(dāng)選定的閾值R0（比如 R02 ），是 DNA 序列上編碼序列片段 （外顯子）通常滿足的特性，而內(nèi)含子則一般不具有該性質(zhì)6 。在 DNA序列 S n , n 0,1,2,N 1 中，若 N為 3的倍數(shù)，將核苷酸符號(hào)b I A,T , G ,C 出現(xiàn)在該序列的 0,3,6,. N 3與1,4,7,N 2以及 2,5,8, N 1等位置上的頻數(shù)分別記為xb, yb 和 zb ，則 N 處的總功率譜值即為3632N22N2N 1nN 12PNj3jnU bb I n 0 ub n eNb I n 0 ub n e33b I3j 2j 22(xb2yb2zb2xb yb xb zb yb zb

11、)xb yb e 3zb e 3b Ib I易見，當(dāng)四種核苷酸符號(hào)b （ bI ）在序列的上述第一、第二、第三個(gè)子序列上出現(xiàn)的頻數(shù) xb , yb , zb 越接近相等時(shí)，N 處的譜值也就越接近于零。所以，基因外顯子序列的功率譜3曲線，在N 頻率處具有較大的頻譜峰值 (Peak Value)，反映了在基因外顯子片段上，四種核3苷酸符號(hào)在序列的三個(gè)子序列上分布的“非均衡性”。通常認(rèn)為這種現(xiàn)象源于編碼基因序列“密碼子”（ coden）使用的偏向性（bias）。雖然目前對(duì)此現(xiàn)象產(chǎn)生的“機(jī)理”還不是十分地清楚，但是頻譜的3-周期性被普遍認(rèn)為是可用于識(shí)別基因編碼序列（外顯子）的一個(gè)重要的特征信息。3.

12、基因識(shí)別頻譜峰值特征的發(fā)現(xiàn)，或者頻譜與信噪比概念的引入，其最終目的是要探測(cè)、預(yù)報(bào)一個(gè)尚未被注釋的完整的DNA 序列的所有基因編碼序列（外顯子）片段。4閾值DNA 序列數(shù)值化DFT功率譜或外顯子預(yù)測(cè)結(jié)果映射變換信噪比計(jì)算判別分類圖 4基于序列頻譜 3周期性的的基因預(yù)測(cè)方法流程圖已經(jīng)有一些研究者提出了識(shí)別基因的算法（如參見6 及其后面的文獻(xiàn)） 。目前利用信噪比的基因識(shí)別算法通常有兩種：一是固定長度窗口滑動(dòng)法2 3；另一是移動(dòng)信噪比曲線識(shí)別法6 ?；诠潭ㄩL度滑動(dòng)窗口上頻譜曲線的基因識(shí)別方法:對(duì)一個(gè) DNA 序列S和它的指示序列 u n ，1。取長度（通bb I n 0,1,2, NM常取為 3 的

13、倍數(shù)，例如M=99, 129, 255, 513 等）作為固定窗口長度。對(duì)任意 n（ 0 nN 1 ），在以 n 為中心的長度為M 的序列片段 nM 1M 1，n22上（當(dāng) n 接近序列的兩端時(shí)，窗口實(shí)際有效長度可能會(huì)小于M ），作四個(gè)指示序列的離散Fourier 變換 (DFT )inUb kinM 12M 12j 2 ikub ieM,k0,1, M1并求出它在 M 處總頻譜(;)，即M3p n3MM2M2M2M2MP3U A 3UT 3U G 3U C 3p(n;3)把這樣得到的頻譜值( ; M)， n0,1,2,N1，經(jīng)過標(biāo)準(zhǔn)化處理（即除以最大頻譜值p n3M)），并畫出其頻譜曲線ma

14、x p(n;0 n N 13)n(pmutrcepsAND0350040004500500055006000650070003000nucleotide position n5M圖 5 固定長度滑動(dòng)窗口的頻譜p p n( ; ) 曲線（人類線粒體基因，NC_012920_1.fasta ）3M圖中紅色水平細(xì)線條是DNA 序列實(shí)際的基因外顯子的區(qū)間?；瑒?dòng)窗口頻譜( ; )曲線的p n3峰與基因外顯子區(qū)間具有“對(duì)應(yīng)”關(guān)系?；?DNA 序列上“移動(dòng)序列”信噪比曲線的基因識(shí)別方法：設(shè)已知 DNA 序列 S 和它的指示序列 ub n， bI ，

15、n0,1,2,N 1。對(duì)任意 n（ 0 n N 1），通常 n 取 3 的倍數(shù)并逐漸增大。 在 n 的左邊一個(gè)長度為n 的序列片段 0，n-1 上，相應(yīng)的子序列 S0 n1 稱為 DNA 序列 S 的“ 移動(dòng)子序列 ”，作該移動(dòng)子序列對(duì)應(yīng)的四個(gè)指示序列的離散Fourier變換 (DFT )i n1j2 ikU b kub i eM,k0,1, n 1i 0并求出移動(dòng)子序列S0 n1 ， n0,1, N1 上的信噪比 R n2222P nU A nU T n U G nU C nR n33333， 0n N 1E nE nn1其中 En 為移動(dòng)子序列P kS0 n1 的功率譜的平均值E nk0。

16、在坐標(biāo)系中畫出移動(dòng)序n列 S0 n 1 的信噪比曲線R n （稱為信噪比移動(dòng)曲線（ SNR walk curve ），見圖 6）141210)n8(RRNS6420300035004000450050005500600065007000nucleotide position n圖 6 DNA 移動(dòng)序列其指示序列的信噪比曲線。（人類線粒體基因， NC_012920_1.fasta）圖中紅色水平細(xì)線條是DNA 序列實(shí)際的基因外顯子的區(qū)間。DNA 序列的信噪比移動(dòng)曲線6的峰、谷與基因外顯子區(qū)間的端點(diǎn)也具有較“明顯的”的對(duì)應(yīng)關(guān)系。三、請(qǐng)研究的幾個(gè)問題：1. 功率譜與信噪比的快速算法對(duì)于很長的DNA 序

17、列， 在計(jì)算其功率譜或信噪比時(shí)，離散 Fourier 變換 (DFT )的總體計(jì)算量仍然很大，會(huì)影響到所設(shè)計(jì)的基因識(shí)別算法的效率。大家能否對(duì)Voss 映射，探求功率譜與信噪比的某種快速計(jì)算方法？在基因識(shí)別研究中，為了通過引入更好的數(shù)值映射而獲取DNA 序列更多的信息，除了上面介紹的Voss 映射外，實(shí)際上人們還研究過許多不同的數(shù)值映射方法。例如，著名的Z-curve 映射（參見 5 或者附件1）。試探討 Z-curve 映射的頻譜與信噪比和Voss 映射下的頻譜與信噪比之間的關(guān)系；此外，能否對(duì)實(shí)數(shù)映射，如：A0, C1, G2, T3 ，也給出功率譜與信噪比的快速計(jì)算公式？2.對(duì)不同物種類型基

18、因的閾值確定對(duì)特定的基因類型的 DNA 序列，將其信噪比 R 的判別閾值取為 R02 ，帶有一定的主觀性、 經(jīng)驗(yàn)性。 對(duì)不同的基因類型， 所選取的判別閾值也許應(yīng)該是不同的。附件中給出了來自于著名的生物數(shù)據(jù)網(wǎng)站：/guide/的幾個(gè)基因序列數(shù)據(jù)， 另外也給出了帶有編碼外顯子信息的100 個(gè)人和鼠類的， 以及 200個(gè)哺乳動(dòng)物類的基因序列的樣本數(shù)據(jù)集合。大家還可以從生物數(shù)據(jù)庫下載更多的數(shù)據(jù)，找你們認(rèn)為具有代表性的基因序列，并對(duì)每類基因研究其閾值確定方法和閾值結(jié)果。此外，對(duì)按照頻譜或信噪比特征將編碼與非編碼區(qū)間分類的有效性，以及分類識(shí)別時(shí)所產(chǎn)生的分

19、類錯(cuò)誤作適當(dāng)分析。3. 基因識(shí)別算法的實(shí)現(xiàn)我們的目的是要探測(cè)、預(yù)報(bào)尚未被注釋的、完整的DNA 序列的所有基因編碼序列（外顯子）。目前基因識(shí)別方面的多數(shù)算法結(jié)果還不是很充分。例如前面所列舉的某些基因識(shí)別算法，由于DNA 序列隨機(jī)噪聲的影響等原因，還很難“精確地”確定基因外顯子區(qū)間的兩個(gè)端點(diǎn)。對(duì)此，你的建模團(tuán)隊(duì)有沒有更好的解決方法？請(qǐng)對(duì)你們所設(shè)計(jì)的基因識(shí)別算法的準(zhǔn)確率做出適當(dāng)評(píng)估，并將算法用于對(duì)附件中給出的6 個(gè)未被注釋的DNA 序列（ gene6）的編碼區(qū)域的預(yù)測(cè)。74. 延展性研究在基因識(shí)別研究中，還有很多問題有待深入探討。比如（ 1）采用頻譜或信噪比這樣單一的判別特征，也許是影響、限制基因識(shí)

20、別正確率的一個(gè)重要原因。 人們發(fā)現(xiàn)， 對(duì)某些 DNA 序列而言， 其部分編碼序列 （外顯子），尤其是短的 （長度小于 100bp）的編碼序列，就可能不具有頻譜或者信噪比顯著性。你們團(tuán)隊(duì)能否總結(jié)，甚至獨(dú)自提出一些識(shí)別基因編碼序列的其它特征指數(shù)，并對(duì)此做相關(guān)的分析？（ 2）“基因突變 ”是生物醫(yī)學(xué)等方面的一個(gè)關(guān)注熱點(diǎn)?；蛲蛔儼―NA 序列中單個(gè)核苷酸的替換， 刪除或 者插入等。 那么，能否利用頻譜或信噪比方法去發(fā)現(xiàn)基因編碼序列可能存在的突變呢？上面提出的基于頻譜3-周期性的基因預(yù)測(cè)四個(gè)方面問題中，“快速算法” 與“閾值確定”是為設(shè)計(jì)基因預(yù)測(cè)算法做準(zhǔn)備的。此外， 在最后的延展性研究中，各隊(duì)也可以對(duì)你們自己認(rèn)為有價(jià)值的其它相關(guān)問題展開探討。參考文獻(xià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 se