MAKER训练SNAP基因模型

准备阶段

训练SNAP模型,需要准备三个文件,分别是参考基因组序列,组装的转录本序列和同源蛋白序列。

对于参考基因组序列,我们要保证N50需要超过预期基因长度的中位数,否则注释效果不好。不过目前的基因组在三代测序的加持下,基本上都不是问题。

对于同源蛋白, 建议只用UniProt/Sprot的人工检查过的高质量蛋白序列,而不是盲目参考文献,使用近源物种的所有蛋白。除非你的近源物种都是模式物种,蛋白序列可靠性高,否则用错误的输入进行训练,数据越多反而错的越多。

我们可以在ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions/选择合适的uniprot_sprot数据, 然后将其输出为fasta格式。以植物为例

wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions/uniprot_sprot_plants.dat.gz
zcat uniprot_sprot_plants.dat.gz |\
    awk '{if (/^ /) {gsub(/ /, ""); print} else if (/^AC/) print ">" $2}' |\
    sed 's/;$//'> protein.fa

对于转录本,我们通常会测一些转录组数据,有三种策略可以得到转录本。(这里暂时不考虑三代全长转录本)

  1. Trinity重头组装转录本
  2. 使用STAR + Trinity 获取转录本
  3. 使用STAR + StringTie + gffread 获取转录本

对于这三种策略,不推荐策略一,因为在有参考基因组的情况下,策略二不但计算效率高,而且能避免组装错误(多倍体等位基因之间相似度高)。对于策略二和策略三,我会推荐策略三。因为对于靠的比较近的基因,Trinity很可能会把这两个基因装成一个。

Fig1

并且利用该转录本作为输入训练SNAP模型,之后以SNAP模型作为输入,将转录组和同源蛋白作为证据而不是直接用作模型,我们再检查maker的结果, 也会发现使用StringTie进行组装的结果才是对的。

Fig2

因此使用策略二不但计算量大,而且有些情况下还会导致过近的转录本错误融合,反而影响了最终效果,因此我最终推荐策略三。当然,这是我定性通过IGV浏览结果得出的结论,样本小,结论未必可靠,仅供参考。

训练阶段

假设我们准备的三个文件分别命名为, genome.fa, Trinity-GG.fasta 和 protein.fa

接着使用maker -CTL新建配置文件, 设置如下选项

genome=genome.fa
est=Trinity-GG.fasta
protein=protein.fa
est2genome=1
protein2genome=1

然后使用mpiexec -n 线程数 maker &> run.log运行程序。

处理结果后,我们新建一个snap目录训练模型

mkdir snap && cd snap
gff3_merge -d ../genome.maker.output/genome_master_datastore_index.log

使用makerzff构建输入文件

maker2zff -c 0.8 -e 0.8 -o 0.8 -x 0.2 genome.all.gff

maker2zff会根据AED(-x)和QI值进行过滤,其中QI值一共有9项,每一项的含义如下

  1. Length of the 5' UTR
  2. Fraction of splice sites confirmed by an EST alignment (-c)
  3. Fraction of exons that overlap an EST alignmetn(-e)
  4. Fraction of exons that overlap EST or Protein alignments(-o)
  5. Fraction of splice site confrimed by a ab-initio prediction(-a)
  6. Fraction of exons that overlap a ab-initio prediction(-t)
  7. Number of exons in the mRNA
  8. length of the 3' UTR
  9. Length of the protein sequence produced by the mRNA (-l)

如果QI值第二项为-1,表示没有支持该剪切位点的EST证据.

接着构建模型

fathom -categorize 1000 genome.ann genome.dna
fathom -export 1000 -plus uni.ann uni.dna
forge export.ann export.dna
hmm-assembler.pl snap . > ../snap.hmm

修改配置,然后重新运行。MAKER会自动处理冲突的部分,避免重复序列屏蔽等的一些重复计算。

genome=genome.fa
est=Trinity-GG.fasta
protein=protein.fa
snap=snap.hmm
est2genome=0
protein2genome=0

根据输出结果再一次训练模型

mkdir snap2 && cd snap2
gff3_merge -d ../genome.maker.output/genome_master_datastore_index.log
fathom -categorize 1000 genome.ann genome.dna
fathom -export 1000 -plus uni.ann uni.dna
forge export.ann export.dna
hmm-assembler.pl snap . > ../snap2.hmm

通常迭代2-3次就够了,毕竟我们可能还会训练AUGUSTUS和GeneMark模型,通过比较多个模型来得到最终结果。

参考资料: http://weatherby.genetics.utah.edu/MAKER/wiki/index.php/MAKER_Tutorial_for_WGS_Assembly_and_Annotation_Winter_School_2018#Training_ab_initio_Gene_Predictors

# 注释  MAKER 

评论

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×