准备阶段
训练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
对于转录本，我们通常会测一些转录组数据，有三种策略可以得到转录本。（这里暂时不考虑三代全长转录本）
- Trinity重头组装转录本
- 使用STAR + Trinity 获取转录本
- 使用STAR + StringTie + gffread 获取转录本
对于这三种策略，不推荐策略一，因为在有参考基因组的情况下，策略二不但计算效率高，而且能避免组装错误（多倍体等位基因之间相似度高）。对于策略二和策略三，我会推荐策略三。因为对于靠的比较近的基因，Trinity很可能会把这两个基因装成一个。
Fig1
并且利用该转录本作为输入训练SNAP模型，之后以SNAP模型作为输入，将转录组和同源蛋白作为证据而不是直接用作模型，我们再检查maker的结果, 也会发现使用StringTie进行组装的结果才是对的。
Fig2
因此使用策略二不但计算量大，而且有些情况下还会导致过近的转录本错误融合，反而影响了最终效果。尽管你可以通过调整参数，或者使用链特异性测序的方式来来优化结果，但STAR + StringTie 的组合同样也可以这样子做，因此最终我还是推荐策略三。
当然，这是我定性通过IGV浏览结果得出的结论，样本小，结论未必可靠，仅供参考。
训练阶段
假设我们准备的三个文件分别命名为, genome.fa, protein.fa 和不同组织来源的stringtie组装后的转录本或GFF文件。
注，如果是StringTie组装的GTF文件，需要做如下的转换(感谢上海欧易生物-鲍志贵提供的帮助）
gffread -E sample.gtf -o - | sed -e "s#transcript#match#g" -e "s#exon#match_part#g" > sample.gff
接着使用maker -CTL新建配置文件, 设置如下选项
genome=genome.fa
est=组织1.fa,组织2.fa,组织3.fa
est_gff=组织1.gff,组织2.gff,组织3.gff
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项，每一项的含义如下
- Length of the 5' UTR
- Fraction of splice sites confirmed by an EST alignment (
-c) - Fraction of exons that overlap an EST alignmetn(
-e) - Fraction of exons that overlap EST or Protein alignments(
-o) - Fraction of splice site confrimed by a ab-initio prediction(
-a) - Fraction of exons that overlap a ab-initio prediction(
-t) - Number of exons in the mRNA
- length of the 3' UTR
- 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模型，通过比较多个模型来得到最终结果。
关于SNAP软件细节可以参考使用MAKER进行基因注释(高级篇之SNAP模型训练）
推荐阅读:
网友评论