使用EVM整合基因预测结果

1. EVM下载和安装

EVM-sourceforge-download中下载EVM,并解压缩后直接用于EVM的使用。

2. EVM的使用步骤

2.1 收集gff3格式的基因注释内容

在EVM中,基因注释分为3种,并将所有的注释信息(无论使用那种方法进行注释的)归入到3个文件中:gene_predictions.gff3、protein_alignments.gff3 和 transcript_alignments.gff3 这三个文件中。当然这三个文件的名称可以是其它的名字,但是EVM只是接受这3种文件。

$ cat genemark_hmm.gff3 snap.gff3 aug2.out.gfff3 > gene_predictions.gff3
$ cat LEdodesGGTrinity.pasa_assemblies.gff3 > transcript_alignments.gff3
$ perl -p -i -e 's/^#.*//s' gene_predictions.gff3 transcript_alignments.gff3

2.2 制作Evidence Weights File

3种注释信息有这不同的可信度。依靠自己的直觉来设置不同注释信息的权重。由于PASA是依赖于转录信息的,其权重高;而ab initio预测的则要低。

比如:我使用转录组测序的数据进行了PASA预测;通过PASA的结果,使用Augustus构建HMM模型,结合转录组数据构建hints来进行基因预测;通过PASA结果,使用SNAP来构建HMM模型,进行ab initio基因预测;使用GeneMark_ES进行自我训练建模来ab initio基因预测。

以上4中方法中,分别的权重值则是:PASA > Augustus > SNAP > GeneMark_ES。

一般的建议是: weight(pasa) >> weight (protein) >= weight(prediction)

由于使用了转录组的数据,我对以上4个预测方法的权重是:

echo "ABINITIO_PREDICTION   AUGUSTUS    6
ABINITIO_PREDICTION GeneMark.hmm    1
ABINITIO_PREDICTION SNAP    2
TRANSCRIPT  assembler-LEdodesGGTrinity  10" > weights.txt

该权重文件格式为:三列;第一列有三种,分别是ABINITIO_PREDICTION、PROTEIN 和 TRANSCRIPT;第二列是gff3文件中对应的type一列,即gff3文件的第3列;第3列则是权重大小。

3.3 制作PASA-supported terminal exons supplement

该终止外显子的信息文件,是可选的,并推荐使用。但是可能由于PASA最新版本的程序和EVM使用文档中的描述的不一致,该文件的制作,有点难度。最后制作出来了,使用该文件的时候,却是程序不能运行。需要再次摸索。就不阐述了。略过该步骤即可。

3.4 运行EVM的分块程序

使用 partition_EVM_inputs.pl 程序将输入的需要预测的scaffolds[或chromosomes、contigs]分开到单独的文件夹中。结果得到很多个文件夹,每个文件夹对应1条序列,以及相应的gff3注释信息;同时也得到了一个分块的信息文件。该程序的运行和参数为:

$ $EVMHome/EvmUtils/partition_EVM_inputs.pl --genome genome.fasta --gene_predictions gene_predictions.gff3 --transcript_alignments transcript_alignments.gff3 --segmentSize 500000 --overlapSize 10000 --partition_listing partitions_list.out

--genome                * :fasta file containing all genome sequences
--gene_predictions      * :file containing gene predictions
--protein_alignments      :file containing protein alignments
--transcript_alignments   :file containing transcript alignments
--pasaTerminalExons       :file containing terminal exons based on PASA long-orf data.
--repeats                 :file containing repeats to be masked
--segmentSize           * :length of a single sequence for running EVM
--overlapSize           * :length of sequence overlap between segmented sequences
--partition_listing     * :name of output file to be created that contains the list of partitions
-h                        : help message with more detailed information.
(*required options)

3.5 生成EVM命令

上一步生成了很多文件夹,这一步对每个文件夹生成一个EVM运行的命令,并将这些命令放入的文件commands.list中。

$ $EVMHome/EvmUtils/write_EVM_commands.pl --genome genome.fasta --gene_predictions gene_predictions.gff3 --transcript_alignments transcript_alignments.gff3 --weights `pwd`/weights.txt --output_file_name evm.out partitions_list.out > commands.list

3.6 执行commands.list中的命令

对于这些命令,可以使用计算机集群进行多线程的运算;当然也能将这些命令拆成多份,一起运算,也就是多线程的了;可以使用EVM自带的一个程序来运行,这样有好的结果展示。多线程运行能节约时间。每个命令的执行后能在对应的文件夹生成一个名为 evm.out 的结果文件(上一步的–output_file_name参数的设定)。

$EVMHome/EvmUtils/execute_EVM_commands.pl commands.list | tee run.log

3.7 和并EVM的文件

将每个文件夹中的EVM输出结果文件进行合并,得到最后的gff3文件。

$EVMHome/
EvmUtils/convert_EVM_outputs_to_GFF3.pl --partitions partitions_list.out --output_file_name evm.out --genome genome.fasta

GBrowse的安装和使用

1. GBrowse的安装

GBrowse安装说明文档:http://gmod.org/wiki/GBrowse_2.0_Install_HOWTO
GBrowse的安装很少有能顺利安装成功的。需要不断的摸索,看文档,并搜索相关错误,google看别人是怎么解决的。有管一些我安装过程遇到的困难如下:

1.1 安装 io-lib和Bio-SCF, io-lib是安装Bio-SCF所必须的。在这一步容易出问题,导致无法安装。

$ wget http://downloads.sourceforge.net/project/staden/io_lib/1.13.0/io_lib-1.13.0.tar.gz
$ make -j 8
$ sudo make install

$ wget http://search.cpan.org/CPAN/authors/id/L/LD/LDS/Bio-SCF-1.03.tar.gz
$ tar zxf Bio-SCF-1.03.tar.gz
$ cd Bio-SCF-1.03
$ perl Makefile.PL
$ make
$ sudo make install

1.2 安装Bio::Graphics

当使用CPAN安装的时候test错误,使用手工安装,不用进行test。

2. Gbrowse的配置

2.1 Gbrowse文件

安装好Gbrowse后,须知道几个主要的文件存放位置,默认如下:

GBrowse-2.54/ 解压的GBrowse安装目录,其中bin目录下有很多重要的程序,比如创建Gbrowse2的用户的程序等。

/etc/gbrowse2/ 存放Gbrowse的配置文件,有全局的配置文件 GBrowse.conf 和 自己建立的一个物种的配置文件 species.conf;

/var/www/html/gbrowse2/ 存放gbrowse2的一些网页文件,比如gbrowse2的使用教程等;

/var/www/cgi-bin/gb2/ gbrowse2的web程序可执行文件,

/var/lib/gbrowse2/ gbrowse2的数据库目录和用户目录等。需要修改数据库目录的用户拥有者,以便于导入数据。

$ sudo chown chenlianfu /var/lib/gbrowse2/databases/

2.2 /etc/gbrowse2/GBrowse.conf 几个可能需要需要定制的配置:

1. max_render_processes   = 12 设置渲染需要的最大CPU线程数
2. image widths        = 450 640 800 1024 1280 1440 设置基因组浏览器显示的宽度,可选的值,在perference项中进行设置时可选的值
   default width       = 1024  默认的值
3. show sources        = 1  默认下开启下拉菜单 数据源 ,以利于转移到其它物种的基因组浏览中。
4. #include "themes/warm_colors"  基因组浏览器的3个主题,此为默认的主题。
   # #include "themes/transparent_colors"
   # #include "themes/solid_gray_colors"
5. default source = yeast  设置基因组浏览器中默认的物种,即http://host/cgi-bin/gb2/gbrowse/默认所访问的物种。
6. [lentinula_edodes]      设置所要建立的物种的基因组浏览器的配置文件路径
   description = Lentinula edodes Genome
   path        = lentinula_edodes.conf
7. 在GBrowse.conf同目录下有个languages的文件夹,其中gbrowse2默认使用zh.pm模块,于是GBrowse的显示文字默认是繁体中文,可以使用其目录下的其它语言模块代替。

2.3 Data Source Sections的配置

对所需要浏览基因组的物种,则需要建立一个专门的配置文件,比如本文中的 /etc/gbrowse2/lentinula_edodes.conf 文件。

[GENERAL]                                             全局设置
restrict      = require user chenlianfu hzaumycology  设置该基因组浏览器的访问权限,只有chenlianfu和hzaumycology能访问。使用命令gbrowse_create_account.pl来创建gbrowse用户
description   = Lentinula edodes Genome Database      对数据库的描述
database      = gene_Prediction_EVM                   数据库的名字,这个一定要有,不然会提示错误。当然如果将只有一个数据库,并放入到GENERAL中,则不需要该项
initial landmark = scaffold_1:20000..40000            初始访问数据库时候显示的区域
default tracks   = Genes_EVM                          初始访问数据库时候显示的tracks
                   Genes_Augustus
metadata      =                                       对数据库的一些描述
        -description    Lentinula edodes Genome (strain: W1-26)
        -maintainer     Lianfu Chen 
        -created        2013-05-28
        -modified       2013-05-28
        -authority      hzaumycology
        -coordinates    http://www.hzaumycology.com/
        -coordinates_version    1
        -source         Scaffold
        -testrange      scaffold_1:103361..105454
        -species        Lentinula edodes W1-26
example       = scaffold_1                           给出的几个例子用于选择
                scaffold_1000:2164..4463              
#################################
# database definitions                               数据库设置。
#################################
[gene_Prediction_EVM:database]                       第一个数据库。该数据库稍微特殊些,在全局中使用该数据库,需要将fasta文件和scaffold信息导入到该数据库。
db_adaptor    = Bio::DB::SeqFeature::Store           数据库的读取方法
db_args       = -adaptor DBI::mysql                  使用mysql数据库
                -dsn lentinula_edodes_EVM            mysql数据库中的数据库名
                -user chenlianfu                     mysql数据库的可写用户
                -password 1234567                    用户的密码
search options = default                             该数据库中的搜索选项

[gene_Prediction_Augustus:database]                  另外的一个数据库。将不同的基因注释放入不同的数据库,然后放入不同的tracks,有利于阅读和使用基因组。
db_adaptor    = Bio::DB::SeqFeature::Store
db_args       = -adaptor DBI::mysql
                -dsn lentinula_edodes_Augustus
                -user chenlianfu
                -password 1234567
search options = default
########################
# Default glyph settings                             默认的glyph设置
########################
[TRACK DEFAULTS]
glyph         = generic           glyph的默认类型
height        = 10                glyph的高度
r       = black
font2color    = blue
label density = 25                当labels比该数目多的时候,则关闭labels的显示以节约空间
bump density  = 100               当features的数目多于该值的时候, 则不在垂直方向上显示features,它们都被限制在了一个水平线上。
link          = AUTO              点击feature的时候,链接到feature的信息文档中。
################## TRACK CONFIGURATION ####################
# the remainder of the sections configure individual tracks   设置track
###########################################################
#[Genes_EVM:overview]                                      将该track置于overview中,而不是detail中,此时,分类自动成为overview。好处是在整体上看到该track的特征,但是细节上就没法放大了。
#[Genes_EVM:region]                                        同上。
[Genes_EVM]                                                track名
database           = gene_Prediction_EVM                   track所用到的数据库
feature            = gene                                  track所用的feature
glyph              = gene                                  track的glyph
starnd_arrow       = 1                                     glyph具有方向性;有些glyph内在就具有或不具方向性,设置该值则不影响。
bgcolor            = peachpuff                             颜色
decorate_introns   = 1                                     intron显示方法
label_transcripts  = 1
draw_translation   = 1
category           = Genes                                 track所属的分类,对应着gbrowse2的“Select Tracks"的分类
label_transcripts  = 1
visible            = 1                                     初始访问数据库时候显示该tracks
key                = Genes Predictions Intergrated by EVM  track在浏览器中的名称
citation           = EVM was used to integrate the genes prediction results of Augustus, SNAP and GeneMarkES ;PASA was used to add UTR annotations and Alternatively spliced isoforms.  该track的介绍。

[Genes_Augustus]
database           = gene_Prediction_Augustus
#feature            = gene:AUGUSTUS                         feature为gff3文件的type:source。这样做的话,就可以不必建多个mysql数据库,只需要把source设置好即可。
feature            = gene
glyph              = gene
bgcolor            = peachpuff
decorate_introns   = 1
label_transcripts  = 1
draw_translation   = 1
category           = Genes
label_transcripts  = 1
key                = Genes Predicted by Augustus

3. Gbrowse的数据的导入

3.1 EVM将多个基因组预测结果进行融合后,使用PASA加上5’和3’端非翻译区后,得到最终基因组注释结果。将该文件导入到mysql数据库 gene_Prediction_EVM 中。由于该最终的注释结果文件中mRNA的Nama的值包含的字符过长,需要进行缩短,以利于阅读;该文件中的genes也需要进行排序。因此编写程序提取出适合于gbrowse导入的gff3文件。可用于导入到gbrowse2的gff3文件有一些特点:

首先,每个scaffold、chromosome或contig之前要有一行指定其feature和name。比如:
scaffold_1 . scaffold 1 322871 . . . Name=scaffold_1
这样gbrowse才能识别scaffod_1是属于scaffold类型,并有个Name是scaffold_1。才会将其在基因组浏览器中显现出来。而正常的gff3文件是没有这样一行的。

gff3文件中mRNA中的Name有些太长,在基因组浏览器中的图片中占空间太大,需要重命名得简洁些;gff3文件中intron的可以去掉;gff3文件中将feature为transcript的改变为mRNA;去掉注释行等。以上这些都会影响gborows的显示结果。

使用如下命令来将基因预测信息导入到数据库中:

$ perl parse_evm_pasa_gff3.pl LEdodesGGTrinity.gene_structures_post_PASA_updates.26576.gff3 genome.fasta        该perl程序生成适合于导入gbrowse2的gff3文件gbrowse.gff3 和 protein.fasta文件.后者为预测的蛋白组文件。
$ mysql -h localhost -u root -p
mysql > CREATE DATABASE gene_Prediction_EVM;     创建一个名为 gene_Prediction_EVM 的 mysql数据库
EOF
mysql > Bye
$ /usr/local/bin/bp_seqfeature_load.pl -c -a DBI::mysql -d gene_Prediction_EVM -u root -p password genome.fasta gbrowse.gff3  
该程序能将gff3文件或fasta文件导入到数据库。其参数:
-c 清空数据库
-a 导入的数据库类型
-d mysql数据库对应的数据库名称
-u mysql数据库用户名
-p mysql数据库密码
该程序导入的时间有点长,依据feature的数目,时间长短不一。对于1.2万个基因的注释,需要约10分钟导入完成。

3.2 将Augustus的基因预测信息导入

perl prepare_Augustus_gff3_for_gbrowse2.pl Agustus.gff3 > gbrowse2.gff3
$ mysql -h localhost -u root -p
mysql > CREATE DATABASE gene_Prediction_Augustus;     创建一个名为 gene_Prediction_Augustus 的 mysql数据库
EOF
mysql > Bye
$ /usr/local/bin/bp_seqfeature_load.pl -c -a DBI::mysql -d gene_Prediction_Augustus-u root -p password genome.fasta gbrowse2.gff3

3.3 其它SNAP和PASA的基因预测信息导入和上面2中一致。不赘述。

4. NGS数据的导入

4.1 安装Bio::DB::Sam。需要有samtools安装,并且该samtools的安装和正常安装不一样。

$ wget http://garr.dl.sourceforge.net/project/samtools/samtools/0.1.19/samtools-0.1.19.tar.bz2
$ tar jxf samtools-0.1.19.tar.bz2
$ cd samtools-0.1.19
$ perl -p -i -e 's/CFLAGS.*/CFLAGS=     -g -Wall -O2 -fPIC #-m64 #-arch ppc/' Makefile
$ make -j 8
$ make clean        如果之前已经安装过samtools的话,需要修改makefile,再重新安装
$ make install

$ wget http://search.cpan.org/CPAN/authors/id/L/LD/LDS/Bio-SamTools-1.38.tar.gz
$ tar zxf Bio-SamTools-1.38.tar.gz
$ cd Bio-SamTools-1.38
$ perl Build.pl
$ ./Build
$ sudo ./Build install

4.2 修改data source配置文件

[NGS_Genome:database]
db_adaptor     = Bio::DB::Sam        # 数据库的读取方法                                       
db_args        = -fasta genome.fasta # 基因组的fasta文件
                 -bam bowtie2.bam    # NGS reads的比对结果
search options = default

[GenomeReadCoverageXyplot]
feature        = coverage            # 基因组测序的reads的coverage
glyph          = wiggle_xyplot       # 使用峰图来显示reads的覆盖度
database       = NGS_Genome
height         = 50
fgcolor        = black
bicolor_pivot  = 20                  # 设定一个颜色变换的coverage值
pos_color      = blue                # 当coverage > 以上设置的数值,s使用蓝色
neg_color      = red                 # 当coverage < 以上设置的该数值,s使用红色
category       = GenomeReads
label          = 0                   # labels on wiggle tracks are redundant
key            = Coverage (xyplot) of Genome NGS data

[GenomeReadCoverageDensity]
feature        = coverage
glyph          = wiggle_density      # 使用密度来显示reads的覆盖度,覆盖度越高则线条颜色越深
database       = NGS_Genome
height         = 30
bgcolor        = blue
bicolor_pivot  = 5                   # 小于该值,则线条无颜色
pos_color      = blue
neg_color      = red
category       = GenomeReads
label          = 0
key            = Coverage (density plot) of Genome NGS data

[GenomeReads]                       # reads比对到基因组的图形显示
feature        = match
glyph          = segments
draw_target    = 1
show_mismatch  = 1
mismatch_color = red
database       = NGS_Genome
bgcolor        = blue
fgcolor        = white
height         = 5
label density  = 50
bump           = fast
category       = GenomeReads
key            = Reads of Genome NGS data

[GenomeReadsPair]                  # reads pair比对到基因组的图形显示
feature       = read_pair
glyph         = segments
database      = NGS_Genome
draw_target   = 1
show_mismatch = 1
bgcolor       = sub {
                my $f = shift;
                return $f->attributes('M_UNMAPPED') ? 'red' : 'green';
                }
fgcolor       = green
height        = 3
label         = sub {shift->display_name}
label density = 50
bump          = fast
connector     = dashed
balloon hover = sub {
                my $f     = shift;
                return '' unless $f->type eq 'match';
                return 'Read: '.$f->display_name.' : '.$f->flag_str;
                }
category      = GenomeReads
key           = Read Pairs of Genome NGS data

[GenomeReadsMappingQuality]        # reas比对到基因组的Mapping质量图,高质量使用深蓝色表示,低质量使用浅蓝色显示
feature        = match
glyph          = segments
draw_target    = 1
show_mismatch  = 1
mismatch_color = red
database       = NGS_Genome
bgcolor        = sub {
        my $feature = shift;
        my $blueness = 255 - int($feature->qual * 2.40);
        my $colour = chr(35) . sprintf("%X", $blueness) .
                               sprintf("%X", $blueness) . "FF";
        return $colour;
        }
fgcolor        = black
height         = 5
label density  = 50
bump           = fast
category       = GenomeReads
key            = Reads' Mapping Quality of Genome NGS data

4.3 将Bam文件放置到配置文件中对应的位置。

若在track中提示该错误:
Track rending error: No index file for bam file; try opeing file with -autoindex at /usr/local/lib64/perl5/Bio/DB/Sam.pm line 2064
则表示缺少bam文件对应的index文件。该文件以bai为后缀,使用samtools生成

$ samtools index geonme.bam

GFF3格式

GFF3的官方介绍:Generic Feature Format Version 3 (GFF3)

1. GFF3文件格式描述

GFF3格式文件为文本文件,分为9列,以TAB分开。控制符使用 RFC 3986 Percent-Encoding 编码。比如:%20 代表着ASCII的空格。

9列文件依次是:

1. seqid:参考序列的id。该id的取名不能以’>’开头,不能包含空格。

2. source :注释的来源。如果未知,则用点(.)代替。一般指明产生此gff3文件的软件或方法。

3. type :属性的类型。建议使用符合SO惯例的名称(sequence ontology,参看[[Sequence Ontology Project]]) ,如gene,repeat_region,exon,CDS等。

4. start position :属性对应片段的起点。从1开始计数。

5. end position :属性对应片段的终点。一般比起点的数值要大。

6. score :得分,对于一些可以量化的属性,可以在此设置一个数值以表示程度的不同。如果为空,用点(.)代替。

7. strand :“+”表示正链,“-”表示负链,“.”表示不需要指定正负链。

8. phase :步进。对于编码蛋白质的CDS来说,本列指定下一个密码子开始的位置。可以是0,1或2,表示到达下一个密码子需要跳过的碱基个数。
对于其它属性,则用点(.)代替。

9. attributes :属性
一个包含众多属性的列表。格式为“标签=值”(tag=value)。不同属性之间以分号相隔。可以存在空格,不过若有“,=;”则用URL转义(URL escaping rule),同时TAB也需要转换为“%09”表示。所有以大写字幕开头的标签被保留,用于大众认可的用途,而以小写字母开头的标签则根据自己安排随意应用。
常用的标签有:
ID
Feature的标识。该ID具有唯一性。
Name
Feature的展示名称。Name的值在可视化的时候得到展示。因此,Name可以根据自己展示的需要随意取值。
Alias
Feature的第2个Name。
Parent
指明feature所从属的上一级ID。用于将exons聚集成transcript,将transripts聚集成gene。
Target
指明比对的目标区域,一般用于表明序列的比对结果。格式为”target_id start end [strand]”,其中strand是可选的(“+”或”-“), target_id中如果包含空格,则要转换成’%20’。
Gap
比对结果的gap信息,和Target一起,用于表明序列的比对结果。
Note
文本描述
Is_circular
表明featrue是否为环化的。用于环状基因组序列。

同一个tag如果有多个值,则多个值之间使用逗号隔开,比如:
Parent=AF2312,AB2812,abc-3
Alias=M19211,gna-12,GAMMA-GLOBULIN
能够使用多个值的tag有:Parent, Alias, Note, Dbxref and Ontology_term。

2. GFF3文件检测

检验GFF3格式文件: GFF3 Validator

Octave的安装

GNU Octave is a high-level interpreted language, primarily intended for numerical computations. It provides capabilities for the numerical solution of linear and nonlinear problems, and for performing other numerical experiments. It also provides extensive graphics capabilities for data visualization and manipulation. Octave is normally used through its interactive command line interface, but it can also be used to write non-interactive programs. The Octave language is quite similar to Matlab so that most programs are easily portable.

CentOS安装octave,使用源码包安装会出很多问题,使用yum安装简单容易。

# vim /etc/yum.repos.d/naulinux-school.repo

[naulinux-school]
name=NauLinux School
baseurl=http://downloads.naulinux.ru/pub/NauLinux/6.2/$basearch/sites/School/RPMS/
enabled=0
gpgcheck=1
gpgkey=http://downloads.naulinux.ru/pub/NauLinux/RPM-GPG-KEY-linux-ink

# yum --enablerepo=naulinux-school install octave*
# yum install hdf5*

# ln -s /usr/lib64/mpich2/lib/libhdf5.so.6.0.4 /usr/lib64/libhdf5.so.6

SNAP的安装和使用

1. SNAP简介

SNAP是Ian Korf独自开发的软件,简单易用。

2. SNAP的安装

$ wget http://korflab.ucdavis.edu/Software/snap-2013-02-16.tar.gz
$ tar zxf snap-2013-02-16.tar.gz
$ cd snap
$ make

3. SNAP的Parameter Estimation

3.1 需要 genes that are not too related to each other。和AUGUSTUS一致,基因两两之间的identity不要超过80%。Gene structures must be in ZFF format.

ZFF格式是Ian Korf自行使用的一个格式,有长短两种格式。In the short format, there are 4 fields: Label, Begin, End, Group.if Begin > End, the feature is on the minus strand.

>sequence-1
Einit    201    325   Y73E7A.6
Eterm   2175   2319   Y73E7A.6
>sequence-2
Einit    201    462   Y73E7A.7
Exon    1803   2031   Y73E7A.7
Exon    2929   3031   Y73E7A.7
Exon    3467   3624   Y73E7A.7
Exon    4185   4406   Y73E7A.7
Eterm   5103   5280   Y73E7A.7

不管是在正链,还是负链,从上到下,起始位置是逐渐变大的;如果是在正链上,Begin > End, Einit在Eterm前面;如果是在负链上,则 Begin < End, Eterm在Einit前面。 3.2 做Parameter Estimation需要一定数目的genes,这些genes的ZFF文件和相应的genome序列文件。

注意geneome序列文件只是包含有这些基因的序列,如果含有其它的序列的话,程序运行会出问题(core dump等)。

使用gff3_to_zff.pl来将gff3文件转换成zff文件; 再使用order.pl来调整Einit, Exon, Eterm的顺序。

$ ./gff3_to_zff.pl genes.gff3 > genes.zff
$ ./order.pl genes.zff > species.ann

提取出相应的genome的序列,并对其进行排序:

$ grep '^>' species.ann | tr -d '>' > species.seqs2keep
$ ./fasta_sort.pl species.seqs2keep < geome.fasta > species.dna

3.3 检测genes中的错误和警示,然后对genes进行修正或丢弃

$ $SnapHome/fathom species.ann species.dna -gene-stats &> gene-stats.log
$ $SnapHome/fathom species.ann species.dna -validate &> validate.log

$ grep OK validate.log > species.zff2keep
$ perl -p -i -e 's/.*:\s+(\S+)\s+OK/$1/' species.zff2keep
$ ./filterGenes.pl species.zff2keep species.ann > tmp
$ mv tmp species.ann

3.4 将序列打断成一个序列一个gene的片段,在CDS两端各加1000bp长度序列,并将所有的genes转换到正义链上。

$ $SnapHome/fathom species.ann species.dna -categorize 1000
$ $SnapHome/fathom species.ann species.dna -export 1000

3.5 run the parameter estimation program

$ mkdir params; cd params
$ $SnapHome/forge ../export.ann ../export.dna 
$ cd ..

3.6 Last is to build an HMM

$ $SnapHome/hmm-assembler.pl species params > species.hmm

SNAP的Parameter Estimation能很快地完成。不像AUGUSTUS那样耗时。

4. 运行snap进行基因预测

$ $SnapHome/snap species.hmm species.geonme > species.zff
$ $SnapHome/zff2gff3.pl species.zff > species.gff3

5. 将预测的zff结果转换成gff3结果

在4中转换成的gff3文件的最后一列attributes中,只有Name标签。而使用EVM软件将基因预测结果融合的时候,需要的gff3文件中该列的标签有ID和Parent,因此需要其它的方法来将zff文件转换成gff3文件。可以使用EVM包含的一个perl程序来解决。

使用该程序需要先将EVM自带的一些perl模块export到其路径;同时需要安装cdbfasta

$EVMHome/OtherGeneFinderTrainingGuid/SNAP/SNAP_output_to_gff3.pl species.zff species.geome > species.gff3

SNAP基因ab initio的基因预测速度很快。

Augustus training

Augustus training

1. Convert GFF file to Genbank format file

$ $AugustusHome/scripts/gff2gbSmallDNA.pl PASA.gff genome.fa 1000 genes.raw.gb

将gff文件转换成genebank格式,左右侧翼各加1000bp序列。gff文件可以由PASA将RNA-Seq的转录子比对到genome得到。而PASA得到的gff文件是有5’端非翻译区注释的,这样的信息会被trainig忽略。it is sufficient to have only the coding parts of the gene structure (CDS).

当然,genebank文件也可以使用NCBI的nucleotide数据库进行检索得到。

2. remove these problematic locis from genes.raw.gb

$ $AugustusHome/bin/etainig --species=SPECIES --stopCodonExcludedFromCDS=false genes.raw.gb 2> train.err
$ cat train.err | perl -pe 's/.*in sequence (\S+): .*/$1/' > badgenes.lst
$ $AugustusHome/scripts/filterGenes.pl badgenes.lst genes.raw.gb > genes.gb

第一条命令用于输出trainig过程中的错误信息,根据错误信息找到 badgenes,然后在去掉这些badgenes,剩下的genes用于training。

值得注意的是:

1. 至少有200个gene structures用于training,才能得到不错的结果。越多的gene,则training的效果越好;当然,达到1000个genes的时候,提升的效果就很小了。

2. 当有多于1000个基因的时候,则需要注重基因的质量,而不是数量了。要保证multi-exon genes的数目要多,这样用于tain the introns。并且gene structures越精确越好。

3. gene set should be non-redundant.如果2个不同的基因序列绝大部分的amino acid sequence是一致的,则去掉其中一个。推荐的条件是:gene set里面任意两个gene在amino acid level上的identity要不高于80%。可以使用blast来解决,由于80%的阈值算是比较高的,一般也就需要去除掉20多个基因。

3. Split gene structure set into training and test set

$ $AugustusHome/scripts/randomSplit.pl genes.gb 100

将genes.gb分隔成了genes.gb.test和genes.gb.train两个文件。其中前者为genes.gb中随机取出的100个genes,后者为剩下的genes。后者将用于不停地traning。

4. CREATE A META PARAMETERS FILE FOR YOUR SPECIES

$ $AugustusHome/scripts/new_species.pl --species=lentinula_edodes

假如我们要建立香菇物种的traning参数,则上命令建立了其参数文件和文件夹,不过文件内容是初始的。

注意的是,用于training的gene的最后一个CDS的最后3个碱基若不是终止密码子,则需要手动修改Lentinula_edodes_parameters.cfg文件,将其中的stopCodonExcludedFromCDS由默认的false改为true。

5. MAKE AN INITIAL TRAINING

$ $AugustusHome/bin/etrainig --species=lentinula_edodes genes.gb.train
$ $AugustusHome/bin/augustus --species=lentinula_edodes genes.gb.test | tee firsttest.out
$ grep -A 22 Evalustion firsttest.out

使用genes.gb.train做一次trainig,然后使用genes.gb.test来检测training的精确性。分别在nucleotide,exon和gene level上检测其sensitivity和specificity。

sensitivity表示被被检测出来的百分率;specificity表示检测出来的nucleotide,exon或gene和test set中的完全一致的百分率。

6. RUN THE SCRIPT optimize_augustus.pl

$ $AugustusHome/scripts/optimize_augustus.pl --species=lentinula_edodes --cpus=8 genes.gb.train
$ $AugustusHome/bin/etrainig --species=lentinula_edodes genes.gb.train
$ $AugustusHome/bin/augustus --species=lentinula_edodes genes.gb.test

1. optimize_augustus.pl所做的事情:

默认情况下,optimize_augustus.pl将genes.gb.train中的genes随机分成8等份,然后使用其中的7个等份的genes做training,另外的1个做精确性评估。这样相互下来,共有8个方案,每个方案取1个等份用于精确性评估,另外7个用于training。

进行一次随机分配后再运行10次training和精确性评估,即为一次预测,得到一个target value。该值是 base,exon和gene level上sensitivities和specificities的权重值。

每次预测,如果得到更高的target value,则修正参数文件中的值:lentinula_edodes_parameters.cfg。

默认下参数文件中有28项参数需要按一定顺序进行优化;一般情况下每个参数最多设置5个值各进行一次预测(即对一项参数而言,这设置的5个值其中可能有1个值是用于之前的预测,故每个参数优化需要运行最多5次预测),取最大的target value对应的值为参数的值;对所有的参数进行优化一次是一轮,则5轮参数优化完毕后程序会停止运行(以1800左右个genes来进行training,则每次augustus对200多个gene进行预测需要1min,那么每个参数优化需要28*4*8*1min=896min=15h,5轮参数的优化总共需要75h,即3天),或如果在一轮参数优化中没有improvements则提前停止运行。当然,如果等不及,也能手动停止程序运行。由于optimize_augustus.pl运行时间太长,最好使用screen来运行。

如果了解了上述运行原理,则可以视情形终止其运行,或保存配置文件后接着运行。

2. 在optimize_augustus.pl完成或中断之后,需要(re)train AUGUSTUS with genes.gb.train。然后在使用genes.test.gb进行预测的精确性检测,如果gene level sensitivity低于20%,则表明training set不够大,或者质量不够好,或者物种somehow special。

7. Training AUGUSTUS UTR parameters

这部分的Training则需要5’和3’端的UTR都存在的gene structure。

Augustus使用技巧

1. Predict Genes ab initio

Ab initio prediction means that no other input is used than the target genome itself. Predict the genes in the range 7,000,001-7,500,000 of chr2R of D. melanogaster. Use the FASTA format file chr2R.fa, which includes the whole chromosome 2R.

augustus --species=fly --predictionStart=7000001 --predictionEnd=7500000 chr2R.fa > augustus.abinitio.gff   # takes ~1m

If you want the protein sequences you can retrieve them with

getAnnoFasta.pl augustus.abinitio.gff

2. Predic Genes Using Hints

augustus --species=fly --predictionStart=7000001 --predictionEnd=7500000 chr2R.fa \
  --extrinsicCfgFile=extrinsic.bug.cfg --hintsfile=hints.gff > augustus.hints.gff

3. Creating hints from ESTs or assembled RNAseq transcripts with BLAT

4. Creating hints from proteins with Exonerate

5. Run AUGUSTUS predictions parallel

6. RNAseq integration (raw reads)

6.1 Incorporating RNAseq data into AUGUSTUS predictions with BLAT (including iterative mapping)

6.2 Incorporating RNAseq data with GSNAP (including iterative mapping)

6.3 Incorporating RNAseq data with Bowtie/Tophat (including iterative mapping)

6.4 Incorporating RNAseq data with Palmapper?

7. Training AUGUSTUS

Incorporating Illumina RNAseq into AUGUSTUS with Tophat

Incorporating Illumina RNAseq into AUGUSTUS with Tophat

#!/bin/bash

###### Step 1 ---------------------

# Aligning reads with Bowtie/Tophat
bowtie2-build genome.fasta.masked genome
tophat -r 20 -mate-std-dev 40 --coverage-search --microexon-search -p 24 --library-type fr-unstranded --phred64-quals genome reads_1.fq reads_2.fq 

# Filtering raw alignments
samtools sort -n tophat_out/accepted_hits.bam tophat_out/accepted_hits.s
augustus.2.7/auxprogs/filterBam/bin/filterBam --uniq --paired --in tophat_out/accepted_hits.s.bam --out tophat_out/accepted_hits.sf.bam

samtools view -H tophat_out/accepted_hits.sf.bam > header.txt

# Creating intron hints
samtools sort tophat_out/accepted_hits.sf.bam both.ssf
augustus.2.7/auxprogs/bam2hints/bam2hints --intronsonly --in=both.ssf.bam --out=hints.gff

# Run Augustus
augustus.2.7/bin/augustus --species=SPECIES --extrinsicCfgFile=extrinsic.cfg --alternatives-from-evidence=true --hintsfile=hints.gff --allow_hinted_splicesites=atac --introns=on --genemodel=complete genome.fasta > aug1.out
######-----------------------------

###### Step 2 ---------------------

# Create an exon-exon junction database
cat aug1.out | tee aug.prelim.gff | grep -P "\tintron\t" > aug1.introns.gff
cat hints.gff aug1.introns.gff | perl -ne 'split; print "@_[0]:@_[3]-@_[4]\n";' | sort -u > introns.lst
augustus.2.7/scripts/intron2exex.pl --flank=90 --introns=introns.lst --seq=genome.fasta.masked --exex=exex.fa --map=map.psl
bowtie2-build exex.fa genome_exex1

# Aligning reads with Bowtie
bowtie2 -q --phred64 --fr -p 24 --reorder -x genome_exex1 -1 reads_1.fq -2 reads_2.fq -S bowtie.sam
samtools view -S -F 4 bowtie.sam > bowtie.F.sam
augustus.2.7/scripts/samMap.pl bowtie.F.sam map.psl > bowtie.global.sam 2> bowtie.samMap.err
# 需要将samMap.pl中的75修正为需要的reads长度。注意使用 2> 来重定向STDERR数据,以提高运行速度。

# Join data from step 1 and step 2
bamtools filter -in tophat_out/accepted_hits.bam -out tophat_out/accepted_hits.noN.bam -script augustus.2.7/auxprogs/auxBamFilters/operation_N_filter.txt
cat header.txt bowtie.global.sam > bowtie.global.h.sam
samtools view -bS -o bowtie.global.h.bam bowtie.global.h.sam
bamtools merge -in bowtie.global.h.bam -in tophat_out/accepted_hits.noN.bam -out both.bam
samtools sort -n both.bam both.s

# Filtering raw alignments
augustus.2.7/auxprogs/filterBam/bin/filterBam --uniq --paired --in both.s.bam --out both.sf.bam

# Creating intron hints
samtools sort both.sf.bam both.ssf
augustus.2.7/auxprogs/bam2hints/bam2hints --intronsonly --in=both.ssf.bam --out=hints.2.gff

#Run Augustus
augustus --species=yourSpecies --extrinsicCfgFile=extrinsic.cfg --alternatives-from-evidence=true --hintsfile=hints.2.gff --allow_hinted_splicesites=atac genome.fasta > aug2.out

Augustus的安装和使用参数

AUGUSTUS is a program that predicts genes in eukaryotic genomic sequences.

1. Augustus的安装

Augustus下载:http://bioinf.uni-greifswald.de/augustus/binaries/

$ wget http://bioinf.uni-greifswald.de/augustus/binaries/augustus.2.7.tar.gz
$ tar zxf augustus.2.7.tar.gz
$ cd augustus.2.7
$ cd src
$ make -j 8
$ export AUGUSTUS_CONFIG_PATH=$PWD/../config/ (可以加入到.bashrc中)

2. Augustus使用方法

2.1 基因预测例子

$ augustus --strand=both --genemode=partial --singlestrand=false --hintsfile=hints.gff --extrinsicCfgFile=extrinsic.cfg --protein=on --introns=on --start=on --stop=on --cds=on --codingseq=on --alternatives-from-evidence=true --gff3=on --UTR=on ----outfile=out.gff --species=human genome.fa
$ augustus --noprediction=true --species=SPECIES sequences.gb

2.2 Augustus使用参数

Usage:

augustus [parameters] --sepcies=SPECIES queryfilename

重要参数:

--strand=both, --strand=forward or --strand=backward
    report predicted genes on both strands, just the forward or 
just the backward strand.default is 'both'

--genemodel=partial, --genemodel=intronless, --genemodel=complete, 
--genemodel=atleastone or --genemodel=exactlyone
    partial : allow prediction of incomplete genes at the sequence boundaries (default)
    intronless : only predict single-exon genes like in prokaryotes and some eukaryotes
    complete : only predict complete genes
    atleastone : predict at least one complete gene
    exactlyone : predict exactly one complete gene

--singlestrand=true
    predict genes independently on each strand, allow overlapping
 genes on opposite strands. This option is turned off by default.

--hintsfile=hintsfilename
    When this option is used the prediction considering hints (ex
trinsic information) is turned on. hintsfilename contains the hints
 in gff format.

--extrinsicCfgFile=cfgfilename
    Optional. This file contains the list of used sources for the 
hints and their boni and mali. If not specified the file "extrin
sic.cfg" in the config directory $AUGUSTUS_CONFIG_PATH is used.

--maxDNAPieceSize=n
    This value specifies the maximal length of the pieces that the 
sequence is cut into for the core algorithm (Viterbi) to be run. 
Default is --maxDNAPieceSize=200000.
    AUGUSTUS tries to place the boundaries of these pieces in the 
intergenic region, which is inferred by a preliminary prediction. 
GC-content dependent parameters are chosen for each piece of DNA 
if /Constant/decomp_num_steps > 1 for that species. This is why 
this value should not be set very large, even if you have plenty 
of memory.

--protein=on/off
--introns=on/off
--start=on/off
--stop=on/off
--cds=on/off
--codingseq=on/off
    Output options. Output predicted protein sequence, introns, 
start codons, stop codons. Or use 'cds' in addition to 'initial', 
'internal', 'terminal' and 'single' exon. The CDS excludes the 
stop codon (unless stopCodonExcludedFromCDS=false) whereas the 
terminal and single exon include the stop codon.

--AUGUSTUS_CONFIG_PATH=path
    path to config directory (if not specified as environment var
iable)

--alternatives-from-evidence=true/false
    report alternative transcripts when they are suggested by hints

--alternatives-from-sampling=true/false
    report alternative transcripts generated through probabilistic 
sampling

--sample=n
--minexonintronprob=p
--minmeanexonintronprob=p
--maxtracks=n

--proteinprofile=filename
Read a protein profile from file filename. See section 7 below.

--predictionStart=A, --predictionEnd=B
    A and B define the range of the sequence for which predictions 
should be found. Quicker if you need predictions only for a small 
part.

--gff3=on/off
    output in gff3 format.

--UTR=on/off
    predict the untranslated regions in addition to the coding 
sequence. This currently works only for human, galdieria, toxopl
asma and caenorhabditis.

--outfile=filename
    print output to filename instead to standard output. This is 
useful for computing environments, e.g. parasol jobs, which do 
not allow shell redirection.

--noInFrameStop=true/false
    Don't report transcripts with in-frame stop codons. Otherwise, 
intron-spanning stop codons could occur. Default: false

--noprediction=true/false
    If true and input is in genbank format, no prediction is made. 
Useful for getting the annotated protein sequences. Augustus也可以以
genebank格式文件为输入文件,进行基因预测,并将预测结果和genebank的结果进行比较后
得出一个精确性的统计结果。
    当然,由于genebank格式文件中有些sequences没有cds的注释结果,因此可以使用该
参数进行检测,从而得到没有cds的序列号,在人为去去除这些没有cds注释的序列,再去进行
预测准确性的评估。

--contentmodels=on/off
    If 'off' the content models are disabled (all emissions unif
ormly 1/4). The content models are; coding region Markov chain 
(emiprobs), initial k-mers in coding region (Pls), intron and int
ergenic regin Markov chain. This option is intended for special 
applications that require judging gene structures from the signal 
models only, e.g. for predicting the effect of SNPs or mutations 
on splicing. For all typical gene predictions, this should be 
true. Default: on

--paramlist
    For a complete list of parameters, type "augustus --paramlist"