基因注释软件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.