MAKER的使用

1. MAKER 简介

MAKER 是进行基因组结构注释的一整套流程。它综合了ab initio方法和 evidence 的方法进行基因组注释。使用起来非常方便快捷,一个命令能搞定全套基因组结构注释。

2. MAKER下载与安装

http://yandell.topaz.genetics.utah.edu/cgi-bin/maker_license.cgi页面提交信息并下载最新版本MAKER。

此外,安装 maker 需要先安装一些软件

必须安装:
Perl 5.8.0 or higher
BioPerl 1.6 or higher
WU-BLAST 2.0 or higher or NCBI-BLAST 2.2.X or higher 
RepeatMasker 3.1.6 or higher
    RepeatMasker requires a repeat library, available from Repbase. 
Exonerate 1.4 or higher

可选安装:
SNAP version 2009-02-03 or higher (for eukaryotic genomes).
Augustus 2.0 or higher (for eukaryotic genomes) 
GeneMark-ES (for eukaryotic genomes)
FGENESH 2.6 or higher - requires licence (for eukaryotic genomes)
GeneMarkS (for prokaryotic genomes) 
snoscan
MPICH2 (install MPICH2 with the --enable-sharedlibs flag)
安装perl模块
$ cpan -i BioPerl Bit::Vector DBD::SQLite DBI Error Error::Simple File::NFSLock File::Which forks forks::shared Inline Inline::C IO::All IO::Prompt PerlIO::gzip Perl::Unsafe::Signals Proc::ProcessTable Proc::Signal threads URI::Escape

安装mpich2
$ sudo yum install mpich2*
使用 MPI 运行 Maker 需要先运行服务 mpd。
正常运行 mpd 命令,需要正确的主机名。否则可能会遇到报错。需要修改 /etc/hosts 文件,将主机名的 IP 对应正确。

安装maker
$ wget http://yandell.topaz.genetics.utah.edu/maker_downloads/2D3B/22C6/BE99/41DBB7F9B1D181B2CDCE1E49DFE6/maker-2.31.8.tgz
$ tar zxf maker-2.31.8.tgz 
$ cd maker/src/
$ perl Build.PL 
$ ./Build install
$ echo 'PATH=$PATH:/opt/biosoft/maker/bin/' >> ~/.bashrc
$ source ~/.bashrc 

3. MAKER 的使用

3.1 准备 MAKER 的输入文件

MAKER 需要的输入文件如下:

1. genome.fasta
    基因组序列文件。
2. inchworm.fasta
    转录组组装结果文件。此文件内容代表转录子序列或 EST 序列。为了增加其准确性,我个人选用了 trinity 中 inchworm 的结果文件,而不是转录组最终的组装结果文件。
3. proteins.fasta
    临近物种的蛋白质序列。可以从 NCBI 的 Taxonomy 数据库中下载。
4. consensi.fa.classified
    适合本物种的重复序列数据库。该文件是一个fasta文件,其序列代表着本物种的高重复序列。该文件使用 RepeatModeler 制作。
5. Augustus hmm files [species]
    适合本物种的 Augustus 的 HMM 文件。该文件以文件夹形式存放在 /opt/biosoft/augustus-3.0.3/config/species/ 目录下。将 genome.fasta 和 inchworm.fasta 输入 PASA,得到用于 trainning 的基因,然后使用 Augustus 进行 training,得到 HMM 文件。
6. Snap hmm file [species.hmm]
    适合本物种的 SNAP 的 HMM 文件。该文件是一个单独的文件。将 genome.fasta 和 inchworm.fasta 输入 PASA,得到用于 trainning 的基因,然后使用 SNAP 进行 training,得到 HMM 文件。
7. Genemark_es hmm file [es.mod]
    适合本物种的 Genemark_es 的 HMM 文件。将 genome.fasta 输入 genemark_es 软件进行自我训练得到的 HMM 文件。
8. rRNA.fa
    rRNA 的序列文件。该文件使用 RNAmmer 对基因组进行分析得到。

得到上述输入文件后,使用 MAKER 命令得到初始化的 3 个配置文件。

$ maker -CTL

得到3个配置文件: maker_exe.ctl, maker_bopts.ctl 和 maker_opts.ctl。
maker_exe.ctl: 配置 Maker 需要调用的程序的路径
maker_bopts.ctl: 配置 blast 的各种阈值
maker_opts.ctl: MAKER 的主要配置文件,配置输入文件的路径,以及 MAKER 的使用流程。

3.2 使用 MAKER 进行基因组注释

先修改 maker_opts.ctl 配置文件,需要修改的内容如下:

genome=genome.fasta            # 基因组序列文件
est=inchworm.fasta             # EST 序列文件
protein=proteins.fasta         # 临近物种的蛋白质序列文件
model_org=fungi                # 选择repebase数据库的物种,对真菌则用fungi,也可以使用 all
rmlib=consensi.fa.classified   # 适合本物种的重复序列数据库
snaphmm=species.hmm            # snap traing 的 HMM 文件
gmhmm=es.mod                   # genemark_es 的 HMM 文件
augustus_species=species       # augustus 的 HMM 文件
est2genome=1                   # 使用 EST 序列进行基因预测
protein2genome=1               # 使用蛋白质序列进行基因预测
trna=1                         # 使用 tRNAscan 进行 tRNA 预测
snoscan_rrna=rRNA.fa           # 使用 Snoscan 利用 rRNA 序列进行 C/D box 类型的 snoRNAs 预测,不过我加入了这个信息后,maker运行不成功,不知道为什么。
correct_est_fusion=1           # 进行融合 EST 的校正

并行化运行 maker。由于基因组序列比较大,MAKER 的计算量大,使用并行化能极大节约时间。

$ mpd &
$ maker                  # 不并行化运行
$ mpiexec -n 20 maker    # 使用20个线程并行计算

JBrowse 的下载和安装

1. JBrowse 简介

2. JBrowse 下载和安装

下载 JBrowse 安装包,并安装:

$ wget http://jbrowse.org/releases/JBrowse-1.11.5.zip
$ wget http://jbrowse.org/releases/JBrowse-1.11.5-dev.zip
$ unzip JBrowse-1.11.5.zip -d /opt/biosoft
$ unzip JBrowse-1.11.5-dev.zip -d /opt/biosoft
$ cd /opt/biosoft/JBrowse-1.11.5
$ sudo ./setup.sh
$ echo 'PATH=$PATH:/opt/biosoft/JBrowse-1.11.5' >> ~/.bashrc

配置 apache:

# echo 'Alias /JBrowse "/opt/biosoft/JBrowse-1.11.5/"
<Directory "/opt/biosoft/JBrowse-1.11.5/">
    Options MultiViews ExecCGI Indexes FollowSymlinks
    AllowOverride None
    Order allow,deny
    Allow from all
</Directory>' > /etc/httpd/conf.d/JBrowse.conf
# /etc/init.d/httpd restart

打开软件自带的示例 Volvox JBrowse

3. JBrowse 配置

基因预测流程

1. 获取准确的转录子序列

利用 RNA-seq 数据进行转录组 de novo 组装。 推荐使用 trinity 软件进行有基因指导的 de novo 组装,只进行 inchworm 组装,这样能获得比较准确的转录组序列。

2. 获取完整的基因模型

使用 PASA 将上一步骤得到的转录子序列比对到基因组上。根据比对结果,提取完整并且长度较长的基因模型。

3. HMM training

根据上一步骤的基因模型进行 AUGUSTUS 和 SNAP 的 training,得到这 2 个基因预测软件的 HMM 参数文件。

4. 由 RNA-seq 数据得到 hints

对基因组进行重复序列分析,对 DNA 转座子和逆转录转座子进行 hardmask (即对该区域使用 N 屏蔽),对 low-complexity 区域进行 softmask(即将该区域碱基换成小写)。 将 RNA-seq 数据比对到屏蔽了重复序列的基因组,得到 hints 信息。

5. 获取低表达区域来自其它物种的保守蛋白

在屏蔽重复序列的基础上,对上一步骤得到 hints 区域进行hardmask。 从 NCBI 的 Taxonomy 数据库下载蛋白质序列,并将这些序列比对到屏蔽了重复序列和基因表达区域的基因组上,设定阈值筛选出保守蛋白。这些蛋白比对到了表达量低的基因区,有利于这些基因的预测。

6. 使用 MAKER 进行第一遍基因预测

使用 MAKER 进行基因预测,输入文件有:

基因组序列: 基因组装的结果文件
重复序列文库: 使用 ReapeterModeler 得到
HMM 文件: AUGUSTUS, SNPA 和 GeneMark-ES 的 HMM 文件
转录子序列: PASA 的结果文件
蛋白质序列: 筛选出来的蛋白序列

7. 使用 MAKER 进行第二遍基因预测

根据上一步的预测结果,提取 AED 值较低的基因模型,进行第二遍 HMM training,然后再次使用 MAKER 进行基因预测。

8. 进行基因模型的人工校正

Augustus Training and Prediction

1. Augustus training

首选,需要有至少 200 个完整基因模型的数据。 例如: 使用 genome-guided 方法进行 trinity 有参考基因组的 de novo 组装;再使用 PASA 将组装出来的 inchworm 序列比对到基因组; 再提取完整基因模型数据,得到文件 trainingSet_CompleteBest.gff3 。

然后,将 GFF3 文件转换为 GeneBank 文件:

$ /opt/biosoft/augustus-3.0.3/scripts/gff2gbSmallDNA.pl trainingSet_CompleteBest.gff3 ../../genome.fasta 50 trainingSetComplete.gb

使用 genebank 格式的文件预先进行一次 Augustus training,一般会得到一些错误提示
$ etraining --species=generic --stopCodonExcludedFromCDS=false trainingSetComplete.gb 2> train.err

根据错误提示,提取出有错误的基因模型
$ cat train.err | perl -ne 'print "$1\n" if /in sequence (\S+):/' > badlist

去除错误的基因模型,得到能用于 training 的基因模型
$ /opt/biosoft/augustus-3.0.3/scripts/filterGenes.pl badlist trainingSetComplete.gb > training.gb

在进行 Augustus training 之前,最好保证这些基因模型两两之间在氨基酸水平的 identity < 70% :

先获取这些基因模型的 Proteins 序列
$ /opt/biosoft/augustus-3.0.3/scripts/gbSmallDNA2gff.pl training.gb > training.gff2
$ perl -ne ‘print “$1\n” if /gene_id \”(\S+?)\”/’ training.gff2 | uniq > trainSet.lst
$ perl extract_genes.pl trainingSet_CompleteBest.gff3 trainSet.lst > training.gff3
$ /opt/biosoft/EVM_r2012-06-25/EvmUtils/gff3_file_to_proteins.pl training.gff3 genome.fasta prot > training.protein.fasta
$ perl -p -i -e ‘s/(>\S+).*/$1/’ training.protein.fasta
$ perl -p -i -e ‘s/\*//’ training.protein.fasta

对这些 proteins 序列构建 blast 数据库,并将 proteins 序列比对到此数据库
$ makeblastdb -in training.protein.fasta -dbtype prot -title training.protein.fasta -parse_seqids -out training.protein.fasta
$ blastp -db training.protein.fasta -query training.protein.fasta -out training.protein.fasta.out -evalue 1e-5 -outfmt 6 -num_threads 8

提取 blast 结果中 identity >= 70% 的比对信息,identity 高的 proteins 序列仅保留一条
$ grep -v -P “\t100.00\t” training.protein.fasta.out | perl -ne ‘split; print if $_[2] >= 70’ > blast.out
$ perl delete_high_identity_proteins_in_training_gff3.pl training.protein.fasta blast.out training.gff3 > training.new.gff3

获得去除了冗余的基因模型
$ /opt/biosoft/augustus-3.0.3/scripts/gff2gbSmallDNA.pl training.new.gff3 genome.fasta 50 training.gb

进行 Augustus Training without CRF

初始化本物种的 HMM 文件
$ /opt/biosoft/augustus-3.0.3/scripts/new_species.pl --species=my_species

如果有 RNA-Seq 数据,获得能用于 training 的基因模型一般会有好几千个。我们将基因模型随机分成两份: 第一份 300 个基因,用于检测 training 的精确性; 另外一份含有更多的基因,用于进行 Augustus training。
$ /opt/biosoft/augustus-3.0.3/scripts/randomSplit.pl training.gb 300

使用大份的 traing.gb.train 进行 Augustus Training
$ etraining --species=my_species training.gb.train >train.out

根据输出结果 train.out 来修正参数文件 species_parameters.cfg 中终止密码子的频率
$ tag=$(perl -ne 'print "$1" if /tag:\s+\d+\s+\((\S+)\)/' train.out)
$ perl -p -i -e "s#/Constant/amberprob.*#/Constant/amberprob                   $tag#" /opt/biosoft/augustus-3.0.3/config/species/lentinula_edodes/lentinula_edodes_parameters.cfg
$ taa=$(perl -ne 'print "$1" if /taa:\s+\d+\s+\((\S+)\)/' train.out)
$ perl -p -i -e "s#/Constant/ochreprob.*#/Constant/ochreprob                   $taa#" /opt/biosoft/augustus-3.0.3/config/species/lentinula_edodes/lentinula_edodes_parameters.cfg
$ tga=$(perl -ne 'print "$1" if /tga:\s+\d+\s+\((\S+)\)/' train.out)
$ perl -p -i -e "s#/Constant/opalprob.*#/Constant/opalprob                    $tga#" /opt/biosoft/augustus-3.0.3/config/species/lentinula_edodes/lentinula_edodes_parameters.cfg

根据 training 的结果,进行第一遍精确性评估
$ augustus --species=my_species training.gb.test > test.1.out

再将 training.gb.train 分成两份: 第一份 800 个基因,剩下基因为另外一份。这两份基因都用于 Optimizing meta parameters of AUGUSTUS
$ /opt/biosoft/augustus-3.0.3/scripts/randomSplit.pl training.gb.train 800

使用上面的 2 份基因模型,进行 Augustus training 的优化
$ /opt/biosoft/augustus-3.0.3/scripts/optimize_augustus.pl --species=my_species --rounds=5 --cpus=16 --kfold=16 training.gb.train.test --onlytrain=training.gb.onlytrain --metapars=/opt/biosoft/augustus-3.0.3/config/species/my_species/my_species_metapars.cfg > optimize.out
按如上参数,则程序会对 my_species_metapars.cfg 文件中的 28 个参数进行优化,总共优化 5 轮或有一轮找不到可优化的参数为止。每进行一个参数的优化时: 将 training.gb.train.test 文件中 800 个基因随机分成 16 等份,取其中 15 等份和 training.gb.onlytrain 中的基因模型一起进行 training,剩下的 1 等份用于精确行评估; 要对每个等份都进行一次精确性评估;使用 16 个 CPU 对 16 个等份并行进行 training 和 精确性评估,得到平均的精确值;优化的每个参数会分别 3~4 个值,每个值得到一个 training 的精确值,对参数的多个设定值进行比较,找到最佳的值。
此外, training 的精确值的算法: accuracy value = (3*nucleotide_sensitivity + 2*nucleotide_specificity + 4*exon_sensitivity + 3*exon_specificity + 2*gene_sensitivity + 1*gene_specificity) / 15 。

优化参数完毕后,需要再此使用 training.gb.train 进行一遍 training
$ etraining --species=my_species training.gb.train

再进行第二遍的精确性评估,一般进行优化后,精确性会有较大幅度的提高
$ augustus --species=my_species training.gb.test > test.2.withoutCRF.aout

进行 Augustus Training with CRF(Conditional Random Field)

在进行 Training with CRF 之前,先备份非 CRF 的参数文件
$ cd /opt/biosoft/augustus-3.0.3/config/species/species/
$ cp species_exon_probs.pbl lentinula_edodes_exon_probs.pbl.withoutCRF
$ cp species_igenic_probs.pbl lentinula_edodes_igenic_probs.pbl.withoutCRF
$ cp species_intron_probs.pbl lentinula_edodes_intron_probs.pbl.withoutCRF
$ cd -

在 training 的时候,加入 --CRF 参数,这样进行 training 比较耗时
$ etraining --species=my_species training.gb.train --CRF=1 1>train.CRF.out 2>train.CRF.err

再次进行精确行评估
$ augustus --species=my_species training.gb.test > test.2.CRF.out
比较 CRF 和 非 CRF 两种情况下的精确性。一般情况下,CRF training 的精确性要高些。若 CRF training 的精确行低些,则将备份的参数文件还原回去即可。

2. Training UTR parameters for Augustus

Training UTR 有利于利用 exonpart hint。 当使用 exonpart hint 时,对基因结构预测有显著提高。

首先,需要获得用于 training 的基因模型,这些基因模型需要同时带有 5’UTR 和 3’UTR 。

$ mkdir utr
$ cd utr

从用于 Augustus Training 的数据中提取同时带有 5'UTR 和 3'UTR 的基因模型。
$ perl -ne 'print "$2\t$1\n" if /.*\t(\S+UTR)\t.*ID=(\S+).utr/' ../training.new.gff3 | sort -u | perl -ne 'split; print "$_[0]\n" if ($g eq $_[0]); $g = $_[0];' > bothutr.lst
$ perl -e 'open IN, "bothutr.lst"; while () {chomp; $keep{$_}=1} $/="//\n"; while (<>) { if (/gene=\"(\S+?)\"/ && exists $keep{$1}) {print} }' ../training.gb.test > training.gb.test
$ perl -e 'open IN, "bothutr.lst"; while () {chomp; $keep{$_}=1} $/="//\n"; while (<>) { if (/gene=\"(\S+?)\"/ && exists $keep{$1}) {print} }' ../training.gb.train > training.gb.train

使用 traing.gb.train 中的基因模型进行 Training UTR parameters
$ /opt/biosoft/augustus-3.0.3/scripts/randomSplit.pl training.gb.train 400
$ mv training.gb.train.train training.gb.onlytrain
$ optimize_augustus.pl --species=my_species --cpus=16 --rounds=5 --kfold=16 --UTR=on --trainOnlyUtr=1 --onlytrain=training.gb.onlytrain --metapars=/opt/biosoft/augustus-3.0.3/config/species/lentinula_edodes/lentinula_edodes_metapars.utr.cfg training.gb.train.test > optimize.out

进行 without CRF 和 with CRF 的精确性比较,选取其中精确性较高的参数文件
$ etraining --species=species --UTR=on training.gb.train
$ augustus --species=species --UTR=on training.gb.test > test.withoutCRF.out
$ etraining --species=species --UTR=on training.gb.train --CRF=1
$ augustus --species=species --UTR=on training.gb.test > test.CRF.out

3. Creating hints from RNA-Seq data with Tophat

先使用 tophat 将 RNA-Seq 数据比对到屏蔽了重复序列的基因组,得到 bam 文件。对 bam 文件进行转换,得到 intron 和 exonpart hints。

得到 intron 的 hints 信息
$ /opt/biosoft/augustus-3.0.3/bin/bam2hints --in=tophat.bam --out=hints.intron.gff --maxgenelen=30000 --intronsonly

得到 exonpart 的 hints 信息
$ /opt/biosoft/augustus-3.0.3/bin/bam2wig tophat.bam tophat.wig
$ cat tophat.wig | /opt/biosoft/augustus-3.0.3/scripts/wig2hints.pl --width=10 --margin=10 --minthresh=2 --minscore=4 --prune=0.1 --src=W --type=ep --UCSC=unstranded.track --radius=4.5 --pri=4 --strand="." > hints.exonpart.gff

cat hints.intron.gff hints.exonpart.gff > hints.rnaseq.gff

4. Creating hints from Protiens with Exonerate

要先得到临近物种的 pritein 序列,然后使用 tblastn 将这些 protein 序列比对到基因组,再将相似性较高的序列使用 exonerate 比对到基因组,对 exonerate 结果进行分析,得到 hint 信息。

将 protein 序列比对到屏蔽了重复序列和转录组序列匹配区域的基因组。这样得到的 hints 主要位于转录组中未表达的基因处,以免和转录组的 hint 有冲突,或得到连接 2 个相邻基因的错误 hint。
$ makeblastdb -in genome.estMasked.fa -dbtype nucl -title genome -parse_seqids -out genome
$ blast.pl tblastn genome proteins.fasta 1e-8 24 blast 5

对 blast 的结果进行分析,挑选相似性高的 pritein 序列用于 exonerate 分析
$ perl tblastn_xml_2_exonerate_commands.pl --exonerate_percent 50 --exonerate_maxintron 20000 --cpu 24 blast.xml proteins.fasta genome.estMasked.fa

将 exonerate 结果转换为 hint 文件
$ /opt/biosoft/augustus-3.0.3/scripts/exonerate2hints.pl --in=exonerate.out --maxintronlen=20000 --source=P --minintron=31 --out=hints.protein.gff

5. Creating hints from RepeatMasker output

在转座子重复区,不倾向于存在编码区,该区域可以作为 nonexonpart hint。将 RepeatMakser 结果转换为此类 hint,有力于 exon 的预测。这样比使用屏蔽了重复序列的基因组进行 Augustus gene prediction 要好。

不对 Simple_repeat, Low_complexity 和 Unknown 这 3 类进行屏蔽,因为这些区域存在 CDS 的可能性相对较大。
$ grep -v -P "position|begin|Unknown|Simple_repeat|Low_complexity" genome.repeat.out | perl -pe 's/^\s*$//' | perl -ne 'chomp; s/^\s+//; @t = split(/\s+/); print $t[4]."\t"."repmask\tnonexonpart\t".$t[5]."\t".$t[6]."\t0\