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

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

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

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

4、是目前生物信息學領域的一個研究熱點。二、數(shù)字序列映射與頻譜3-周期性：對給定的DNA序列，怎么去識別出其中的編碼序列（即外顯子），也稱為基因預測，是一個尚未完全解決的問題，也是當前生物信息學的一個最基礎、最首要的問題。基因預測問題的一類方法是基于統(tǒng)計學的1。很多國際生物數(shù)據(jù)網(wǎng)站上也有“基因識別”的算法。比如知名的數(shù)據(jù)網(wǎng)站/GENSCAN.html提供的基因識別軟件GENSCAN（由斯坦福大學研究人員研發(fā)的、可免費使用的基因預測軟件）,主要就是基于隱馬爾科夫鏈（HMM）方法。但是，它預測人的基因組中有45000個基因，相當于現(xiàn)在普遍認可數(shù)目的兩倍。另外，統(tǒng)計

5、預測方法通常需要將編碼序列信息已知的DNA序列作為訓練數(shù)據(jù)集來確定模型中的參數(shù)，從而提高模型的預測水平。但在對基因信息了解不多的情況下，基因識別的準確率會明顯下降。因此在目前基因預測研究中，采用信號處理與分析方法來發(fā)現(xiàn)基因編碼序列也受到廣泛重視 4。1. 數(shù)字序列映射在DNA序列研究中，首先需要把A、T、G、C四種核苷酸的符號序列，根據(jù)一定的規(guī)則映射成相應的數(shù)值序列，以便于對其作數(shù)字處理。令，長度（即核苷酸符號個數(shù)，又稱堿基對（Base Pair）長度，單位記為bp）為N的任意DNA序列，可表達為即A、T、G、C的符號序列：?，F(xiàn)對于任意確定的，令， 稱之為Voss映射5，于是生成相應的0-1序

6、列（即二進制序列）：， ()。例如，假設給定的一段DNA序列片段為S = ATCGTACTG，則所生成的四個0-1序列分別為：； ：；：； ：。這樣產(chǎn)生的四個數(shù)字序列又稱為DNA序列的指示序列(indicator Sequence)。2. 頻譜3-周期性為研究DNA編碼序列（外顯子）的特性，對指示序列分別做離散Fourier變換(DFT) (1)以此可得到四個長度均為N的復數(shù)序列，。計算每個復序列的平方功率譜，并相加則得到整個DNA序列的功率譜序列： (2)對于同一段DNA序列，其外顯子與內(nèi)含子序列片段的功率譜通常表現(xiàn)出不同的特性圖3 編號為BK006948.2的酵母基因DNA序列的功率譜（因

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

8、稱為DNA序列的“信噪比”（Signal Noise Ratio，SNR），即 (4)DNA序列的信噪比值的大小，既表示頻譜峰值(Peak Value)的相對高度，也反映編碼或非編碼序列3-周期性的強弱。信噪比大于某個適當選定的閾值（比如），是DNA序列上編碼序列片段（外顯子）通常滿足的特性，而內(nèi)含子則一般不具有該性質(zhì)6。在DNA序列, 中，若N為3的倍數(shù)，將核苷酸符號出現(xiàn)在該序列的0,3,6,. N3與1,4,7,N2以及2,5,8,N1等位置上的頻數(shù)分別記為和，則處的總功率譜值即為36 易見，當四種核苷酸符號（）在序列的上述第一、第二、第三個子序列上出現(xiàn)的頻數(shù)越接近相等時，處的譜值也就越接

9、近于零。所以，基因外顯子序列的功率譜曲線，在頻率處具有較大的頻譜峰值(Peak Value)，反映了在基因外顯子片段上，四種核苷酸符號在序列的三個子序列上分布的“非均衡性”。通常認為這種現(xiàn)象源于編碼基因序列“密碼子”（coden）使用的偏向性（bias）。雖然目前對此現(xiàn)象產(chǎn)生的“機理”還不是十分地清楚，但是頻譜的3-周期性被普遍認為是可用于識別基因編碼序列（外顯子）的一個重要的特征信息。3. 基因識別頻譜峰值特征的發(fā)現(xiàn)，或者頻譜與信噪比概念的引入，其最終目的是要探測、預報一個尚未被注釋的完整的DNA序列的所有基因編碼序列（外顯子）片段。預測結(jié)果DNA序列數(shù)值化映射DFT變換 功率譜或信噪比計算

10、外顯子判別分類閾值圖4 基于序列頻譜3周期性的的基因預測方法流程圖已經(jīng)有一些研究者提出了識別基因的算法（如參見6及其后面的文獻）。目前利用信噪比的基因識別算法通常有兩種：一是固定長度窗口滑動法2 3；另一是移動信噪比曲線識別法6。基于固定長度滑動窗口上頻譜曲線的基因識別方法: 對一個DNA序列和它的指示序列，。取長度M（通常取為3的倍數(shù)，例如M=99, 129, 255, 513等）作為固定窗口長度。對任意n（），在以n為中心的長度為M的序列片段n，n上（當n接近序列的兩端時，窗口實際有效長度可能會小于M），作四個指示序列的離散Fourier變換(DFT)并求出它在處總頻譜，即把這樣得到的頻譜

11、值，經(jīng)過標準化處理（即除以最大頻譜值），并畫出其頻譜曲線圖5 固定長度滑動窗口的頻譜曲線（人類線粒體基因，NC_012920_1.fasta）圖中紅色水平細線條是DNA序列實際的基因外顯子的區(qū)間?；瑒哟翱陬l譜曲線的峰與基因外顯子區(qū)間具有“對應”關系?；贒NA序列上“移動序列”信噪比曲線的基因識別方法：設已知DNA序列和它的指示序列，。對任意n（），通常n取3的倍數(shù)并逐漸增大。在n的左邊一個長度為n的序列片段0，n-1上，相應的子序列稱為DNA序列的“移動子序列”，作該移動子序列對應的四個指示序列的離散Fourier變換(DFT)并求出移動子序列，上的信噪比，其中為移動子序列的功率譜的平均值。

12、在坐標系中畫出移動序列的信噪比曲線（稱為信噪比移動曲線（SNR walk curve），見圖6）圖6 DNA移動序列其指示序列的信噪比曲線。（人類線粒體基因，NC_012920_1.fasta）圖中紅色水平細線條是DNA序列實際的基因外顯子的區(qū)間。DNA序列的信噪比移動曲線的峰、谷與基因外顯子區(qū)間的端點也具有較“明顯的”的對應關系。三、請研究的幾個問題：1. 功率譜與信噪比的快速算法對于很長的DNA序列，在計算其功率譜或信噪比時，離散Fourier變換(DFT)的總體計算量仍然很大，會影響到所設計的基因識別算法的效率。大家能否對Voss映射，探求功率譜與信噪比的某種快速計算方法？在基因識別研究

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

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

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

16、突變包括DNA序列中單個核苷酸的替換，刪除或者插入等。那么，能否利用頻譜或信噪比方法去發(fā)現(xiàn)基因編碼序列可能存在的突變呢？上面提出的基于頻譜3-周期性的基因預測四個方面問題中，“快速算法”與“閾值確定”是為設計基因預測算法做準備的。此外，在最后的延展性研究中，各隊也可以對你們自己認為有價值的其它相關問題展開探討。參考文獻： 【1】Burge, C., Karlin, S., 1997. Prediction of complete gene structures in human genomic DNA. J. Mol. Biol. 268, 7894. 【2】Anastassiou, D.,

17、2000. Frequency-domain analysis of biomolecular sequences. Bioinformatics 16, 10731081.【3】Kotlar, D., Lavner, Y., 2003. Gene prediction by spectral rotation measure: a new method for identifying protein-coding regions. Genome Res. 13, 19301937.【4】Berryman, M. J., Allison, A., 2005. Review of signal processing in genetics. Fluctuation and Noise Letters. 5(4), 13-35.【5】Sharma, S. D., Shakya，K., Sharma，S. N., 2011. Evaluation of DNA Mapping Schemes for Exon Detection