基因注釋軟件GeneMarks和RAST

看過好多大神的博客，對自己的學習幫助很大，這是額的第一篇博客，其實是額的生物信息學作業，感覺還是有用的，分享給大家。

基因組注釋是在得到全基因組序列后首先要做的。它是利用生物信息學方法，對基因組所有基因的生物學功能進行功能注釋，包括基因預測和基因功能注釋兩個方面。目前已經有許多的基因預測工具或者在線注釋網站。基因預測的方法主要有3 種：(1)分析mRNA和EST數據直接得到結果；(2)通過相似性比對從已知基因和蛋白質序列得到間接證據；(3)基于各種統計模型和算法從頭預測，比如隱馬可夫模型。其中通過相似性比對得到預測基因的方法最常見。例如，現在流行的做法是先通過Glimmer、GeneMarks等軟件預測出基因組的ORF。然后通過Blast方法將ORF同其他物種的基因進行比對。有同源基因的ORF被注釋為同樣功能的基因，沒有同源性的ORF被舍去或注釋為假說蛋白（hypothetical protein）。由于注釋需要大量的數據庫，為了使注釋變得簡單，一些研究機構將不同功能的注釋軟件整合在一起，提供在線的注釋服務。如RAST，Xbase等，NCBI的PGAAP能提供人工的注釋服務。這些網站只需要用戶將序列和序列的所屬物種分類信息提交即可。注釋好的結果為gbk 格式文件（包含序列和注釋信息）

GeneMarks軟件的原理都是使用統計學模型的從頭預測(ab initio)方法，不依賴任何先驗知識和經驗參數，通過描述DNA序列中核苷酸的離散模型，利用編碼區和非編碼區的核苷酸分布概率不同來進行基因預測。GeneMarks是不需要人為干預和相關DNA或rRNA基因的資料即可對新的細菌基因組進行預測，測試表明GeneMarks對GeneBank數據庫中已注釋的枯草芽孢桿菌的預測準確度達到82.9%，而對已通過實驗方法證實注釋功能的大腸桿菌的預測高達93.8%，其對新測序基因組的預測與Glimmer存在同樣問題，即相當一部分基因在數據庫并不能發現同源，只能作為假蛋白基因存在。

如何在沒有明確實驗證據的前提下鑒定此類基因預測的準確性，切實可行的方法就是綜合利用多個預測軟件對預測結果進行比較，分析其中的異同點。

本研究主要以A.baumanniiACICU染色體序列為例對基因預測與注釋的方法進行分析，以找到合適的基因預測與注釋的方法。

2.   材料與方法（Methods and Materials）

下面利用從NCBI上下載的A.baumanniiACICU全基因組染色體序列（不包含質粒序列）（.fasta格式）為例，分別使用GeneMarks（http://topaz./GeneMark/genemarks.cgi）進行ORF（開放閱讀框）基因預測，RAST（http://rast./）進行功能基因（CDS）注釋，對比原結果進行分析。

2.1.使用GeneMarks進行ORF預測

（1）第一步是上傳A.baumaniiACICU染色體序列，并設置合適的參數，填加自己的郵箱。全部設置好之后，點擊[StartGeneMarks]開始注釋。如下圖所示：

（2）第一步上傳結束序列之后，會出現如下界面，提示序列已成功提交，注釋好的文件會發到所填郵箱。

2.2.使用RAST進行功能基因注釋

（1）上傳A.baumaniiACICU（.fasta格式）序列，上傳結束后點擊[Usethis data and go to step 2]進行下一步。如下圖所示：

（2）第二步填加必須的的參數，Domain選擇[Bacteria]，GeneticCode選擇[11]，然后點擊[Usethis data and go to step 3]進行下一步操作。如下圖所示：

（3）如下圖所示，選擇好合適的參數后點擊[Finishthe upload]，即可等待結果，注釋結束后，其會發郵件告知

3.   結果與討論（Results and Discussion）

3.1. 使用GeneMarks預測ORF的結果以及分析

使用GeneMarks進行預測后，生成了gms.out  gms.out.faa gms.out.fnn gms.out.ps四個文件：

其中gms.out文件如下顯示（其中一部分，使用linux系統cat或者head命令查看）：

  Gene      Strand    LeftEnd   RightEnd       Gene        Class

    #                                         Length

    1       -          76         468          393        1

    2       -         506        2974         2469        1

    3       -        3027        4109         1083        1

    4       -        4124       5272         1149        1

    5       -        5370        6767         1398        1

    6       +        7438        7572          135        1

    7       +        7602        7994          393        1

    8       +        8005        8325          321        1

    9       +        8331       10091         1761        1

   10       +       10182       11537         1356        1

 …………

 3711       +     3894879     3896006         1128        1

 3712       +     3896134     3896979          846        1

 3713       -     3897035     3897370          336        1

 3714       -     3897495     3898499         1005        1

 3715       -     3898842     3899849         1008        1

 3716       -     3900105    3901109         1005        1

 3717       +     3901366     3903297         1932        1

 3718       +     3903549     3904106          558        1

其中gms.out.faa氨基酸序列文件顯示如下（其中之一）：

>gene_3718|GeneMark.hmm|185_aa|+|3903549|3904106    >gi|184156320|ref|NC_010611.1|Acinetobacter baumannii ACICU, complete genome

MNFIDFITNFEQFLPILIQEYGAWVYAILFLIIFSETAFVFMFFLPGDSLLLTVGALCSV

VELMHLGYMITLLTVAATLGYIVNYSIGRHFGNRIFEAKSRFIKKEYLNKTNRYFLQHGG

KTILLARFIPFARSFAPLAAGSSNMSYGKFLIYNVAGAILWICILLTAGYLFGHALIQVT

DFVEN

其中gms.out.fnn核苷酸序列如下所示，起始密碼子為ATG，終止密碼子為TAATGA和TAG（其中之一）：

>gene_3718|GeneMark.hmm|558_nt|+|3903549|3904106    >gi|184156320|ref|NC_010611.1|Acinetobacter baumannii ACICU, complete genome

ATGAATTTTATTGATTTTATTACTAATTTTGAACAATTTTTACCTATTTTGATTCAGGAG

TATGGTGCATGGGTTTATGCCATACTCTTTTTGATTATTTTTTCTGAAACTGCTTTTGTG

TTTATGTTCTTTTTACCTGGAGATAGCTTACTTTTAACTGTAGGTGCACTGTGCTCGGTG

GTTGAACTGATGCATCTTGGTTATATGATTACTCTGCTCACCGTTGCAGCAACATTAGGC

TATATCGTCAATTATTCTATTGGCCGCCATTTTGGAAACCGTATTTTTGAAGCAAAATCA

CGTTTTATTAAAAAAGAATATTTGAATAAAACGAACCGCTATTTCTTGCAACATGGCGGTAAAACTATTCTTTTAGCACGTTTTATTCCTTTCGCACGTTCTTTTGCACCCCTCGCTGCCGGCTCAAGCAATATGAGCTATGGAAAATTTTTGATTTACAATGTGGCAGGAGCTATTTTGTGGATCTGCATCCTTTTAACGGCTGGCTACCTATTTGGCCATGCACTCATTCAAGTTACAGATTTTGTTGAAAATTAA

由此可知A.baumanniiACICU全基因組經GeneMarks預測到了3718個基因。

3.2.使用RAST進行功能基因注釋結果以及分析

     以上兩圖是使用RAST對A.baumannii ACICU染色體序列進行注釋的結果菌株A.baumanniiACICU染色體基因組經RAST功能基因注釋，共注釋到3683個功能基因。其中分布于不同功能子系統（457）的有1831個，確定的基因（non-hypothetical）有1736個，不確定（hypothrtical）的有95個；其余的編碼基因不分布于這些不同功能的子系統中，共有1852個，其中確定的有908個，不確定的有944個。

3.3.       綜合分析

對于A.baumaniiACICU染色體序列，由GeneMarks預測到3718個基因，由RAST注釋到3683個編碼蛋白基因，與原文獻結果含有預測基因數（ORF）為3758個，其中編碼蛋白質的基因數為3670個相比有所不同。其中預測基因數比原文獻少了有40個，差別較大，原文獻聯合使用GeneMarks與Glimmer對比預測，效果較好；注釋基因數相差比原文獻多13個，差別不大，原文獻中綜合使用COG與KEGG數據庫對預測到的蛋白序列進行注釋，說明RAST注釋結果還是比較可靠的。整個過程只是基因注釋的初始工作，要想得到完整準確的基因注釋結果，需要使用多個軟件進行注釋，對于不能準確注釋的基因還需要單獨進行注釋，最后綜合分析得到結果。

參考文獻：

1.  黃勇: 基于高通量測序的微生物基因組學研究. 中國人民解放軍軍事醫學科學院, 2013.

2.  AzizRK, Bartels D, Best AA, Dejongh M, Disz T, Edwards RA, Formsma K, Gerdes S,Glass EM, Kubal M: The RAST Server:Rapid Annotations using Subsystems Technology. Bmc Genomics 2008,9::75.

3.  夏偉: Gluconobacter oxydans 621H全基因組自動注釋結果的分析評估. 江南大學, 2013.

4.   BesemerJ, Lomsadze A, Borodovsky M: GeneMarkS:a self-training method for prediction of gene starts in microbial genomes.Implications for finding sequence motifs in regulatory regions. American Banker 2001,29:2607-2618.

5.    IaconoM, Villa L, Fortini D, Bordoni R, Imperi F, Bonnal RJP, Sicheritz-Ponten T, DeBellis G, Visca P, Cassone A, Carattoli A:Whole-genomepyrosequencing of an epidemic multidrug-resistant Acinetobacter baumanniistrain belonging to the European clone II group. Antimicrobial Agents and Chemotherapy 2008,52:2616-2625.