当前位置: 首页 > 工具软件 > Augustus > 使用案例 >

使用MAKER进行基因注释(高级篇之AUGUSTUS模型训练)

虞裕
2023-12-01

准备训练集和测试集

根据Augutus的官方教程,可靠的基因结构序列的要求如下:

  • 提供基因的编码部分,包含上游几KB。通常而言,基因越多,效果越好,至少准备200个基因以上。还得保证这些基因中要有足够多的外显子,这样子才能训练内含子。

  • 这些基因的基因结构一定要足够的准确。不过,也不需要百分百的正确,甚至注释都不需要特别的完整,只要保证起始密码子和终止密码子的准确是准确的即可。

  • 需要保证这些基因没有冗余,也就是说不同序列如果有几乎相同的注释后氨基酸序列,那么仅仅取其中一个(AUGUSTUS教程的建议是:保证任意两个基因在氨基酸水平上低于70%的相似度),这一步既可以避免过度拟合现象,也能用于检验预测的准确性

  • 一条序列允许有多个基因,基因可以在正链也可以在负链,但是这些基因间不能有重叠,每个基因只要其中一个转录本,存放格式是GenBank

之后随机将注释数据集分成训练集和测试集,为了保证测试集有统计学意义,因此测试集要足够多的基因(100~200个),并且要足够的随机。

基因结构集的可能来源有:

  • Genbank
  • EST/mRNA-seq的可变剪切联配, 如PASA
  • 临近物种蛋白的可变剪切联配,如GeneWise
  • 相关物种的数据
  • 预测基因的迭代训练

由于目前的转录组数据比较多,我先用Trinity对转录组数据进行从头组装,然后用PASA将Trinity组装的转录本回贴到参考基因组上

Launch_PASA_pipeline.pl \
      -c alignAssembly.config \
      -C \ # 创建数据库
      -R \ # 运行alignment/assembly 流程
      -g reference.fasta \
      -t Trinity.fasta \
      --ALIGNERS blat,gmap \ # 可以只用一个
      --CPU 12 &> &> pasa_$(date +%Y-%m-%d-%H-%M).log &

alignAssembly.config可以复制PASApipeline-v2.3.3/pasa_conf的pasa.alignAssembly.Template.txt, 修改如下内容

DATABASE=<__DATABASE__> # MySQL中的数据库名, 也是最后文件名的前缀

# 如下调整的是PASApipeline-v2.3.3/scripts/validate_alignments_in_db.dbi参数
script validate_alignments_in_db.dbi
validate_alignments_in_db.dbi:--MIN_PERCENT_ALIGNED=<__MIN_PERCENT_ALIGNED__>
validate_alignments_in_db.dbi:--MIN_AVG_PER_ID=<__MIN_AVG_PER_ID__>

这一步得到的<prefix>.assemblies.fasta<prefix>.pasa_assemblies.gff3

从PASA组装中提取ORF(开放阅读框)

PASApipeline-v2.3.3/scripts/pasa_asmbls_to_training_set.dbi \
    --pasa_transcripts_fasta <prefix>.assemblies.fasta.transdecoder.genome.gff3\
    --pasa_transcripts_gff3 <prefix>.pasa_assemblies.gff3

最后会生成一系列文件,如下

  • <prefix>.assemblies.fasta.transdecoder.cds/pep/gff3/bed: 虽然不再基因组上,但是根据转录本信息,有可能是编码区的结果
  • <prefix>.assemblies.fasta.transdecoder.genome.bed/gff3: 对应基因组序列的基因模型

我们需要的是后者。PASA组装的GFF3文件的属性一栏分为:

  • complete:
  • 5primer_partial: 缺少起始密码子,翻译到5'端为止
  • 3primer_partial: 缺少终止密码子,翻译到3'端为止
  • internal: 缺少起始密码子和终止密码子

格式转换

我们基于PASA的结果,将GFF3格式转换成Genbank格式

augustus/scripts/gff2gbSmallDNA.pl \
    <prefix>.assemblies.fasta.transdecoder.genome.bed/gff3 \ #gff3 from pasa_asmbls_to_training_set.dbi
    reference.fa \ # reference genome sequence
    1000 \ # flank
    <prefix>.gene.raw.gb # output

这一步的gff2gbSmallDNA会忽略掉UTR序列。

过滤可能错误的基因结构

转录组数据不像以前的EST序列,最后组装的基因数不多。如果同时使用多个组织的转录组,就会组装出一万条以上的基因,所以我们更加关注于质量。

先创建初始化的物种HMM文件

augustus/scripts/new_species.pl \
    --species=for_bad_genes_removing \
    --AUGUSTUS_CONFIG_PATH=/path/to/augustus/config

然后尝试训练,捕捉错误

augustus/bin/etraining \
    --species=for_bad_genes_removing \
    --stopCodonExcludedFromCDS=false \
    <prefix>.gene.raw.gb 2> train.err

提取错误,并进行过滤

egrep -o 'ctg[[:digit:]]+_[[:digit:]]+-[[:digit:]]+' train.err > badgenes_list.txt
augustus/scripts/filterGenes.pl badgenes_list.txt <prefix>.gene.raw.gb > <prefix>.genes.gb

初步训练

将上一步过滤后的文件随机分成两份,测试集和训练集。其中训练集的数目根据gb的LOCUS数目决定,至少要有200.

augustus/scripts/randomSplit.pl <prefix>.gene.gb 200 # 200为测试集的基因数

创建初始化HMM参数文件,并进行训练

augustus/scripts/new_species.pl --species aha_train1
augustus/bin/etraining --species=<prefix> <prefix>.genes.gb.train | tee train.out

用测试数据集检验预测效果。这里可以比较我们训练的结果,和近缘已训练物种的训练效果

augustus --species=<prefix> <prefix>.genes.gb.test | tee train_test1.out
augustus --species=arabidopsis <prefix>.genes.gb.test | tee ara_test.out

可以用tail -n 40 xx_test.out查看预测结果的统计。要看着最后的结果,需要解释几个统计学概念:

  • TP(True Positive): 预测为真,事实为真
  • FP(False Positive): 预测为真,事实为假
  • FN(False Negative): 预测为假,事实为真
  • TN(True Negative): 预测为假,事实为假

基于上述,引出下面两个概念。"sensitivity"等于TP/(TP+FP), 是预测为真且实际为真的占你所有认为是真的比例."specificity"等于TN/(TN+FN), 是预测为假且实际为假的占你所有认为是假的比例。我们希望在预测中,尽可能地不要发生误判,也就是没有基因的地方不要找出基因,有基因的地方不要漏掉基因。

优化训练参数

很有可能的一种情况是,我们第一次的训练结果没有已有训练的效果好,所以我们需要进行循环训练找到最优参数

augustus/scripts/optimize_augustus.pl --species=<prefix> \
    --cpus=12 \ #受限于--kfold=k
    --rounds=12 \
    <prefix>.genes.gb.train

评估优化训练结果

使用optimize_augustus.pl训练的参数未必会比前一步的训练结果好,因此需要

augustus/bin/etraining --species=<prefix> <prefix>.genes.gb.train
augustus/bin/augustus --species=<prefix> <prefix>.genes.gb.test | tee train_test2.out

比较前后的"sensitivity"和"specificity"参数,选择其中表现比较好的一个,作为最终用于预测的HMM文件。

 类似资料: