Augustus指南(Trainning部分)

Augustus指南


官方 Tutorial Index

Augustus是一个真核生物基因预测软件,目前有网页服务端本地版,它基于Hidden-Markov Model(隐马尔科夫链模型HMM)(一个不错的HMM介绍博客)的预测方法,因此需要一个已经研究清楚的物种进行training(学习)之后再对新物种进行prediction(预测),用于trainning的物种应该和需要预测的物种具有较近的亲缘关系

特点:官方介绍

Input:

Trainning:

已知物种序列

核酸序列:fasta格式 ( genome.fa )
蛋白序列:fasta格式 ( proteins.aa )

Prediction

未知基因组序列

核酸序列:fasta格式 ( NewGenome.fa )

Output:

可能用到的其他软件

  1. Scipio
    下载地址中下载Scipio Version 1.4,压缩包中包括perl文件和说明文件

Scipio是利用蛋白序列标记出基因组结构的软件,提供核酸序列蛋白序列后,它可以为你解析出基因的结构。它基于BLAT的序列比对功能。scipio可以很好地跨contigs比对。
用途包括:

  1. 检测序列assembly的正确性
  2. 识别外显子区域,可变剪切
  3. 用已知物种蛋白序列描述未知物种的基因组结构
  1. Scipio依赖:perl中的YAML module
    Ubuntu安装yaml-perl命令:
sudo apt-get install libyaml-perl
  1. Scipio依赖:BLAT 下载地址,下载后的文件中有数个可执行文件,请添加进环境变量中,方便Scipio直接调用

Blat,全称 The BLAST-Like Alignment Tool,可以称为"类BLAST 比对工具",对于DNA序列,BLAT是用来设计寻找95%及以上相似至少40个碱基的序列。对于蛋白序列,BLAT是用来设计寻找80%及以上相似至少20个氨基酸的序列。速度快,共线性输出结果简单易读。Blat 把相关的呈共线性的比对结果连接成为更大的 比对结果,从中也可以很容易的找到 exons 和 introns。


Trainning Augustus

重要提示:如果是第一次training,按照下面所述的步骤可以正常完成,但是如果效果不佳,官方介绍中有一些小经验,本文没有包含这部分内容

  • Augustus基于Hidden-Markov Model(隐马尔科夫链模型)的预测方法,因此需要一个已经研究清楚的物种进行training(学习)之后再对新物种进行prediction(预测),用于trainning的物种应该和需要预测的物种具有较近的亲缘关系或者就是本物种。
  • Trainning的实质是用已知的知识调整模型中的种种参数,因此一般已知的知识必须包含各种可能的情况。具体地说,在预测基因的模型中,用于training的数据就必须拥有是基因的序列以及不是基因的序列(基因间隔区,内含子等)。
  • Trainning的数据会被分为两部分,一部分用于Trainning,另一部分用于检测调整参数后的模型。前者称为training set,后者称为test set

作为用于training Augustus的已知信息,一般可以有以下几种:

  1. 该物种已经存在的基因结构信息(比如:GenBank中的数据)
  2. 该物种ESTs信息加上相对于的基因组信息(可使用PASA
  3. 该物种de novo组装的RNA-seq数据加上对应的基因组数据
  4. 该物种已知的蛋白序列和对应的基因组数据
  5. 近源物种的基因结构信息
  6. 已经Trainning过的模型保存参数集后可以反复使用

1. 对于没有直接基因结构的情况,可以使用Scipio制作基因结构文件

如果已经有基因结构信息,且文件格式为Genbank格式的文件(*.gb)可以跳过此步。这种文件包含了:

  1. 核酸序列
  2. CDS编码区域的序列坐标信息
  3. UTR信息

1.1 Run Scipio

scipio.1.4.pl --blat_output=prot.vs.genome.psl genome.fa proteins.aa > scipio.yaml # takes ~7m

scipio.yaml包含了每个蛋白alignment的细节,可用于下一步提取

  • 输入的基因组5M,耗时约7min
  • --blat_output=prot.vs.genome.psl的含义就是记录时间。
  • genome.fa和proteins.aa为fasta格式。

1.2 Extract a GFF file,从scipio.yaml中提取GFF (General Feature Format) 基因结构文件

cat scipio.yaml | yaml2gff.1.4.pl > scipio.scipiogff
scipiogff2gff.pl --in=scipio.scipiogff --out=scipio.gff
cat scipio.yaml | yaml2log.1.4.pl > scipio.log     #产生log文件

scipio.gff文件是我们所希望得到的GFF文件,log文件可以帮助检测是否每个蛋白都得到标注。
scipio.gff长这样:

chr2R   Scipio  CDS     900562  900621  1.000   +       0       transcript_id "392"
chr2R   Scipio  CDS     904518  904880  1.000   +       0       transcript_id "392"
chr2R   Scipio  CDS     904940  905131  1.000   +       0       transcript_id "392"
chr2R   Scipio  CDS     905195  905263  1.000   +       0       transcript_id "392"
chr2R   Scipio  CDS     3595076 3596041 1.000   +       0       transcript_id "2517"
...

1.3 可选步骤 可视化基因结构(GBrowse)

GBrowse 下载并配置GBrowse,将GFF文件转化成GBrowse可读的文件。

1.4 将GFF文件转化成Genbank格式的文件

Augustus的etrainning学习软件需要输入Genbank格式的文件,这种文件包含了:

  1. 核酸序列
  2. CDS编码区域的序列坐标信息
  3. UTR信息

由于大部分的基因结构信息都是GFF或GTF格式,Augustus提供了gff2gbSmallDNA.pl进行格式转化

gff2gbSmallDNA.pl scipio.gff genome.fa 1000 genes.raw.gb

整合gff文件和fasta文件,genes.raw.gb文件中将会有每个基因序列和其1000bp的上下游基因间隔序列

1.5 去除不适合trainning的基因序列

genes.raw.gb文件已经符合Augustus的要求,但是其中可能包含:

  1. 可变剪切的基因
  2. 缺失启示密码子的基因
  3. 缺失终止密码子的基因
  4. 非正常的终止密码子(in-frame stop codons)

这些序列不适合用于软件的trainning,所以最好去除

etraining --species=generic --stopCodonExcludedFromCDS=true genes.raw.gb 2> train.err

得到问题序列文件train.err
--stopCodonExcludedFromCDS=true选项指将gb文件中的终止密码子视为cds中的一部分,这个选项需要视gb文件的来源而定,非scipio来源的文件可能需要设置为false
train.err长这样:

gene 186 transcr. 1 in sequence chr2R_549753-555284: Initial exon has length < 3!
gene 461 transcr. 1 in sequence chr2R_1034318-1036751: Initial Exon doesn't begin with start codon.
gene 567 transcr. 1 in sequence chr2R_1198437-1201521: Initial Exon doesn't begin with start codon.
gene 4537 transcr. 1 in sequence chr2R_1354359-1361857: Initial Exon doesn't begin with start codon.
gene 4783 transcr. 1 in sequence chr2R_1669860-1673241: Terminal exon doesn't end in stop codon. Variable stopCodonExcludedFromCDS set right?
gene 5161 transcr. 1 in sequence chr2R_2043765-2046183: Single exon gene doesn't begin with atg codon but with ccc
gene 6319 transcr. 1 in sequence chr2R_3734070-3735386: Initial Exon doesn't begin with start codon.
gene 3577 transcr. 1 in sequence chr2R_4767472-4770826: Initial Exon doesn't begin with start codon.

使用下面的命令来筛选掉这些问题序列:

cat train.err | perl -pe 's/.*in sequence (\S+): .*/$1/' > badgenes.lst
filterGenes.pl badgenes.lst genes.raw.gb > genes.gb
grep -c "LOCUS" genes.raw.gb genes.gb
# genes.raw.gb:594    原来序列有594个基因
# genes.gb:586           筛选后还剩586个

badgenes.lst为问题基因序列
genes.gb为没有问题的基因,依然是gb格式

2. 将用于Trainning的基因结构信息分为Trainning set和test set

随机的将genes.gb文件分为两个部分

randomSplit.pl genes.gb 100

该命令会生成genes.gb.test文件,它包含100条序列。
同时会生成genes.gb.trains文件,包含剩下的序列。

  • 为了满足test的统计学意义,test set必须足够大(100-200 genes),才能起到检测作用。
  • test set的选择必须满足随机,不能只是选择文件最前面100个或者最后面100个genes
  • randomSplit.pl是Augustus中包含的软件

3. 为你的物种产生meta参数文件

  • 模型中有两类参数:meta parametersparameters(一般参数)。
  • 前者包括拼接位点模型窗口的大小(the size of the window of the splice site models)马尔科夫模型顺序(the order of the Markov model)等;
  • 后者则是比如拼接位点模式的分布(the distribution of splice site patterns)编码区非编码区k-mer概率(the k-mer probabilities of coding and noncoding regions)等。
  • 一般来说meta parameters决定了parameters的算法,而且meta的数量比较少,而parameters则非常多
  • Training本质是调整模型中的parameters,但是不会改变meta parameters
new_species.pl --species=bug

假设我们的物种叫bug
这条命令会在环境变量AUGUSTUS_CONFIG_PATH指定的位置建立一系列文件和文件夹
就像这样:

creating directory /home/mario/augustus/trunk/config/species/bug/ ...
creating /home/mario/augustus/trunk/config/species/bug/bug_parameters.cfg ...
creating /home/mario/augustus/trunk/config/species/bug/bug_weightmatrix.txt ...
creating /home/mario/augustus/trunk/config/species/bug/bug_metapars.cfg ...
...

bug_parameters.cfg文件包含了meta parameters和parameters,还包括一些对output文件输出的控制选项。

编辑bug_parameters.cfg文件,将stopCodonExcludedFromCDS选项改为ture

4. 第一次Training

4.1 Trainning

这里将进行第一次Training,training所使用的meta parameters是默认值。

etraining --species=bug genes.gb.train

这条命令将会在环境变量AUGUSTUS_CONFIG_PATH指定的位置$AUGUSTUS_CONFIG_PATH/species/bug建立或者更新一些文件,这些文件是exon, intron和intergenic region的parameter文件(training的结果文件)

ls -ort $AUGUSTUS_CONFIG_PATH/species/bug/

展示出与模型相关的各个参数保存文件,
其中:bug_{intron,exon,igenic}.pbl是新生成的文件

4.2 Testing

那么,我们现在就可以使用之前划分出来的test set来检测Trainning的效果
Test set的文件为genes.gb.train

augustus --species=bug genes.gb.test | tee firsttest.out # takes ~1m
  • 这条命令对genes.gb.test中的所有核酸序列信息进行了预测,这里不会使用其中的结构信息的。
  • 之后又将预测信息和gb文件给出的结构信息进行了比较,从而评价trainning的模型。
  • 产生文件为firsttest.out*
grep -A 22 Evaluation firsttest.out

该命令展示了评价部分,输出这样的表格:

*******      Evaluation of gene prediction     *******

---------------------------------------------\
                 | sensitivity | specificity |
---------------------------------------------|
nucleotide level |       0.975 |        0.89 |
---------------------------------------------/

----------------------------------------------------------------------------------------------------------\
           |  #pred |  #anno |      |    FP = false pos. |    FN = false neg. |             |             |
           | total/ | total/ |   TP |--------------------|--------------------| sensitivity | specificity |
           | unique | unique |      | part | ovlp | wrng | part | ovlp | wrng |             |             |
----------------------------------------------------------------------------------------------------------|
           |        |        |      |                115 |                 76 |             |             |
exon level |    511 |    472 |  396 | ------------------ | ------------------ |       0.839 |       0.775 |
           |    511 |    472 |      |   40 |    5 |   70 |   43 |    3 |   30 |             |             |
----------------------------------------------------------------------------------------------------------/

----------------------------------------------------------------------------\
transcript | #pred | #anno |   TP |   FP |   FN | sensitivity | specificity |
----------------------------------------------------------------------------|
gene level |   118 |   100 |   53 |   65 |   47 |        0.53 |       0.449 |
----------------------------------------------------------------------------/

表格中有三个子表格:

  1. nucleotide level,sensitivity(预测到的百分率),specificity(其中正确的百分率)
  2. exon level, #pred total/unique(预测得到unique外显子总数),#anno total/unique(实际unique外显子总数),TP(正确的预测),FP(假阳性),FN(假阴性)
  3. gene level
  4. 100个基因中,预测到53个
  5. 83.9%的外显子被预测到
  6. 77.5%的外显子预测成功

5. optimize_augustus.pl的优化

  • 脚本optimize_augustus.pl可以通过不断迭代Tranning和Testing的过程,根据模型评价自动修改*_parameters.cfg文件中的meta parameters,使meta parameters的取值最优。
  • 这个脚本的运行会非常费时间,而且最终的效果一般只能提高准确度几个百分点,慎重使用
  • 在本例中Trainning基因组越5M,耗时越1d
    使用此脚本augustus和etraining命令必须在环境变量中

5.1 自动调整meta parameters

optimize_augustus.pl --species=bug genes.gb.train  # takes ~1d

本句命令最好加上nohup
命令完成后,meta parameters调整完毕。

5.2 retraining Augustus

meta parameters调整好后必须重新Trainning模型,否则没有任何意义。

etraining --species=bug genes.gb.train

同样你可以进行检验

augustus --species=bug genes.test.gb
  • 如果此时你的gene level的sensitivity还是低于20%说明Trainning set不够大,请添加数据。

  • 如果你获得了满意的Trainning结果,请开始prediction

posted @ 2015-05-12 15:35  南方丙丁  阅读(10676)  评论(0编辑  收藏  举报