使用MAKER进行基因组注释

1. 实验目的

在本实验中,我们将使用MAKER进行基因组注释,它的工作原理是将尽可能多的证据沿着基因组序列排列,然后协调所有这些信号,以确定可能的基因结构。证据可以是来自同一(或密切相关)生物体的转录物或蛋白质序列。这些序列可以来自公共数据库或来自自己的实验数据(例如RNASeq实验的转录组组装)。Maker还能够考虑到重复的元素。

Maker使用ab-initio预测器(如Augustus或SNAP)来改进其预测:这些软件工具能够通过使用统计模型分析基因组序列来进行基因结构预测。

在本实验中,将学习如何执行基因组注释,以及如何评估其质量,学习训练从头算预测器产生良好结果的重要步骤,最后,将学习如何使用IGV来可视化结果。

参考:https://training.galaxyproject.org/training-material/topics/genome-annotation/tutorials/annotation-with-maker/tutorial.html

2. 实验准备

2.1 实验平台

Linux JSvr01 3.10.0-1160.80.1.el7.x86_64 #1 SMP Tue Nov 8 15:48:59 UTC 2022 x86_64 x86_64 x86_64 GNU/Linux

2.2 数据简述

数据 数据来源
maker_tutorial.tgz /work/cb-data/annotation2019
pyu_rnd1.maker.output.tgz /work/cb-data/annotation2019

2.3 软件配置

软件 版本 来源
FileZilla 3.65.0 https://filezilla-project.org/
IGV 2.6.3 https://data.broadinstitute.org/igv/projects/downloads/2.6/
BuildDatabase 2.0.1 /work/bin/RepeatModeler-2.0.1/
RepeatModeler 2.0.1 /work/bin/RepeatModeler-2.0.1/
maker 2.31.10 /work/bin/maker
mpiexec 1.5 /work/bin/maker/exe/mpich2/bin

3. 实验内容

在本实验中,将学习如何执行基因组注释,以及如何评估其质量,学习训练从头算预测器产生良好结果的重要步骤,最后,将学习如何使用IGV来可视化结果。

对于第一轮,我们将Maker配置为仅通过将EST和蛋白质与基因组对齐来构建基因模型。这将生成一个初稿注释,我们将在接下来的步骤中对其进行改进。

第二轮,我们使用SNAP进行从头算预测因子首次训练。在 Maker 中使用它之前,我们需要使用我们在前面步骤中生成的第一个注释草稿来训练它们。只有最好的评分基因(即那些有最强证据的基因)才会被保留来训练预测因子。“Ab-initio”意味着这些预测因子能够仅根据基因组的序列和物种特异性统计模型来预测基因组中基因的结构。他们不使用任何证据(例如EST或蛋白质)来预测基因,但他们需要用一组已经预测的基因进行训练。Maker能够使用EST和蛋白质证据,并将它们与几个从头预测因子的结果相结合,以预测共识基因模型。它允许检测没有EST或蛋白质对齐的区域的基因,还可以在有EST和/或蛋白质证据以及从头预测的区域改进基因结构。

第三轮,我们需要运行新一轮的Maker。由于证据在第一次运行时已经在基因组上对齐,我们可以按原样重用这些对齐。 但这一次,启用从头的基因预测,禁用直接从所有 EST 和蛋白质推断基因预测:现在我们希望 Maker 通过协调证据比对从头基因预测来推断基因预测。为了获得更好的结果,我们将对 SNAP进行第二次训练,然后第三次(最后一次)运行 Maker。

3.1 Prepare working directory

准备文件,创建工作目录。

将“ /work/cb-data/annotation2019/” 目录下的 data 文件拷贝到“ ~/$USER” 目录下,并解压。

1
2
3
4
5
6
7
mkdir $USER
mkdir $USER/tmp
cd $USER
cp ../annotation2019/maker_tutorial.tgz ./
tar -zxf maker_tutorial.tgz
cd maker_tutorial
ls -1

3.2 Maker round 1 - Map known genes to the genome

将已知基因映射到基因组中,包括:重复掩蔽和将已知的转录组/蛋白质序列与基因组对齐 。

(1) 构建自定义的重复数据库

这个步骤对于这个练习来说是可选的,因为它是一个非常小的基因组,不需要重复掩蔽就可以了。 当在做一个真实的项目时,可以从 RepBase,需要许可下载一个数据库,也可以用自己的基因组序列构建一个自定义的重复数据库。

RepeatModeler 是一个用于构建自定义数据库的软件。这里提供了用于构建重复数据库的命令。

1
2
3
4
5
6
7
8
cd ~/$USER/maker_tutorial/example_02_abinitio
# 可以删除路径环境变量
export PATH=$(echo $PATH | tr ':' '\n' | grep -v '/work/bin/RepeatModeler-2.0.1' | tr '\n' ':' | sed 's/:$//')
# 重新定义路径
export PATH=/work/bin/RepeatModeler-2.0.1:$PATH
BuildDatabase -name pyu pyu_contig.fasta
RepeatModeler -pa 4 -database pyu >& repeatmodeler.log
nano repeatmodeler.log

运行结束后获得pyu-families.fa,将其提供给maker_opts.ctl文件的“rmlib= ”选项

(2) 设置环境运行 Maker 并创建 Maker 控制文件

Maker 中的每一步都是由 Maker 控制文件指定的。 “ maker -CTL” 命令将创建三个控制文件: maker_bopts.ctl, maker_exe.ctl, maker_opts.ctl.by

1
2
3
4
export PATH=/work/bin/maker/bin:/work/bin/RepeatMasker:/work/bin/maker/exe/snap:$PATH
export ZOE=/work/bin/maker/exe/snap/Zoe
cd ~/$USER/maker_tutorial/example_02_abinitio
maker -CTL

(3) 修正控制文件 maker_opts.ctl

将修改后的文件放在相同的目录“example_02_abinitio” 中。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#  在 maker_opts.ctl修改以下变量,分为两部分:

## 运行RepeatMasker的配置
model_org=simple # tells RepeatMasker to mask the low complexity sequence (e.g.“AAAAAAAAAAAAA”)
rmlib= # tells RepeatMasker not to mask repeat sequences like transposon elements. If you have a repeat fasta file (e.g. output fromRepeatModeler) that you need to mask, put the fasta file name next to “rmlib=”
softmask=1 # tells RepeatMasker to do soft-masking which converts repeats to lower case, instead of hard-masking which converts repeats to “N”. "Soft-masking" is important so that short repeat sequences within genes can still be annotated as part of gene.

## 匹配文件pyu_est.fasta中转录序列与文件sp_protein.fasta中蛋白质序列到基因组以及推断证据支持基因模型中
genome=pyu_contig.fasta
est=pyu_est.fasta
protein=sp_protein.fasta # specify the fasta file names of the EST and protein sequences. In general, the EST sequence file contains the assembled transcriptome from RNA-seq data. The protein sequence file include proteins from closely related species or swiss-prot. If you have multiple protein or EST files, separate file names with ",".
est2genome=1
protein2genome=1 # tell MAKER to align the transcript sequences from the pyu_est.fasta file and protein sequences from the sp_protein.fasta file to the genome. These two files are used to define evidence supported gene model.
TMP=~/$USER/tmp #important for big genome, as the default /tmp is too small
(4) 利用screen执行重复序列屏蔽

重复掩蔽和对齐。 这一步需要一个小时。在screen中运行。在命令中:“ mpiexec -n 2” 意味着你将使用 MPI 并行 Maker,并且一次使用两个线程。

1
2
3
4
5
6
screen
export PATH=/work/bin/maker/bin:/work/bin/RepeatMasker:/work/bin/maker/exe/snap:$PATH
export ZOE=/work/bin/maker/exe/snap/Zoe

cd ~/$USER/maker_tutorial/example_02_abinitio
/work/bin/maker/exe/mpich2/bin/mpiexec -n 2 maker -base pyu_rnd1 >& log1 &

Maker is now finished!!!

3.3 Maker round 2 - Gene prediction using SNAP

使用 SNAP 进行基因预测

(1) 训练一个SNAP基因模型

SNAP是从基因组进行从头算基因预测的软件。使用 SNAP做基因预测,将首先使用上一个生成的对齐结果训练 SNAP 模型。

1
2
3
cd ~/$USER/maker_tutorial/example_02_abinitio
cp ~/annotation2019/pyu_rnd1.maker.output.tgz ./
tar xvfz pyu_rnd1.maker.output.tgz

设置环境变量,如果是新session

1
2
export PATH=/work/bin/maker/bin:/work/bin/RepeatMasker:/work/bin/maker/exe/snap:$PATH
export ZOE=/work/bin/maker/exe/snap/Zoe

接下的命令将MAKER round 1的结果传递到输入文件,构造一个SNAP模型

1
2
3
4
mkdir snap1
cd snap1
gff3_merge -d ../pyu_rnd1.maker.output/pyu_rnd1_master_datastore_index.log
maker2zff -l 50 -x 0.5 pyu_rnd1.all.gff # specify that only gene models with AED score>0.5 and protein length>50 are used for building models.

生成了三个文件:pyu_rnd1.all.gff & genome.ann & genome.dna

maker2zff生成一个ZFF格式文件(genome.ann)和一个FASTA格式文件(genome.dna),过滤用于再次训练的高置信度基因,共有7个选项:

1
2
3
4
5
6
7
8
-c 由EST/mRNA-Seq比对确定的剪接位点的比例,默认0.5
-e 与EST/mRNA-Seq比对重叠的外显子的比例,默认0.5
-o 和任何证据(EST或者蛋白)重叠的外显子的比例,默认0.5
-a 从头预测证实的剪接位点的比例,默认0
-t 和从头预测结果重叠的外显子的比例,默认0
-l mRNA翻译的蛋白质序列的最短长度
-x 最大AED值,默认0.5
-n 不过滤

AED值:maker2使用注释编辑距离(AED)来评估基因组注释的准确性,AED是一个介于 0 和 1 之间的数字,衡量注释与支持它的evidence的拟合优度,0 表示与可用证据完全一致,1 表示缺乏对注释基因模型的支持

现在将运行以下命令来训练SNAP。训练SNAP的基本步骤是:首先对输入的基因模型进行过滤,然后立即捕获周围的基因组序列每个模型轨迹,最后使用这些捕获的片段生成HMM。

1
2
3
4
5
6
fathom -categorize 1000 genome.ann genome.dna
fathom -export 1000 -plus uni.ann uni.dna
forge export.ann export.dna
hmm-assembler.pl pyu . > ../pyu1.hmm
mv pyu_rnd1.all.gff ../
cd ..

在此之后,将在 example_02_abinitio 目录中找到两个新文件:pyu_rnd1.all.gff:第一轮的gff文件,它是基于证据的基因。 pyu1.hmm:一个从基于证据的基因训练出来的隐藏马尔科夫模型。

(2) 用SNAP去预测基因

直接修改 maker_opts。之前修改过的 CTL文件。在此之前,保存一个maker_opts 的备份副本 CTL用于第一轮。

1
cp maker_opts.ctl maker_opts.ctl_backup_rnd1

修改文件maker_opts.ctl 中一些变量的值:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
maker_gff= pyu_rnd1.all.gff
est_pass=1 # use est alignment from round 1
protein_pass=1 #use protein alignment from round 1
rm_pass=1 # use repeats in the gff file
snaphmm=pyu1.hmm
est= # remove est file, do not run EST blast again
protein= # remove protein file, do not run blast again
model_org= #remove repeat mask model, so not running RM again
rmlib= # not running repeat masking again
repeat_protein= #not running repeat masking again
est2genome=0 # do not do EST evidence based gene model
protein2genome=0 # do not do protein based gene model.
pred_stats=1 #report AED stats
alt_splice=0 # 0: keep one isoform per gene; 1: identify splicing variants of the same gene
keep_preds=1 # keep genes even without evidence support, set to 0 if no

用新的控制文件运行 maker。这个步骤需要几分钟。 (一个真正的项目可能需要几个小时才能完成)。将使用选项“ -base pyu_rnd2” ,这样结果将被写入一个新目录“ pyu_rnd2” 。

1
/work/bin/maker/exe/mpich2/bin/mpiexec -n 2 maker -base pyu_rnd2 >& log2 &

Maker is now finished!!!

3.4 Maker round 3 - Retrain SNAP model and do another round of SNAP gene prediction

重新训练 SNAP 模型,并进行另一轮 SNAP 基因预测,可能需要运行两到三轮 SNAP。所以你要再次重复第二部分。确保将 snap1 替换为 snap2,这样就不会覆盖上一轮。

(1) 首先训练一个新的SNAP模型
1
2
3
4
5
6
7
8
9
10
mkdir snap2
cd snap2
gff3_merge -d ../pyu_rnd2.maker.output/pyu_rnd2_master_datastore_index.log
maker2zff -l 50 -x 0.5 pyu_rnd2.all.gff
fathom -categorize 1000 genome.ann genome.dna
fathom -export 1000 -plus uni.ann uni.dna
forge export.ann export.dna
hmm-assembler.pl pyu . > ../pyu2.hmm
mv pyu_rnd2.all.gff ..
cd ..
(2) 使用SNAP预测基因

为第二轮的maker_opts.ctl 保存副本:

1
cp maker_opts.ctl maker_opts.ctl_backup_rnd2

编辑第三轮的maker_opts.ctl

1
2
maker_gff=pyu_rnd2.all.gff
snaphmm=pyu2.hmm
1
/work/bin/maker/exe/mpich2/bin/mpiexec -n 2 maker -base pyu_rnd3 >& log3 &

使用下面的命令创建最终的合并 gff 文件。 " -n "选项将产生一个没有基因组序列的 gff 文件:

1
2
gff3_merge -n -d pyu_rnd3.maker.output/pyu_rnd3_master_datastore_index.log>pyu_rnd3.noseq.gff
fasta_merge -d pyu_rnd3.maker.output/pyu_rnd3_master_datastore_index.log

在这之后,会得到一个新的 gff3 文件:pyu_rnd3.noseq.gff,以及蛋白质和转录 fasta 文件。

(3) 生成AED图
1
2
/work/bin/maker/AED_cdf_generator.pl -b 0.025 pyu_rnd2.all.gff > AED_rnd2
/work/bin/maker/AED_cdf_generator.pl -b 0.025 pyu_rnd3.all.gff > AED_rnd3

可以使用 Excel 绘制 AED_rnd2 和 AED_rnd3 文件的第二列,并使用第一列作为 x轴值。 x 轴标签为“AED” , y 轴标签为“Cumulative Fraction of Annotations”

AED值:maker2使用注释编辑距离(AED)来评估基因组注释的准确性,AED是一个介于 0 和 1 之间的数字,衡量注释与支持它的evidence的拟合优度,0 表示与可用证据完全一致,1 表示缺乏对注释基因模型的支持。结合Cumulative Fraction of Annotations可进行质量评估。

3.5 Visualize the gff file in IGV

以将 gff 文件与fasta文件一起加载到 IGV 中。 用FileZilla下载pyu_contig.fasta和pyu_rnd1.all.gff,pyu_rnd2.all.gff,pyu_rnd3.all.gff文件。导入IGV进行可视化:

例如scf1117875582023:644,587-656,318,在第一轮结束时,Maker预测了该地区的部分短基因模型。 在第二轮之后,Maker能够预测该区域的其他基因模型。这意味着 Maker 主要使用 SNAP 的基因预测来构建此基因模型。 第三轮后,基于第二轮预测基因。训练SNAP可以完善基因结构并完善该区域发现的基因结构。

4. 实验总结

4.1 实验结论

第二个注释与前一个注释相比,具有更多的完整单拷贝,以及更多的基因使用从头预测因子可以在EST或蛋白质比对不足以预测基因的区域找到更多的基因。

第三轮训练中,Maker能够预测更复杂的基因,例如合并一些事先被认为是独立的基因,不需要超过两轮训练即可从从头预测器获得最佳结果。

4.2 实验收获

学习了如何使用 Maker 注释真核基因组、训练SNAP并预测基因、如何评估注释的质量以及如何使用 IGV对结果进行可视化。


使用MAKER进行基因组注释
http://horizongazer.github.io/2023/10/20/使用MAKER进行基因组注释/
作者
HorizonGazer
发布于
2023年10月20日
许可协议