首先是安装这个流程,试了一下可以使用conda进行安装
conda search evidencemodeler
conda install evidencemodeler
安装好以后很多perl脚本是在 anaconda3/envs/EVM/opt/evidencemodeler-2.1.0/EvmUtils/
这个目录下
学习这个流程的参考链接
1、https://www.zhouxiaozhao.cn/2020/11/26/2020-11-26-annotion(4)/
基因结构注释(4):整合预测结果2、使用MAKER进行注释: 训练SNAP基因模型 //www.greatytc.com/p/8c378421af12
3、基因组结构注释 //www.greatytc.com/p/2cfc7638663d 这个写的很详细
我这里的数据就使用拟南芥的一条染色体,这个数据来源于论文
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、生物信息学入门学习资料及自己的学习笔记!