Augustus操作指南

Augustus指南(Trainning部分)
目录

Augustus指南
官方 Tutorial Index
Input:
Trainning:
Prediction
Output：
可能用到的其他软件
Trainning Augustus

1. 对于没有直接基因结构的情况，可以使用Scipio制作基因结构文件
2. 将用于Trainning的基因结构信息分为Trainning set和test set
3. 为你的物种产生meta参数文件
4. 第一次Training
5. optimize_augustus.pl的优化
   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：
可能用到的其他软件
Scipio
下载地址中下载Scipio Version 1.4，压缩包中包括perl文件和说明文件

Scipio是利用蛋白序列标记出基因组结构的软件，提供核酸序列和蛋白序列后，它可以为你解析出基因的结构。它基于BLAT的序列比对功能。scipio可以很好地跨contigs比对。
用途包括：
检测序列assembly的正确性
识别外显子区域，可变剪切
用已知物种蛋白序列描述未知物种的基因组结构
Scipio依赖：perl中的YAML module
Ubuntu安装yaml-perl命令：
sudo apt-get install libyaml-perl
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的已知信息，一般可以有以下几种：

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

1. 对于没有直接基因结构的情况，可以使用Scipio制作基因结构文件
   如果已经有基因结构信息，且文件格式为Genbank格式的文件（*.gb）可以跳过此步。这种文件包含了：

核酸序列
CDS编码区域的序列坐标信息
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格式的文件，这种文件包含了：

核酸序列
CDS编码区域的序列坐标信息
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的要求，但是其中可能包含：

可变剪切的基因
缺失启示密码子的基因
缺失终止密码子的基因
非正常的终止密码子（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 parameters 和 parameters（一般参数）。
前者包括拼接位点模型窗口的大小(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 |
----------------------------------------------------------------------------/
表格中有三个子表格：

nucleotide level，sensitivity（预测到的百分率），specificity（其中正确的百分率）
exon level， #pred total/unique（预测得到unique外显子总数），#anno total/unique（实际unique外显子总数），TP（正确的预测），FP（假阳性），FN（假阴性）
gene level
100个基因中，预测到53个
83.9%的外显子被预测到
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吧