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

1、基因識別問題及其算法實現(xiàn)一、背景介紹dna是生物遺傳信息的載體，只化學名稱為脫氧核糖核酸(deoxyribonucleic acid,縮 寫為dna)cdna分子是一種長鏈聚合物,dna序列由腺喋吟adenine, 鳥瞟吟(guanine, g),胞i密噪(cytosine, c),胸腺n密唳thymine, t)這四種核昔酸(nucleotide)符號按一定 的順序連接1何成。其中帶有遺傳訊息的dna片段稱為基因(gm>)(見圖1第一行)。其他的 dna序列片段，冇些直接以自身構造發(fā)揮作川，冇些則參與調控遺傳訊息的表現(xiàn)。在真核生物的dna序列屮，基因通常被劃分為許多間隔的片段(見圖1第

2、二行)，其屮編碼蛋白質的部分，即編碼序列(coding sequence)片段，稱為外顯子(exon),不編碼的部 分稱為內含子(buron)。外顯子在dva序列剪接(splicing)后仍然會被保存下來，并可在外顯子(exon) 內含子(intron )圖1真核牛物dna用列(基因序列)結構示意圖蛋白質合成過程中被轉錄(transcription)復制(replication)而合成為蛋白質(見圖2)。dna序列通過遺傳編碼來儲存信息，指導蛋白質的合成，把遺傳信息準確無誤地傳遞到蛋白質(protein)上去并實現(xiàn)各種生命功能。dna序列蛋白質序列圖2蛋白質結構示意圖對大最、復雜的基因序列的分

3、析，傳統(tǒng)生物學解決問題的方式是基于分子實驗的方法, 其代價高昂。諾貝爾獎獲得者w.吉爾伯特(waltergilbert, 1932；【美】,第一個制備 出混合脫氧核糖核酸的科學家)1991年曾經指出：“現(xiàn)在，基于全部基因序列都將知曉，并 以電子可操作的方式駐超在數(shù)據(jù)庫中，新的生物學研究模式的出發(fā)點應是理論的。一個科學 家將從理論推測出發(fā)，然后再回到實驗中去，追蹤或驗證這些理論假設?！彪S著世界人類基 因組工程計劃的順利完成，通過物理或數(shù)學的方法從人最的dna序列屮獲取豐富的生物信息，對生物學、醫(yī)學、藥學等諸多方而都具有重要的理論意義和實際價值，也是h前生物信 息學領域的-個研究熱點。二、數(shù)字序列

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

5、已知的dna序列作為訓練數(shù)據(jù) 集來確定模型中的參數(shù)，從而提高模型的預測水平。但在對基因信息了解不多的悄況下，基 因識別的準確率會明顯下降。因此在目前基因預測研究屮，采用信號處理與分析方法來發(fā)現(xiàn)基因編碼序列也受到廣泛 重視。1. 數(shù)字序列映射在lwa序列研究中，首先需要把a、7 g、c四種核昔酸的符號序列，根據(jù)一定的規(guī) 則映射成相應的數(shù)值序列，以便于對其作數(shù)字處理。令z = a,t,g,c,長度（即核廿酸符號個數(shù)，又稱堿基對（base paid長度，單位 記為切）為n的任意qm4序列，可表達為s = sn i s“， = 0丄2,n-1即a、t、g、c的符號序列s： s0,su,sw l?，F(xiàn)對于

6、任意確定的bwl ,令sn = bsn工 b稱之為vb刖映射，于是生成相應的01序列（即二進制序列）uhn：地0,旳j1,，,uhn -i (b e i )o例如，假設給定的一段dma序列片段為s = atcgtactg,則所牛成的四個0-1序列分 別為:uan: 1,0,0,0,0,1,0,0,0 ；%“： (),0,0丄0,0,0,0,1；ucn: 0,0,1,0,0,0,1,0,0 ；utn: 0,1,0,0,1,0,0,1,0。這樣產生的四個數(shù)7序列乂稱為dna序列的指示序列indicator sequence)。2. 頻譜3周期性為研究dn4編碼序列(外顯子)的特性，對指示序列分別做

7、離散fourier變換(dft)n-1_ 2兀nk5幻=工坷川幺7 n ,a=0,l,n-l(1)n=0以此可得到四個長度均為n的復數(shù)序列uhk, bel.計算每個復序列uhk的平方功率譜，并相加則得到整個dna序列s的功率譜序列pr：pk = ua*|2 + utkf + ugkf 4-uckf, k = o,hn-(2)對于同一段dm4序列，其外顯子與內含子序列片段的功率譜通常表現(xiàn)出不同的特性圖3編號為bk006948.2的酵母基因dva序列的功率譜(因為對稱性，實際這里只給出了功率譜圖 的一半)。 上圖是基因上一段外顯了(區(qū)間為81787, 82920,長1134切)對應的指示序列映射的

8、功率 譜,它具有3周期性;(b)下圖是基因上一段內含子(區(qū)間為96361,97551，長1191切)的指示序列的功 率譜，它不具有3-周期性。n可以看到：外顯子序列的功率譜曲線在頻率k=處，具有較人的頻譜峰值3value),而內含了則沒冇類似的峰值。這種統(tǒng)計現(xiàn)彖被稱為堿壟的3周期(3-base periodicity) 3o記dna序列s的總功率譜的平均值為n-1_工p幻e =旦(3)而將dna序列在特定位置，即k =n二丄處的功率譜值，與整個序列s的總功率譜的平均值的3e比率稱為dna序列的“信噪比” (signalnoise ratio, snr),即dna序列的信噪比值的大小，既表示頻譜

9、峰值(peak value)的相對高度，也反映編碼或 非編碼序列3周期性的強弱。信噪比/?大于某個適當選定的閾值他(比如/?0=2), ina序列上編碼序列片段(外 顯了)通常滿足的特性，而內含了則一般不具冇該性質叫在dm4序列sn, n=0,l,2,n 1中，若n為3的倍數(shù)，將核廿酸符號b ei = a,t,g,c出現(xiàn)在該序列的 0,3,6,. n-3 與 l,4,7,.n-2 以及 2,5,&.n-1 等位置上n的頻數(shù)分別記為必和z方，則二處的總功率譜值即為c n 2/r/rn-1_ : *n只寧=工匕卻=z為訕“ =工£訕小£山-也=1be!be/be/w=0

10、w=03戻7易見，當四種核昔酸符號b (/?/)在序列的上述第一、第二、第三個子序列上出現(xiàn)的頻n數(shù)風小心越接近相等吋，一處的譜值也就越接近于零。所以，基因外顯子序列的功率譜n曲線，在二頻率處具有較人的頻譜峰值(peak value),反映了在基因外顯子片段上，四種核 3莒酸符號在序列的三個了序列上分布的“非均衡性”。通常認為這種現(xiàn)象源于編碼基因序列 “密碼子”(code)使用的偏向性(bias)。雖然目前對此現(xiàn)象產生的“機理”還不是十分 地清楚，但是頻譜的3-周期性被普遍認為是可用于識別基因編碼序列(外顯子)的一個重要 的特征信息。3. 基因識別頻譜峰值特征的發(fā)現(xiàn)，或者頻譜與信噪比概念的引入，

11、其最終冃的是要探測、預報一個尚未被注釋的完整的dma序列的所侑基因編碼序列(外顯了)片段。圖4基于序列頻譜3周期性的的基因預測方法流程圖已經冇一些研究者提出了識別基因的算法（如參見及其后面的文獻）。冃前利川信噪 比的基因識別算法通常有兩種：一是固定長度窗ii滑動法叫另一是移動信噪比曲線識別 法。基于固定長度滑動窗口上頻譜曲線的基因識別方法:對一個dna序列s和它的指示序列他川 " w /皿二0丄2,n 1。収長度m （通 常取為3的倍數(shù),例如m=99, 129,255,513等）作為固定窗口長度。m -1 m -對任意n<o<n<n- ）,在以料為中心的長度為m的序

12、列片段w+2 2上（當n接近序列的兩端時，窗ii實際啟效長度可能會小于m）,作四個指示序列的離散m-fourier 變換(dft)i=n+_ .imubk=工 uhie j m , k=o,l,m-l2并求出它在理處總頻譜/9（/7；）,即2r m22m以3+叫+噲)把這樣得到的頻譜值"理），二0,1,2,2-1,經過標準化處理（即除以最大頻譜值m。酬"心；?。┎⒒爻銎漕l譜曲線圖5固定長度滑動窗口的頻譜p = 1心;)曲線(人類線粒體棊因，nc_012920_l.fasta)圖中紅色水平細線條是dna序列實際的基因外顯子的區(qū)間。滑動窗i i頻譜p(/2；)曲線的 峰與皋因外

13、顯了區(qū)間具冇“對應”關系?；赿na序列上“移動序列”信噪比曲線的基因識別方法：設已知dna序列s和它的指示序列uhn, bwl ,=1。對任意n(0<n<7v-l ),通常取3的倍數(shù)并逐漸增大。在的左邊一個長度為77的序列片段(), n-lj h,相應的了序列s稱為dna序列s的“移動子序列”，作該移動了序列對應的四 個指示序列的離散fourier變換(£>ft)i=n-_ j2jrik如幻=£唄弘昕，心0,1,1i=q并求出移動子序列n = 0丄，n 1上的信噪比rnnr r r 1r r1 r r r 1r r r 1p匕才+ “丁匕+如匕+&quo

14、t;c§n-i其中可旳為移動子序列sow_幻的功率譜的平均值£=o在坐標系中呦出移動序n列so*的信噪比曲線/?(稱為信噪比移動曲線(snr walk curve),見圖6)1412c)ccccns7000030003500400045005000550060006500nuclootide position n圖6 dna移動序列其指示序列的信噪比曲線。(人類線粒體基因，nc_012920_l.fasta)圖中紅色水平細線條是dna序列實際的基因外顯子的區(qū)間。dna序列的信噪比移動曲線 的峰、谷與基因外顯子區(qū)間的端點也具有較“明顯的”的對應關系。三、請研究的幾個問題：1.

15、功率譜與信噪比的快速算法對于很長的dna序列，在計算其功率譜或信噪比時，離散fourier變換（dft）的總體計 算量仍然很人，會影響到所設計的基因識別算法的效率。大家能否對映射，探求功率 譜與信噪比的某種快速計算方法？在基因識別研究中，為了通過引入更好的數(shù)值映射而獲取dna序列更多的信息，除了 上面介紹的wss映射外，實際上人們還研究過許多不同的數(shù)值映射方法。例如，著名的 z-curve映射（參見或者附件1）。試探討z-curve映射的頻譜為信噪比和voss映射下的頻 譜與信噪比之間的關系；此外，能否對實數(shù)映射，如：4to,ct1,gt2,tt3,也給出功率譜與信噪比 的快速計算公式？2對不

16、同物種類型基因的閾值確定對特定的基因類型的序列，將其信噪比的判別閾值取為r°=2,帶有一定的主 觀性、經驗性。對不同的基因類型，所選取的判別閾值也許應該是不同的。附件小給出了來 口于苦名的牛:物數(shù)據(jù)網(wǎng)站：hup://guidc/的兒個基因序列數(shù)據(jù)，另外也 給出了帶有編碼外顯子信息的10（）個人和鼠類的，以及200個哺乳動物類的基因序列的樣本 數(shù)據(jù)集合。大家還町以從生物數(shù)據(jù)庫下載更多的數(shù)據(jù)，找你們認為具冇代表性的棊因序列， 并對每類基因研究其閾值確定方法和閾值結果。此外，對按照頻譜或信噪比特征將編碼與非 編碼區(qū)間分類的有效性，以及分類識別時所產住

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

18、識別正確率的一 個重要原因。人們發(fā)現(xiàn)，對某些dn4序列而言，其部分編碼序列（外顯子），尤其是短的（長 度小于1 00/7/2 ）的編碼序列，就可能不具有頻譜或者信噪比顯著性。你們團隊能否總結，甚 至獨口提出一些識別基因編碼序列的其它特征指數(shù)，并對此做相關的分析？（2）“基因突變”是生物醫(yī)學等方面的一個關注熱點。基因突變包括dna序列中單個核 廿酸的替換，刪除或者插入等。那么，能否利用頻譜或信噪比方法去發(fā)現(xiàn)基因編碼序列可能 存在的突變呢？上而捉出的基于頻譜3-周期性的基因預測四個方而問題中，“快速算法”與“閾值確定” 是為設計基因預測算法做準備的。此外，在最后的延展性研究中，各隊也可以對你們口己認 為冇價值的其它相關問題展開探討。參考文獻：【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. bioinformatics 16, 1073-1081.3 kotlar, d., lavner, y., 20