Augustus Training and Prediction

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

Augustus Training and Prediction》上有2条评论

  1. 陈老师:
    您好!
    看了您的博客,对您在NGS方面的深入研究以及相关分析工具的熟识程度深感佩服。我也是被生物信息学的有趣分析理念和方法所吸引,期望能在NGS以及其他生物信息方面做一些浅薄的工作。不得不说,您的博客在我生物信息入门方面起到了至关重要的作用。目前我正在对一个植物的基因组的gene models进行重新的注释,我看到您博客中有一些自己编写Perl语言脚本,比如split_hints_and_scaffolds_for_augustus.pl等,不知您能够提供这些脚本给我,非常感谢!另外,我在网络中看到您关于生物信息培训班的招生信息,若这些脚本的提供对您造成利益上以及其他权益的损害,您可以拒绝提供这些信息给我。
    最后,祝您工作顺利!

发表评论

电子邮件地址不会被公开。 必填项已用*标注