1. Augustus training
首选，需要有至少 200 个完整基因模型的数据。 例如： 使用 genome-guided 方法进行 trinity 有参考基因组的 de novo 组装；再使用 PASA 将组装出来的 inchworm 序列比对到基因组； 再提取完整基因模型数据，得到文件 trainingSet_CompleteBest.gff3 。
然后，将 GFF3 文件转换为 GeneBank 文件：
$ /opt/biosoft/augustus-3.0.3/scripts/gff2gbSmallDNA.pl trainingSet_CompleteBest.gff3 ../../genome.fasta 50 trainingSetComplete.gb 使用 genebank 格式的文件预先进行一次 Augustus training，一般会得到一些错误提示 $ etraining --species=generic --stopCodonExcludedFromCDS=false trainingSetComplete.gb 2> train.err 根据错误提示，提取出有错误的基因模型 $ cat train.err | perl -ne 'print "$1\n" if /in sequence (\S+):/' > badlist 去除错误的基因模型，得到能用于 training 的基因模型 $ /opt/biosoft/augustus-3.0.3/scripts/filterGenes.pl badlist trainingSetComplete.gb > training.gb
在进行 Augustus training 之前，最好保证这些基因模型两两之间在氨基酸水平的 identity < 70% ：
先获取这些基因模型的 Proteins 序列 $ /opt/biosoft/augustus-3.0.3/scripts/gbSmallDNA2gff.pl training.gb > training.gff2 $ perl -ne ‘print “$1\n” if /gene_id \”(\S+?)\”/’ training.gff2 | uniq > trainSet.lst $ perl extract_genes.pl trainingSet_CompleteBest.gff3 trainSet.lst > training.gff3 $ /opt/biosoft/EVM_r2012-06-25/EvmUtils/gff3_file_to_proteins.pl training.gff3 genome.fasta prot > training.protein.fasta $ perl -p -i -e ‘s/(>\S+).*/$1/’ training.protein.fasta $ perl -p -i -e ‘s/\*//’ training.protein.fasta 对这些 proteins 序列构建 blast 数据库，并将 proteins 序列比对到此数据库 $ makeblastdb -in training.protein.fasta -dbtype prot -title training.protein.fasta -parse_seqids -out training.protein.fasta $ blastp -db training.protein.fasta -query training.protein.fasta -out training.protein.fasta.out -evalue 1e-5 -outfmt 6 -num_threads 8 提取 blast 结果中 identity >= 70% 的比对信息，identity 高的 proteins 序列仅保留一条 $ grep -v -P “\t100.00\t” training.protein.fasta.out | perl -ne ‘split; print if $_[2] >= 70’ > blast.out $ perl delete_high_identity_proteins_in_training_gff3.pl training.protein.fasta blast.out training.gff3 > training.new.gff3 获得去除了冗余的基因模型 $ /opt/biosoft/augustus-3.0.3/scripts/gff2gbSmallDNA.pl training.new.gff3 genome.fasta 50 training.gb
进行 Augustus Training without CRF
初始化本物种的 HMM 文件 $ /opt/biosoft/augustus-3.0.3/scripts/new_species.pl --species=my_species 如果有 RNA-Seq 数据，获得能用于 training 的基因模型一般会有好几千个。我们将基因模型随机分成两份： 第一份 300 个基因，用于检测 training 的精确性； 另外一份含有更多的基因，用于进行 Augustus training。 $ /opt/biosoft/augustus-3.0.3/scripts/randomSplit.pl training.gb 300 使用大份的 traing.gb.train 进行 Augustus Training $ etraining --species=my_species training.gb.train >train.out 根据输出结果 train.out 来修正参数文件 species_parameters.cfg 中终止密码子的频率 $ tag=$(perl -ne 'print "$1" if /tag:\s+\d+\s+\((\S+)\)/' train.out) $ perl -p -i -e "s#/Constant/amberprob.*#/Constant/amberprob $tag#" /opt/biosoft/augustus-3.0.3/config/species/lentinula_edodes/lentinula_edodes_parameters.cfg $ taa=$(perl -ne 'print "$1" if /taa:\s+\d+\s+\((\S+)\)/' train.out) $ perl -p -i -e "s#/Constant/ochreprob.*#/Constant/ochreprob $taa#" /opt/biosoft/augustus-3.0.3/config/species/lentinula_edodes/lentinula_edodes_parameters.cfg $ tga=$(perl -ne 'print "$1" if /tga:\s+\d+\s+\((\S+)\)/' train.out) $ perl -p -i -e "s#/Constant/opalprob.*#/Constant/opalprob $tga#" /opt/biosoft/augustus-3.0.3/config/species/lentinula_edodes/lentinula_edodes_parameters.cfg 根据 training 的结果，进行第一遍精确性评估 $ augustus --species=my_species training.gb.test > test.1.out 再将 training.gb.train 分成两份： 第一份 800 个基因，剩下基因为另外一份。这两份基因都用于 Optimizing meta parameters of AUGUSTUS $ /opt/biosoft/augustus-3.0.3/scripts/randomSplit.pl training.gb.train 800 使用上面的 2 份基因模型，进行 Augustus training 的优化 $ /opt/biosoft/augustus-3.0.3/scripts/optimize_augustus.pl --species=my_species --rounds=5 --cpus=16 --kfold=16 training.gb.train.test --onlytrain=training.gb.onlytrain --metapars=/opt/biosoft/augustus-3.0.3/config/species/my_species/my_species_metapars.cfg > optimize.out 按如上参数，则程序会对 my_species_metapars.cfg 文件中的 28 个参数进行优化，总共优化 5 轮或有一轮找不到可优化的参数为止。每进行一个参数的优化时： 将 training.gb.train.test 文件中 800 个基因随机分成 16 等份，取其中 15 等份和 training.gb.onlytrain 中的基因模型一起进行 training，剩下的 1 等份用于精确行评估； 要对每个等份都进行一次精确性评估；使用 16 个 CPU 对 16 个等份并行进行 training 和 精确性评估，得到平均的精确值；优化的每个参数会分别 3~4 个值，每个值得到一个 training 的精确值，对参数的多个设定值进行比较，找到最佳的值。 此外， training 的精确值的算法： accuracy value = (3*nucleotide_sensitivity + 2*nucleotide_specificity + 4*exon_sensitivity + 3*exon_specificity + 2*gene_sensitivity + 1*gene_specificity) / 15 。 优化参数完毕后，需要再此使用 training.gb.train 进行一遍 training $ etraining --species=my_species training.gb.train 再进行第二遍的精确性评估，一般进行优化后，精确性会有较大幅度的提高 $ augustus --species=my_species training.gb.test > test.2.withoutCRF.aout
进行 Augustus Training with CRF(Conditional Random Field)
在进行 Training with CRF 之前，先备份非 CRF 的参数文件 $ cd /opt/biosoft/augustus-3.0.3/config/species/species/ $ cp species_exon_probs.pbl lentinula_edodes_exon_probs.pbl.withoutCRF $ cp species_igenic_probs.pbl lentinula_edodes_igenic_probs.pbl.withoutCRF $ cp species_intron_probs.pbl lentinula_edodes_intron_probs.pbl.withoutCRF $ cd - 在 training 的时候，加入 --CRF 参数，这样进行 training 比较耗时 $ etraining --species=my_species training.gb.train --CRF=1 1>train.CRF.out 2>train.CRF.err 再次进行精确行评估 $ augustus --species=my_species training.gb.test > test.2.CRF.out 比较 CRF 和 非 CRF 两种情况下的精确性。一般情况下，CRF training 的精确性要高些。若 CRF training 的精确行低些，则将备份的参数文件还原回去即可。
2. Training UTR parameters for Augustus
Training UTR 有利于利用 exonpart hint。 当使用 exonpart hint 时，对基因结构预测有显著提高。
首先，需要获得用于 training 的基因模型，这些基因模型需要同时带有 5’UTR 和 3’UTR 。
$ mkdir utr $ cd utr 从用于 Augustus Training 的数据中提取同时带有 5'UTR 和 3'UTR 的基因模型。 $ perl -ne 'print "$2\t$1\n" if /.*\t(\S+UTR)\t.*ID=(\S+).utr/' ../training.new.gff3 | sort -u | perl -ne 'split; print "$_[0]\n" if ($g eq $_[0]); $g = $_[0];' > bothutr.lst $ perl -e 'open IN, "bothutr.lst"; while () {chomp; $keep{$_}=1} $/="//\n"; while (<>) { if (/gene=\"(\S+?)\"/ && exists $keep{$1}) {print} }' ../training.gb.test > training.gb.test $ perl -e 'open IN, "bothutr.lst"; while ( ) {chomp; $keep{$_}=1} $/="//\n"; while (<>) { if (/gene=\"(\S+?)\"/ && exists $keep{$1}) {print} }' ../training.gb.train > training.gb.train 使用 traing.gb.train 中的基因模型进行 Training UTR parameters $ /opt/biosoft/augustus-3.0.3/scripts/randomSplit.pl training.gb.train 400 $ mv training.gb.train.train training.gb.onlytrain $ optimize_augustus.pl --species=my_species --cpus=16 --rounds=5 --kfold=16 --UTR=on --trainOnlyUtr=1 --onlytrain=training.gb.onlytrain --metapars=/opt/biosoft/augustus-3.0.3/config/species/lentinula_edodes/lentinula_edodes_metapars.utr.cfg training.gb.train.test > optimize.out 进行 without CRF 和 with CRF 的精确性比较，选取其中精确性较高的参数文件 $ etraining --species=species --UTR=on training.gb.train $ augustus --species=species --UTR=on training.gb.test > test.withoutCRF.out $ etraining --species=species --UTR=on training.gb.train --CRF=1 $ augustus --species=species --UTR=on training.gb.test > test.CRF.out
3. Creating hints from RNA-Seq data with Tophat
先使用 tophat 将 RNA-Seq 数据比对到屏蔽了重复序列的基因组，得到 bam 文件。对 bam 文件进行转换，得到 intron 和 exonpart hints。
得到 intron 的 hints 信息 $ /opt/biosoft/augustus-3.0.3/bin/bam2hints --in=tophat.bam --out=hints.intron.gff --maxgenelen=30000 --intronsonly 得到 exonpart 的 hints 信息 $ /opt/biosoft/augustus-3.0.3/bin/bam2wig tophat.bam tophat.wig $ cat tophat.wig | /opt/biosoft/augustus-3.0.3/scripts/wig2hints.pl --width=10 --margin=10 --minthresh=2 --minscore=4 --prune=0.1 --src=W --type=ep --UCSC=unstranded.track --radius=4.5 --pri=4 --strand="." > hints.exonpart.gff cat hints.intron.gff hints.exonpart.gff > hints.rnaseq.gff
4. Creating hints from Protiens with Exonerate
要先得到临近物种的 pritein 序列，然后使用 tblastn 将这些 protein 序列比对到基因组，再将相似性较高的序列使用 exonerate 比对到基因组，对 exonerate 结果进行分析，得到 hint 信息。
将 protein 序列比对到屏蔽了重复序列和转录组序列匹配区域的基因组。这样得到的 hints 主要位于转录组中未表达的基因处，以免和转录组的 hint 有冲突，或得到连接 2 个相邻基因的错误 hint。 $ makeblastdb -in genome.estMasked.fa -dbtype nucl -title genome -parse_seqids -out genome $ blast.pl tblastn genome proteins.fasta 1e-8 24 blast 5 对 blast 的结果进行分析，挑选相似性高的 pritein 序列用于 exonerate 分析 $ perl tblastn_xml_2_exonerate_commands.pl --exonerate_percent 50 --exonerate_maxintron 20000 --cpu 24 blast.xml proteins.fasta genome.estMasked.fa 将 exonerate 结果转换为 hint 文件 $ /opt/biosoft/augustus-3.0.3/scripts/exonerate2hints.pl --in=exonerate.out --maxintronlen=20000 --source=P --minintron=31 --out=hints.protein.gff
5. Creating hints from RepeatMasker output
在转座子重复区，不倾向于存在编码区，该区域可以作为 nonexonpart hint。将 RepeatMakser 结果转换为此类 hint，有力于 exon 的预测。这样比使用屏蔽了重复序列的基因组进行 Augustus gene prediction 要好。
不对 Simple_repeat, Low_complexity 和 Unknown 这 3 类进行屏蔽，因为这些区域存在 CDS 的可能性相对较大。 $ grep -v -P "position|begin|Unknown|Simple_repeat|Low_complexity" genome.repeat.out | perl -pe 's/^\s*$//' | perl -ne 'chomp; s/^\s+//; @t = split(/\s+/); print $t[4]."\t"."repmask\tnonexonpart\t".$t[5]."\t".$t[6]."\t0\t.\t.\tsrc=RM\n";' > hints.repeats.gff
6. Run Augustus predictions parallel
推荐使用 hints 进行 Augustus gene prediction，这样在有 hints 区域的基因预测会非常准确。 hints 越准确，覆盖基因组的区域越广泛，则基因预测越准确。
hint 的种类很多： 有人工进行确定的 hints； 使用 RNA-Seq 数据得到 intron 和 exonpart 类型的 hint； 使用 protein 得到 intron 和 cdspart 类型的 hint； 使用 RepeatMasker 得到 rm 类型的 hints。
对这些不同类型的 hint 要采用对应的参数来指导基因预测。 这些参数主要用于指导 augustus 对相应类型的 hints 进行得分奖罚。相关的参数文件位于 /opt/biosoft/augustus-3.0.3/config/extrinsic/ 文件夹下。
一个推荐的配置参数文件如下：
[SOURCES]
M RM P E W
[SOURCE-PARAMETERS]
[GENERAL]
start 1 1 M 1 1e+100 RM 1 1 P 1 1 E 1 1 W 1 1
stop 1 1 M 1 1e+100 RM 1 1 P 1 1 E 1 1 W 1 1
tss 1 1 M 1 1e+100 RM 1 1 P 1 1 E 1 1 W 1 1
tts 1 1 M 1 1e+100 RM 1 1 P 1 1 E 1 1 W 1 1
ass 1 1 M 1 1e+100 RM 1 1 P 1 1 E 1 1 W 1 1
dss 1 1 M 1 1e+100 RM 1 1 P 1 1 E 1 1 W 1 1
exonpart 1 .992 M 1 1e+100 RM 1 1 P 1 1 E 1 1 W 1 1.005
exon 1 1 M 1 1e+100 RM 1 1 P 1 1 E 1 1 W 1 1
intronpart 1 1 M 1 1e+100 RM 1 1 P 1 1 E 1 1 W 1 1
intron 1 .34 M 1 1e+100 RM 1 1 P 1 1000 E 1 1e5 W 1 1
CDSpart 1 1 0.985 M 1 1e+100 RM 1 1 P 1 1e5 E 1 1 W 1 1
CDS 1 1 M 1 1e+100 RM 1 1 P 1 1 E 1 1 W 1 1
UTRpart 1 1 0.985 M 1 1e+100 RM 1 1 P 1 1 E 1 1 W 1 1
UTR 1 1 M 1 1e+100 RM 1 1 P 1 1 E 1 1 W 1 1
irpart 1 1 M 1 1e+100 RM 1 1 P 1 1 E 1 1 W 1 1
nonexonpart 1 1 M 1 1e+100 RM 1 1.01 P 1 1 E 1 1 W 1 1
genicpart 1 1 M 1 1e+100 RM 1 1 P 1 1 E 1 1 W 1 1
对上面配置文件的内容简单讲解：
看此配置文件之前，要看 hint 文件内容中第 9 列有 src=P 这样的信息。 Augustus 通过此信息来确定 hint 的来源，然后根据 extrinsic 配置文件中相应的参数来处理 hints 文件中的内容。
[SOURCES] 设置 hint 的来源
M RM P E W 不同来源的 hint 类型使用空格分隔
M ： 手工锚定的 hint
RM ： RepeatMasker 获得的 hint
P ： 来源于 protein 的 hint
E ： 来源于 EST 的 hint
W ： 使用 RNA-Seq 进行 wiggle track coverage 分析得到的 hint
[SOURCE-PARAMETERS] 设置 hins 来源的参数
E 1group1gene 第 1 列是 hints 的来源， 第 2 列表明对该来源的 hints 进行的处理。有 2 种处理方式：
individual_liability
例如， 1 个 group 在 hint 文件中包含多行，即由多个 hints 组成（比如， 1 个基因有多个 intron），当其中一个 hint 不正确时，默认情况下则会对整个 group 进行忽略。而加入该参数则仅忽略错误的 hints。
1group1gene:
对 1 个 group，试图预测 1 个基因。适合使用比较完整的 transcripts 序列做的 hint。
[GENERAL] 不同类型不同来源 hints 的参数设置
前 3 列是对不同类型的 hints 的奖罚系数的设置：
第 1 列
hint 的类型。 当 hints 文件中没有该类型的 hint 时， 则后面不同来源的 hints 数值都使用 1 。
第 2 列
奖励系数。 当该类型 hint 全部匹配的时候，则进行奖励。奖励系数由该值和相应 hints 来源的设置一起决定。
第 3 列
惩罚系数。 每当有 1 个碱基不匹配，则得分乘以该系数。若有 100 个碱基不匹配，是该列值的 100 次方，因此，该值设置一般略比 1 小。该值越小，越能增加基因预测的 specficity。
后面的列，表示不同的来源的 hints 的奖励惩罚系数，每个来源的 hints 设置分 3 列：
第 1 列
设置 hint 的来源，与 hint 文件中 src 的值对应
第 2 列
对 hint 的得分(对应 hint 文件第 6 列)进行了分级，该值表明分成了多少级。 若该值为 1， 则表示不分级，那么当 hint 所有的碱基都匹配的时候，则进行奖励，奖励系数乘以第 3 列的值；若该值大于 1，则将 hint 文件第 6 列的分数分成了多级；若 hint 文件第 6 列没有得分，则该出需要设置为 1 。
第 3 列
若第 2 列值为 1， 则该列只有一个值，奖励系数乘以该值； 若第 2 列不为 1， 则此列分成 2 列。例如：
D 8 1.5 2.5 3.5 4.5 5.5 6.5 7.5 0.58 0.4 0.2 2.9 0.87 0.44 0.31 7.3
第一列为 D， 表明 DIALIGN 类型的 hint；
第二列为 8， 表明根据其 hint 文件的第 6 列，将奖励分成了 8 个级别；
由于第二列不是 1 ， 第三列分成了 2 列。 前一列是 7 个数值，将 hint 的打分分成了 8 个级别；后一列是这 8 个级别的奖励系数乘积。
在存在 hints.gff, extrinsic.cfg 和 genome.fasta 这 3 个文件，以及 HMM 文件 training 完毕后，即可开始进行 Augustus Training。当基因组比较大时，最好进行并行计算：
将 hints.gff 文件内容和 genome.fasta 内容进行分割。不对完整的序列进行切断。此程序将基因组序列按长度进行排序后，将序列写入到一个个分割的fasta文件中。每当写入到一个fasta文件的序列长度大于设置的值时，则将下一条序列下如下一个fasta文件。同时，也将相应的 hints 信息写入对应的 hint 文件。
$ perl split_hints_and_scaffolds_for_augustus.pl --minsize 1000000 --output split genome.fasta hints.gff
并行计算
$ for x in `ls split/*.fa | perl -pe 's/.*\///; s/.fa//' | sort -k 1.14n`
do
echo "augustus --species=my_species --extrinsicCfgFile=extrinsic.cfg --alternatives-from-evidence=true --hintsfile=split/$x.hints --allow_hinted_splicesites=atac --alternatives-from-evidence=true --gff3=on --UTR=on split/$x.fa > split/$x.out" >> command_augustus.list
done
$ ParaFly -c command_augustus.list -CPU 12
合并结果
$ for x in `ls split/*.out | sort -k 1.20n`
do
cat $x >> aug.out
done
$ join_aug_pred.pl aug.out > aug.gff3