使用MAKER进行注释: 如何避免多轮MAKER时的重复运算

通常而言,我们会运行不只一轮的MAKER。如果参考组序列没有变化,那么有一些计算只需要做一次就行了,例如将EST, Repeat和Protein序列比对到参考基因组,得到它们对应的位置。

我们有三种方法可以避免不必要的运算,第一种方法是直接修改配置文件,让MAKER重复利用之前的运行结果;第二种方式是利用之前输出的GFF文件,通过配置"Re-annotation Using MAKER Derived GFF3"里的选项来跳过对应的计算;第三种方法于是利用之前输出的GFF文件,从中提取EST/Repeat/Protein的位置信息保存为GFF文件,通过配置"est_gff", "protein_gff", "rm_gff"来避免重新计算位置信息。

后续分析建立MAKER高级篇-SNAP模型训练基础上,也就是通过protein和est序列直接输出基因模型,然后训练出初步的HMM模型

方法1

方法1最为简单,我们只需要修改之前的maker_opts.ctl里的参数,然后重新运行即可。运行时会输出如下的警告信息。

警告信息

注意: MAKER是通过对比maker_opt.ctl里的配置信息和自己运行时记录的maker_opts.log来判断哪些参数发生了改变。因此,如果SNAP第二次训练生成的文件,要是和上一次命名相同,那么它会认为你这次输入的模型文件和上次相同,就会跳过SNAP预测这一步。

实际运行时,MAKER会跳过BLAST步骤,但是依旧会调用"exonerate"来处理BLAST结果。

方法2

如果你不小心把maker的输出文件删掉了,但是你保留着之前gff3_merge默认参数输出的文件,那么你可以使用该文件来跳过BLAST和Exonerate运算。

Step1: 配置"Re-annotation Using MAKER Derived GFF3"里的参数

#-----Re-annotation Using MAKER Derived GFF3
maker_gff=round1.gff #MAKER derived GFF3 file
est_pass=1 #use ESTs in maker_gff: 1 = yes, 0 = no
altest_pass=1 #use alternate organism ESTs in maker_gff: 1 = yes, 0 = no
protein_pass=1 #use protein alignments in maker_gff: 1 = yes, 0 = no
rm_pass=1 #use repeats in maker_gff: 1 = yes, 0 = no
model_pass=1 #use gene models in maker_gff: 1 = yes, 0 = no
pred_pass=1 #use ab-initio predictions in maker_gff: 1 = yes, 0 = no
other_pass=1 #passthrough anyything else in maker_gff: 1 = yes, 0 = no

此处的round1.gff通过gff3_merge从上一论的maker输出中提取,代码如下

gff3_merge -d genome.maker.output/genome_master_datastore_index.log -o round1.gff

Step2: 将"EST Evidence"和"Protein Homology Evidence"里的配置清空,如下

#-----EST Evidence (for best results provide a file for at least one)
est= #set of ESTs or assembled mRNA-seq in fasta format
altest= #EST/cDNA sequence file in fasta format from an alternate organismest_gff= #aligned ESTs or mRNA-seq from an external GFF3 file
altest_gff= #aligned ESTs from a closly relate species in GFF3 format

#-----Protein Homology Evidence (for best results provide a file for at least one)
protein=  #protein sequence file in fasta format (i.e. from mutiple organisms)
protein_gff=  #aligned protein homology evidence from an external GFF3 file

#-----Repeat Masking (leave values blank to skip repeat masking)
model_org= #select a model organism for RepBase masking in RepeatMasker
rmlib= #provide an organism specific repeat library in fasta format for RepeatMasker
repeat_protein= #provide a fasta file of transposable element proteins for RepeatRunner
rm_gff= #pre-identified repeat elements from an external GFF3 file
prok_rm=0 #forces MAKER to repeatmask prokaryotes (no reason to change this), 1 = yes, 0 = no
softmask=1 #use soft-masking rather than hard-masking in BLAST (i.e. seg and dust filtering)

Step3: 配置"Gene Prediction",例如SNAP, 同时将"est2genome"和"protein2genome"设置为0

#-----Gene Prediction
snaphmm=snap.hmm #SNAP HMM file
gmhmm= #GeneMark HMM file
augustus_species= #Augustus gene prediction species model
# 略过其他参数
est2genome=0 #infer gene predictions directly from ESTs, 1 = yes, 0 = no
protein2genome=0 #infer predictions from protein homology, 1 = yes, 0 = no
# 略过其他参数

会跳过exonerate步骤,直接从snap预测开始。

方法3

我们还可以通过设置est2gff, protein_gffrm_gff,来避免重复序列屏蔽和BLAST+Exonerate运算

Step1: 从之前的MAKER输出的GFF文件种提取EST/Protein/Repeat的位置信息

# transcript alignment
awk '{ if ($2 ~ "est") print $0 }' round1.gff > est.gff
# protein alignments
awk '{ if ($2 == "protein2genome") print $0 }' round1.gff > protein2genome.gff
# repeat alignments
awk '{ if ($2 ~ "repeat") print $0 }' round1.gff > repeats.gff

Step2: 修改EST Evidence / rotein Homology Evidence /Repeat Masking里的配置参数

#-----EST Evidence (for best results provide a file for at least one)
est= #set of ESTs or assembled mRNA-seq in fasta format
altest= #EST/cDNA sequence file in fasta format from an alternate organismest_gff=est.gff #aligned ESTs or mRNA-seq from an external GFF3 file
est_gff=./est.gff
altest_gff= #aligned ESTs from a closly relate species in GFF3 format

#-----Protein Homology Evidence (for best results provide a file for at least one)
protein=  #protein sequence file in fasta format (i.e. from mutiple organisms)
protein_gff=protein2genome.gff  #aligned protein homology evidence from an external GFF3 file

#-----Repeat Masking (leave values blank to skip repeat masking)
model_org= #select a model organism for RepBase masking in RepeatMasker
rmlib= #provide an organism specific repeat library in fasta format for RepeatMasker
repeat_protein= #provide a fasta file of transposable element proteins for RepeatRunner
rm_gff=repeats.gff #pre-identified repeat elements from an external GFF3 file
prok_rm=0 #forces MAKER to repeatmask prokaryotes (no reason to change this), 1 = yes, 0 = no
softmask=1 #use soft-masking rather than hard-masking in BLAST (i.e. seg and dust filtering)

Step3: 配置"Gene Prediction",例如SNAP, 同时将"est2genome"和"protein2genome"设置为0

#-----Gene Prediction
snaphmm=snap.hmm #SNAP HMM file
gmhmm= #GeneMark HMM file
augustus_species= #Augustus gene prediction species model
# 略过其他参数
est2genome=0 #infer gene predictions directly from ESTs, 1 = yes, 0 = no
protein2genome=0 #infer predictions from protein homology, 1 = yes, 0 = no
# 略过其他参数

同样也会跳过exonerate步骤,直接从snap开始。

结果比较

对于这三种方法,从运行日志中看,三者都会跳过重复序列屏蔽,将EST和蛋白序列回帖到参考基因组的步骤,然而最终预测的基因数却不一致。

分析方法2和方法1的输出GFF文件时,发现方法2输出包括exonerate_protein2genome-geneexonerate_est2genome-gene。推测其原因在第二种方法的model_pass, pred_passs参数在设置为1时会使用之前est2genome和protein2genome输出的基因模型,而由于模型本身就来自于EST和Protein,就变成自我验证,于是输出结果就变多了。当设置model_pass, pred_passs参数为0时,最终保证方法2和方法1输出结果一致。

之后设置model_pass, pred_passs参数为0,然后比较方法1,方法2和方法3输出的GFF。我发现方法1和方法2的第二列信息完全相同,是blastn, blastx, est2genome, maker, protein2genome, repeatmasker, snap_masked, 而方法3的第二列为est_gff:est2genome, maker, protein_gff:protein2genome, repeat_gff:repeatmasker, snap_masked. 目前只能推测是MAKER对这些证据使用方式不同引起了最终输出结果的差异,但具体的原理我没有分析清楚,不过不妨碍使用。

最后,三种方法使用优先级分别是方法1 > 方法2 > 方法3,其中方法2要注意设置model_pass, pred_passs的设置。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,456评论 5 477
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,370评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 150,337评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,583评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,596评论 5 365
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,572评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,936评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,595评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,850评论 1 297
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,601评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,685评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,371评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,951评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,934评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,167评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 43,636评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,411评论 2 342