EVidenceModeler(EVM)流程做基因组注释简单小例子

首先是安装这个流程,试了一下可以使用conda进行安装

conda search evidencemodeler
image.png
conda install evidencemodeler

安装好以后很多perl脚本是在 anaconda3/envs/EVM/opt/evidencemodeler-2.1.0/EvmUtils/这个目录下

学习这个流程的参考链接

我这里的数据就使用拟南芥的一条染色体,这个数据来源于论文

Chromosome-level assemblies of multiple Arabidopsis genomes reveal hotspots of rearrangements with altered evolutionary dynamics

C24.chr.all.v2.0.fasta的一号染色体

evm 流程第一步用到的命令是

time ~/anaconda3/envs/EVM/opt/evidencemodeler-2.1.0/EvmUtils/evidence_modeler.pl \
--genome ../../repeat/chr1.fa.masked \
--weights /data/myan/raw_data/practice/pan.genome/at.nc/version2.5.2019-10-09/genome.annotation/braker2/evm/weights.txt \
--gene_predictions evm_abinitio.gff3 \
--protein_alignments evm_pro.gff3 \
--transcript_alignments transcripts.fasta.transdecoder.genome.gff3 > evm.out

这里需要的输入数据

  • 染色体序列

  • weights.txt

这里的内容

PROTEIN genomeThreader  4
TRANSCRIPT transdecoder 8
ABINITIO_PREDICTION Augustus    1
ABINITIO_PREDICTION GeneMark.hmm    1
  • evm_abinitio.gff3

这个是从头预测的结果用到的是augustus和genemark ,我这里直接用braker2这个流程先跑一遍,就可以得到这两个结果文件 将这两个结果文件合并成了evm_abinitio.gff3

  • evm_pro.gff3 这个是基于同源蛋白的结果使用的是gth这个程序

  • transcripts.fasta.transdecoder.genome.gff3 这个是转录组的数据

接下来介绍如何得到这些输入文件

首先是对重复序列进行屏蔽

RepeatMasker -species Arabidopsis -pa 12 -xsmall -dir repeat chr1.fa

这里拟南芥可以直接指定物种,如果自己研究的物种不是很常见的,需要构建这个重复序列的库

基于同源蛋白的注释

evm_pro.gff3 这个相对简单,这里需要有一个同源蛋白文件和需要注释的染色体基因组,得到gff文件后需要用evm这个流程里的脚本对格式进行转换,gth这个软件的安装,如果braker2这个软件安装好是可以直接用的

~/anaconda3/envs/braker2/bin/gth -gff3out -intermediate -duplicatecheck seq -protein proteins.fa -translationtable 1 -genomic chr1.fa.mask -o homo_protein.gff3 -force
~/anaconda3/envs/EVM/opt/evidencemodeler-2.1.0/EvmUtils/misc/genomeThreader_to_evm_gff3.pl homo_protein.gff3 > evm_pro.gff3

基于转录组数据

transcripts.fasta.transdecoder.genome.gff3 这个文件用到转录组数据

hisat2-build ../repeat/chr1.fa.masked chr1 -p 8
hisat2 -p 8 --dta -x chr1 -1 SRR4420296_R1.fastq.gz -2 SRR4420296_R2.fastq.gz -S SRR4420296.sam
samtools sort -@ 8 -O bam -o SRR4420296.bam SRR4420296.sam

stringtie -p 12 -o SRR4420296.gtf SRR4420296.bam
gtf_genome_to_cdna_fasta.pl SRR4420296.gtf ../repeat/chr1.fa.masked > transcripts.fasta
gtf_to_alignment_gff3.pl SRR4420296.gtf > transcripts.gff3
TransDecoder.LongOrfs -t transcripts.fasta
TransDecoder.Predict -t transcripts.fasta
cdna_alignment_orf_to_genome_orf.pl transcripts.fasta.transdecoder.gff3 transcripts.gff3 transcripts.fasta > transcripts.fasta.transdecoder.genome.gff3

从头预测

braker.pl --cores 48 --species=At01 --genome=../repeat/chr1.fa.masked --softmasking --bam=SRR4420296.bam --gff3 --prot_seq=../proteins.fa --prg=gth

这个从头预测的时候也可以指定蛋白证据,但不知道为什么这个生成的结果文件里没有单独的gth结果,不知道是不是需要单独指定某个参数

这一步会生成 augustus.hints.gtf 和 GeneMark-ET/genemark.gtf

也需要做一个格式转换

~/anaconda3/envs/EVM/opt/evidencemodeler-2.1.0/EvmUtils/misc/augustus_GTF_to_EVM_GFF3.pl augustus.hints.gtf > evm_augustus.gff3

~/anaconda3/envs/EVM/opt/evidencemodeler-2.1.0/EvmUtils/misc/GeneMarkHMM_GTF_to_EVM_GFF3.pl GeneMark-ET/genemark.gtf > evm_genemark.gff3

cat evm_augustus.gff3 evm_genemark.gff3 > evm_abinitio.gff3

有了这些结果文件就可以运行开头提到的命令

time ~/anaconda3/envs/EVM/opt/evidencemodeler-2.1.0/EvmUtils/evidence_modeler.pl \
--genome ../../repeat/chr1.fa.masked \
--weights /data/myan/raw_data/practice/pan.genome/at.nc/version2.5.2019-10-09/genome.annotation/braker2/evm/weights.txt \
--gene_predictions evm_abinitio.gff3 \
--protein_alignments evm_pro.gff3 \
--transcript_alignments transcripts.fasta.transdecoder.genome.gff3 > evm.out

得到evm.out后怎么处理后面再来更新,因为程序还没有跑完,还没有拿到这个结果

这一步是可以并行的,这个并行怎么用还需要仔细研究一下

推文记录的是自己的学习笔记,内容可能会存在错误,请大家批判着看,欢迎大家指出其中的错误

示例数据和代码可以给推文点赞,然后点击在看,最后留言获取

欢迎大家关注我的公众号

小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享:1、R语言和python做数据分析和数据可视化的简单小例子;2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记;3、生物信息学入门学习资料及自己的学习笔记!

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

推荐阅读更多精彩内容