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\t.\t.\tsrc=RM\n";' > hints.repeats.gff

6. Run Augustus predictions parallel

推荐使用 hints 进行 Augustus gene prediction,这样在有 hints 区域的基因预测会非常准确。 hints 越准确,覆盖基因组的区域越广泛,则基因预测越准确。
hint 的种类很多: 有人工进行确定的 hints; 使用 RNA-Seq 数据得到 intron 和 exonpart 类型的 hint; 使用 protein 得到 intron 和 cdspart 类型的 hint; 使用 RepeatMasker 得到 rm 类型的 hints。
对这些不同类型的 hint 要采用对应的参数来指导基因预测。 这些参数主要用于指导 augustus 对相应类型的 hints 进行得分奖罚。相关的参数文件位于 /opt/biosoft/augustus-3.0.3/config/extrinsic/ 文件夹下。

一个推荐的配置参数文件如下:

[SOURCES]
M RM P E W
[SOURCE-PARAMETERS]
[GENERAL]
      start             1             1  M    1  1e+100  RM    1      1    P    1       1    E    1       1    W 1    1
       stop             1             1  M    1  1e+100  RM    1      1    P    1       1    E    1       1    W 1    1
        tss             1             1  M    1  1e+100  RM    1      1    P    1       1    E    1       1    W 1    1
        tts             1             1  M    1  1e+100  RM    1      1    P    1       1    E    1       1    W 1    1
        ass             1             1  M    1  1e+100  RM    1      1    P    1       1    E    1       1    W 1    1
        dss             1             1  M    1  1e+100  RM    1      1    P    1       1    E    1       1    W 1    1
   exonpart             1          .992  M    1  1e+100  RM    1      1    P    1       1    E    1       1    W 1    1.005
       exon             1             1  M    1  1e+100  RM    1      1    P    1       1    E    1       1    W 1    1
 intronpart             1             1  M    1  1e+100  RM    1      1    P    1       1    E    1       1    W 1    1
     intron             1           .34  M    1  1e+100  RM    1      1    P    1    1000    E    1     1e5    W 1    1
    CDSpart             1       1 0.985  M    1  1e+100  RM    1      1    P    1     1e5    E    1       1    W 1    1
        CDS             1             1  M    1  1e+100  RM    1      1    P    1       1    E    1       1    W 1    1
    UTRpart             1       1 0.985  M    1  1e+100  RM    1      1    P    1       1    E    1       1    W 1    1
        UTR             1             1  M    1  1e+100  RM    1      1    P    1       1    E    1       1    W 1    1
     irpart             1             1  M    1  1e+100  RM    1      1    P    1       1    E    1       1    W 1    1
nonexonpart             1             1  M    1  1e+100  RM    1   1.01    P    1       1    E    1       1    W 1    1
  genicpart             1             1  M    1  1e+100  RM    1      1    P    1       1    E    1       1    W 1    1

对上面配置文件的内容简单讲解:

看此配置文件之前,要看 hint 文件内容中第 9 列有 src=P 这样的信息。 Augustus 通过此信息来确定 hint 的来源,然后根据 extrinsic 配置文件中相应的参数来处理 hints 文件中的内容。

[SOURCES] 设置 hint 的来源
M RM P E W  不同来源的 hint 类型使用空格分隔
M : 手工锚定的 hint
RM : RepeatMasker 获得的 hint
P : 来源于 protein 的 hint
E : 来源于 EST 的 hint
W : 使用 RNA-Seq 进行 wiggle track coverage 分析得到的 hint

[SOURCE-PARAMETERS] 设置 hins 来源的参数
E 1group1gene  第 1 列是 hints 的来源, 第 2 列表明对该来源的 hints 进行的处理。有 2 种处理方式:
individual_liability
    例如, 1 个 group 在 hint 文件中包含多行,即由多个 hints 组成(比如, 1 个基因有多个 intron),当其中一个 hint 不正确时,默认情况下则会对整个 group 进行忽略。而加入该参数则仅忽略错误的 hints。
1group1gene:
    对 1 个 group,试图预测 1 个基因。适合使用比较完整的 transcripts 序列做的 hint。

[GENERAL]  不同类型不同来源 hints 的参数设置
前 3 列是对不同类型的 hints 的奖罚系数的设置:
第 1 列
    hint 的类型。 当 hints 文件中没有该类型的 hint 时, 则后面不同来源的 hints 数值都使用 1 。
第 2 列
    奖励系数。 当该类型 hint 全部匹配的时候,则进行奖励。奖励系数由该值和相应 hints 来源的设置一起决定。
第 3 列
    惩罚系数。 每当有 1 个碱基不匹配,则得分乘以该系数。若有 100 个碱基不匹配,是该列值的 100 次方,因此,该值设置一般略比 1 小。该值越小,越能增加基因预测的 specficity。

后面的列,表示不同的来源的 hints 的奖励惩罚系数,每个来源的 hints 设置分 3 列:
第 1 列
    设置 hint 的来源,与 hint 文件中 src 的值对应
第 2 列
    对 hint 的得分(对应 hint 文件第 6 列)进行了分级,该值表明分成了多少级。 若该值为 1, 则表示不分级,那么当 hint 所有的碱基都匹配的时候,则进行奖励,奖励系数乘以第 3 列的值;若该值大于 1,则将 hint 文件第 6 列的分数分成了多级;若 hint 文件第 6 列没有得分,则该出需要设置为 1 。
第 3 列
    若第 2 列值为 1, 则该列只有一个值,奖励系数乘以该值; 若第 2 列不为 1, 则此列分成 2 列。例如:
    D    8     1.5  2.5  3.5  4.5  5.5  6.5  7.5  0.58  0.4  0.2  2.9  0.87  0.44 0.31  7.3
    第一列为 D, 表明 DIALIGN 类型的 hint;
    第二列为 8, 表明根据其 hint 文件的第 6 列,将奖励分成了 8 个级别;
    由于第二列不是 1 , 第三列分成了 2 列。 前一列是 7 个数值,将 hint 的打分分成了 8 个级别;后一列是这 8 个级别的奖励系数乘积。

在存在 hints.gff, extrinsic.cfg 和 genome.fasta 这 3 个文件,以及 HMM 文件 training 完毕后,即可开始进行 Augustus Training。当基因组比较大时,最好进行并行计算:

将 hints.gff 文件内容和 genome.fasta 内容进行分割。不对完整的序列进行切断。此程序将基因组序列按长度进行排序后,将序列写入到一个个分割的fasta文件中。每当写入到一个fasta文件的序列长度大于设置的值时,则将下一条序列下如下一个fasta文件。同时,也将相应的 hints 信息写入对应的 hint 文件。
$ perl split_hints_and_scaffolds_for_augustus.pl --minsize 1000000 --output split genome.fasta hints.gff

并行计算
$ for x in `ls split/*.fa | perl -pe 's/.*\///; s/.fa//' | sort -k 1.14n`
do
    echo "augustus --species=my_species --extrinsicCfgFile=extrinsic.cfg --alternatives-from-evidence=true --hintsfile=split/$x.hints --allow_hinted_splicesites=atac --alternatives-from-evidence=true --gff3=on --UTR=on split/$x.fa > split/$x.out" >> command_augustus.list
done
$ ParaFly -c command_augustus.list -CPU 12
 
合并结果
$ for x in `ls split/*.out | sort -k 1.20n`
do
    cat $x >> aug.out
done

$ join_aug_pred.pl aug.out > aug.gff3

Web Apollo 使用心得

1. Maker 的基因预测结果综合了 Augustus, GeneMarkES 和 SNAP 的 ab initio 基因预测结果和 EST 与 Protein 的比对结果。所以 Maker 的基因预测结果比较准确。 Augustus 的基因预测结果提供了 UTR 预测结果,同时由于输入了 intron 信息,预测的结果对 intron 的预测很准。GeneMarkES 的预测结果对真菌是很准的,比 SNAP 要好。所以,进行手工修正的时候,为了获得最准确的基因结构,若基因有表达,特别有intron证据时候,一般使用 Augustus 的预测结果;若无表达的基因,很多时候选择 Maker 或 GeneMarkES 的预测结果。
2. intron 的长度一般在 1000bp 以内, 30bp 以上。 如果 intron 长度特别长,则要看有少 RNA-Seq 数据的比对结果,有多少 reads 能跨过该 intron 。 没有该类证据,则认为该 intron 是错误的。
3. 将 RNA-Seq 数据比对到基因组,转换得到 BigWig track 的显示中,表达量集中且高的区域一般位于 CDS 区域。若预测的基因在这些区域是 UTR,则需要注意。
4. 在表达量比较低的区域,多个软件的基因预测结果比较混乱。在这些区域选择基因的排序原则是: a) intron 长度太大太小的忽略; b) 按软件选择的顺序是 Augustus > Maker > GeneMarkES > SNAP; c) 若有多个软件给的结果一致,则选择一致的结果。
5. SNAP的基因预测结果中允许有很短的 intron,个人觉得是不对的; SNAP 预测的基因长度比其它软件短,容易打断基因。
6. Web Apollo 在非 GTAG intron模式处有感叹号码显示。 Maker,SNAP 和 GeneMarkES 等基因预测软件一般仅认可 GTAG 的剪接模式。而实际上存在很多 GCAG 和 ATAC 模式的剪接模式。
7. 一般情况下 introns 在基因中的位置分布比较均匀,长度也比较均匀。当 intron 在基因中数目较多同时碱基比例较大时候, 基因结构像一条鱼骨架一样。若 introns 的长度不均匀,分布不均匀,则基因结构极可能错误。
8. Augustus 预测的结果和 GeneMarkES 比较接近; Maker 综合多个基因预测结果的时候,时常选择 SNAP 的预测结果。
9. 若在尾部的属于 UTR 的 exon 非常短,则去掉该 exon, 甚至去掉该端的 UTR 区域。
10. 有时候选择 PASA 的基因预测结果,需要选择 Set longest ORF。 而有时候选择 Augustus 的基因预测结果,添加 UTR 区域后, 由于 Web Apollo 自动识别最长的 ORF,从而使 CDS 区发生了改变,需要放大到碱基水平,重新设定起始和终止密码子(而不是选择最长的 ORF)。
11. 当基因组一个区域,有的软件预测是 1 个基因,有的预测为多个。 若预测的 1 个基因中有个过长的 intron,并且仅仅只有少数预测结果是 1 个基因,则需要打断,认为该区域是多个基因。

Web Apollo 的安装

1. 简介

Web Apollo 属于 JBrowse 的一个插件,可用于基因组注释的手工修正。由于通过在网页上进行手工修正,可以由多个不同地理位置的人员进行修正。
Web Apollo 比较难以安装,新版安装说明: http://webapollo.readthedocs.org/en/latest/; 2014-04-03 版本的安装说明: http://www.gmod.org/wiki/WebApollo_Installation。推荐参考后者。

2. Web Apollo 的安装步骤

2.1 安装 Perl 模块,PostGreSQL 数据库

安装 Web Apollo 之前,需要安装一些 perl 模块、 PostGreSQL 数据库和 Tomcat 7。

安装 Perl 模块比较麻烦耗时,有些 Perl 模块不容易安装上去,需要耐心。在其它安装过程中,可能提示缺少一些 Perl 模块,安装相应的 Perl 模块即可
$ sudo cpan -i BioPerl JSON JSON::XS PerlIO::gzip Heap::Simple Heap::Simple::XS Devel::Size Hash::Merge Bio::GFF3::LowLevel::Parser  Digest::Crc32  Cache::Ref::FIFO

安装 PostGreSQL
$ sudo yum install postgresql postgresql-devel

启动 PostGreSQL, 第一次启动则会生成初始化的文件。
$ sudo /etc/init.d/postgresql start

修改 PostGreSQL 配置文件  /var/lib/pgsql/data/pg_hba.conf,在尾部加入一行内容,表示用户 chenlianfu 通过密码连接 postgres 数据库。
$ sudo echo "local   all         chenlianfu                        md5" >> /var/lib/pgsql/data/pg_hba.conf

2.2 下载并解压缩 Web Apollo 和 Tomcat 7

推荐使用 2014-04-03 版本的 Web Apollo。
使用 CentOS yum 安装的 Tomcat 7 貌似不可用,因此推荐直接下载二进制版本的 Tomcat 7。

下载并解压缩 Web Apollo 2014-04-03
$ wget http://icebox.lbl.gov/webapollo/releases/previous_releases/WebApollo-2014-04-03.tgz
$ tar zxf WebApollo-2014-04-03.tgz -C /opt/biosoft/
$ cd /opt/biosoft/WebApollo-2014-04-03

下载并解压缩示例数据
$ wget http://icebox.lbl.gov/webapollo/data/pyu_data.tgz
$ tar zxf pyu_data.tgz

下载 Tomcat 7 的二进制版本
$ wget http://apache.fayea.com/tomcat/tomcat-7/v7.0.57/bin/apache-tomcat-7.0.57.tar.gz
$ tar zxf apache-tomcat-7.0.57.tar.gz

修改 Tomcat 7 配置文件
$ cd apache-tomcat-7.0.57
修改 Tomcat 7 的报错设置
$ perl -p -i -e 's/(autoDeploy=.*)>/$1\n            errorReportValveClass="org.bbop.apollo.web.ErrorReportValve">/' conf/server.xml
修改 Tomcat 7 的内存设置,推荐设置 heap size 至少 1G, permgen size 至少 256M,基因组越大,设置越大。
$ perl -p -i -e 's/cygwin=false/CATALINA_OPTS="-Xms512m -Xmx1g -XX:+CMSClassUnloadingEnabled -XX:+CMSPermGenSweepingEnabled -XX:+UseConcMarkSweepGC -XX:MaxPermSize=256m"\ncygwin=false/' bin/catalina.sh

2.3 部署 Web Apollo

$ cd apache-tomcat-7.0.57/webapps/
$ mkdir webApollo_Pythium_ultimum
$ cd webApollo_Pythium_ultimum
$ jar -xvf /opt/biosoft/WebApollo-2014-04-03/war/WebApollo.war
$ chmod 755 jbrowse/bin/*

将 JBrowse 数据文件夹链接过来
$ ln -s /opt/biosoft/WebApollo-2014-04-03/data_Pythium_ultimum data

其实,到这一步,Web Apollo 算是安装完毕了,只需要通过 tomcat 访问其 webapps 文件夹下的内容即可。但是现在访问是没用的,还需要设置一些配置文件,告诉 tomcat 使用相应的用户来访问 PostGres 数据库内容和 JBrowse 数据文件。而这些用户与权限设置,以及 JBrowse 数据文件建立是难点。在如上两点准备完毕后,即可修改 webapp 的配置文件来连通相应的数据信息,从而得到网页展示结果。

2.4 设置 PostGres 和 Web Apollo 的用户,数据库和权限

Web Apollo 网页中的用户名和密码,以及用户对某条序列的修改权限存储于 PostGres 数据库中,因此需要建立 PostGres 的用户和数据库。

创建 postgres 用户 chenlianfu ,密码 123456。 和上面配置文件相对应。该用户必须是 Linux 系统用户
$ sudo su  postgres
$ createuser -P chenlianfu
Enter password for new role: 
Enter it again: 
Shall the new role be a superuser? (y/n) n
Shall the new role be allowed to create databases? (y/n) y
Shall the new role be allowed to create more new roles? (y/n) n

创建 PostGres 数据库
$ createdb -U chenlianfu apollo_Pythium_ultimum

建立数据库的表
$ cd /opt/biosoft/WebApollo-2014-04-03/
$ cd tools/user/
$ psql -U chenlianfu apollo_Pythium_ultimum < user_database_postgresql.sql

创建 Web Apollo 用户 apollo,密码 123456
./add_user.pl -D apollo_Pythium_ultimum -U chenlianfu -P 123456 -u apollo -p 123456

提取基因组序列名,将序列名导入到 postgres 数据库,并使 apollo 用户具有访问这些序列的权限
$ cd ../../
$ ./tools/user/extract_seqids_from_fasta.pl -p Annotations- -i pyu_data/scf1117875582023.fa -o seqids.txt
$ ./tools/user/add_tracks.pl -D apollo_Pythium_ultimum -U chenlianfu -P 123456 -t seqids.txt
$ ./tools/user/set_track_permissions.pl -D apollo_Pythium_ultimum -U chenlianfu -P 123456 -u apollo -t seqids.txt -a

2.5 Jbrowse 数据文件制作

$ cd /opt/biosoft/WebApollo-2014-04-03/

创建 2 个文件夹,前者用于存放 Web Apollo 修改的数据与历史痕迹,后者用于存放 JBrowse 的数据
$ mkdir annotations_Pythium_ultimum data_Pythium_ultimum

对 fasta 文件进行 JBrowse Track 设置
$ ./apache-tomcat-7.0.57/webapps/webApollo_Pythium_ultimum/jbrowse/bin/prepare-refseqs.pl \
    --fasta pyu_data/scf1117875582023.fa
添加 webApollo 插件
$ ./apache-tomcat-7.0.57/webapps/webApollo_Pythium_ultimum/jbrowse/bin/add-webapollo-plugin.pl \
    -i data/trackList.json

先将示例文件中的 maker output 进行分割,成为不同 source 的 GFF3 文件
$ mkdir pyu_data/split_gff
$ ./tools/data/split_gff_by_source.pl -i pyu_data/scf1117875582023.gff -d pyu_data/split_gff/

对 GFF3 文件进行 JBrowse Track 设置
$ ./apache-tomcat-7.0.57/webapps/webApollo_Pythium_ultimum/jbrowse/bin/flatfile-to-json.pl \
    --gff pyu_data/split_gff/maker.gff --arrowheadClass trellis-arrowhead --getSubfeatures \
    --subfeatureClasses '{"wholeCDS": null, "CDS":"brightgreen-80pct", "UTR": "darkgreen-60pct", "exon":"container-100pct"}' \
    --className container-16px --type mRNA --trackLabel maker
$ for i in `ls pyu_data/split_gff/ | grep -v maker`
do
    i=${i##*/}
    i=${i/.gff/}
    ./apache-tomcat-7.0.57/webapps/webApollo_Pythium_ultimum/jbrowse/bin/flatfile-to-json.pl \
    --gff pyu_data/split_gff/$i.gff --arrowheadClass webapollo-arrowhead --getSubfeatures \
    --subfeatureClasses '{"match_part": "darkblue-80pct"}' --className container-10px --trackLabel $i
done

JBrowse Track 建立完毕,再创建索引
$ ./apache-tomcat-7.0.57/webapps/webApollo_Pythium_ultimum/jbrowse/bin/generate-names.pl

BAM 数据设置
$ mkdir data/bam
$ cp pyu_data/*.bam* data/bam/
$ ./apache-tomcat-7.0.57/webapps/webApollo_Pythium_ultimum/jbrowse/bin/add-bam-track.pl \
    --bam_url bam/simulated-sorted.bam --label simulated_bam --key "simulated BAM"

BigWig 数据设置
$ mkdir data/bigwig
$ cp pyu_data/*.bw data/bigwig/
$ ./apache-tomcat-7.0.57/webapps/webApollo_Pythium_ultimum/jbrowse/bin/add-bw-track.pl \
    --bw_url bigwig/simulated-sorted.coverage.bw --label simulated_bw --key "simulated BigWig"

rm data_Pythium_ultimu
mv data data_Pythium_ultimu

2.6 Web Apollo 配置

配置 Web Apollo 在 tomcat 中的 webapp 设置,从而展示 JBrowse 数据。同时,连接 PostGres 数据库从而设置得到 Web Apollo 用户和权限信息。
配置文件存放路径: /opt/biosoft/WebApollo-2014-04-03/apache-tomcat-7.0.57/webapps/webApollo_Pythium_ultimum/config

主要配置文件是 config.xml,修改内容:

设置 Web Apollo 的修改结果存放路径
<datastore_directory>/opt/biosoft/WebApollo-2014-04-03/annotations_Pythium_ultimum/</datastore_directory>

设置 intron 最小长度
<default_minimum_intron_size>40</default_minimum_intron_size>

设置 PostGres 数据库登录参数
<database>
    <driver>org.postgresql.Driver</driver>
    <url>jdbc:postgresql://localhost/apollo_Pythium_ultimum</url>
    <username>chenlianfu</username>
    <password>123456</password>
</database>

设置 JBrowse 的 refSeqs.json 文件路径
<refseqs>/opt/biosoft/WebApollo-2014-04-03/apache-tomcat-7.0.57/webapps/webApollo_Pythium_ultimum/jbrowse/data/seq/refSeqs.json</refseqs>

设置物种名,必须为 2 个单词
<organism>Pythium ultimum</organism>

设置基因组序列类型
<sequence_type>sequence:scaffold</sequence_type>

设置自定义的 feature_type, 去掉该段注释,使之生效
<status>
    <status_flag>Approved</status_flag>
    <status_flag>Needs review</status_flag>
</status>

设置 GFF3 文件的 data_adapter
<data_adapter>
    <key>GFF3</key>
    <class>org.bbop.apollo.web.dataadapter.gff3.Gff3DataAdapter</class>
    <permission>publish</permission>
    <config>/config/gff3_config.xml</config>
    <options>output=file&amp;format=gzip</options>
</data_adapter>

在主配置文件中,包含了其它配置文件的路径。对 blat_config.xml 进行修改,是 Web Apollo 的 blat 工具生效。需要 blat 命令和数据库文件。
使用 blat 软件中的 faToTwoBit 命令生成 blat 的数据库文件

$ cd /opt/biosoft/WebApollo-2014-04-03
$ faToTwoBit pyu_data/scf1117875582023.fa pyu_data/scf1117875582023.2bi
t

修改 blat_config.xml 的文件内容

设置 blat 命令路径
<blat_bin>/opt/biosoft/blat/bin/blat</blat_bin>

设置 blat 的输出文件路径
<tmp_dir>/opt/biosoft/WebApollo-2014-04-03/pyu_data/</tmp_dir>

设置数据库文件
<database>/opt/biosoft/WebApollo-2014-04-03/pyu_data/scf1117875582023.2bit</database>

设置 blat 参数
<blat_options>-minScore=100 -minIdentity=60</blat_options>

其它配置文件 canned_comments.xml, gff3_config.xml, hibernate.xml, chado_config.xml 等也要经过一些修改。

2.7 运行 Web Apollo

$ cd /opt/biosoft/WebApollo-2014-04-03
$ ./apache-tomcat-7.0.57/bin/startup.sh

开启 tomcat 服务后,即可访问 web Apollo 了。 由于将 web Apollo 安装在了服务器上,访问服务器 IP 8080端口: www.chenlianfu.com:8080/

使用 SpliceGrapher 进行可变剪接分析

1. SpliceGrapher 简介

SpliceGrapher 主要用于使用 RNA-Seq 数据对已有的基因模型进行可变剪接分析,并能给出图形结果。此外, SpliceGrapher 也能接受 EST 数据作为输入。由于使用 Sam 文件作为输入,其可变剪接分析结果比较全面准确。
SpliceGraper 的使用说明非常详细,具体请见:http://splicegrapher.sourceforge.net/userguide.html

2. SpliceGrapher 下载与安装

安装 SpliceGrapher

$ wget http://cznic.dl.sourceforge.net/project/splicegrapher/SpliceGrapher-0.2.4.tgz
$ tar zxf SpliceGrapher-0.2.4.tgz -C /opt/biosoft/
$ cd /opt/biosoft/SpliceGrapher-0.2.4
$ python setup.py build
$ sudo python setup.py install

检测 SpliceGrapher 是否安装成功以及能否正常运行
$ python
>>> import SpliceGrapher
>>> SpliceGrapher.__version__
'0.2.4'
$ cd examples
$ sh run_tests.sh

此外,SpliceGrapher有较多的步骤,根据需要来决定是否运行相应的步骤。这些步骤的正常运行需要安装一些其他软件。
创建剪接位点(splice sites)模型文件时,需要安装 PyML 0.7.9或以上版本。

$ wget http://downloads.sourceforge.net/project/pyml/PyML-0.7.13.3.tar.gz
$ tar zxf PyML-0.7.13.3.tar.gz
$ cd PyML-0.7.13.3
$ python setup.py build
$ sudo python setup.py install

如果需要使用 Bam 文件,则需要安装 Pysam 0.5或以上版本。由于 Pysam 将代码放置于 google 上,因此国内无法下载。Pysam 只是 python 版本的 samtools,因此,自行使用 samtools 将 Bam 文件转换成 Sam 文件即可。

当使用 PSGInfer pipeline 进行可变剪接转录子预测的时候,需要安装 PSGInfer 1.1.3或以上版本

$ wget deweylab.biostat.wisc.edu/psginfer/psginfer-1.1.3.tar.gz
$ tar zxf psginfer-1.1.3.tar.gz -C /opt/biosoft
$ echo 'PATH=$PATH:/opt/biosoft/psginfer-1.1.3' >> ~/.bashrc

当使用 IsoLasso Pipeline 进行可变剪接转录子预测的时候,需要 UCSC 的 gtfToGenePred genePredToBed 程序,以及 IsoLasso 2.6.1或以上版本。

$ wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/gtfToGenePred
$ wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/genePredToBed

安装 Isolasso 前需要依赖 GSL(http://www.gnu.org/s/gsl/) and CGAL(http://www.cgal.org/) library支持
使用 yum 安装 gsl
$ sudo yum install gsl*
在 Isolasso 中带有 CGAL,但是需要 libboost_thread.so.1.42.0 支持,需要安装 boost C++ library
$ wget http://sourceforge.net/projects/boost/files/latest/download?source=files
$ tar zxf boost_1_57_0.tar.gz
$ cd boost_1_57_0
$ ./bootstrap.sh 
$ ./b2 -j 8
$ sudo ./b2 install
$ sudo ln -s /usr/local/lib/libboost_thread.so.1.57.0 /usr/local/lib/libboost_thread.so.1.42.0
安装 Isolasso
$ wget http://alumni.cs.ucr.edu/~liw/isolasso-2.6.1.tar.gz
$ tar zxf isolasso-2.6.1.tar.gz -C /opt/biosoft
$ cd /opt/biosoft/isolasso/src
$ export CXXFLAGS="$CXXFLAGS -I/opt/biosoft/isolasso/src/isolassocpp/CGAL/include -L/opt/biosoft/isolasso/src/isolassocpp/CGAL/lib"
$ export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/biosoft/isolasso/src/isolassocpp/CGAL/lib/
$ make
$ echo 'PATH=$PATH:/opt/biosoft/isolasso/bin/' >> ~/.bashrc
程序运行时需要 CGAL library 支持,将之copy到系统能识别的路径中。
$ sudo cp isolassocpp/CGAL/lib/libCGAL.so* /usr/local/lib/
$ sudo cp isolassocpp/CGAL/include/* /usr/local/include/ -r

3. SpliceGrapher 的使用

SpliceGrapher 是的使用其实是运行一些 python 程序,这些 python 程序的输入文件大都包含两个文件: 基因组fasta文件(genome.fasta)和基因模型文件(geneModels.gff3)。
此处讲解 SpliceGrapher 的使用,使用了 RNA-Seq 数据用于可变剪接预测,在当前工作目录下有如下文件作为输入:

$ ls
genome.fasta geneModels.gff3 reads.1.fastq reads.2.fastq tophat.bam

有些转录组数据量较大,推荐使用 Trinity 进行标准化后的数据。

由于多个 python 程序都需要 genome.fasta 和 geneModels.gff3 作为输入,因此,软件接受 Linux 环境变量指定这两个文件的路径,或单独使用 -f 和 -m 参数接受基因组文件和基因模型文件的输入。
使用如下命令来进行环境变量设置,有利于后续程序的简单运行。

$ export SG_FASTA_REF=$PWD/genome.fasta
$ export SG_GENE_MODEL=$PWD/geneModels.gff3

3.1 创建剪接位点模型文件 (Create Splice Site Classifiers)

以下命令是运行 build_classifiers.py 来创建剪接位点模型,用于鉴定 GT-AG/GC-AG/AT-AC 位点。一般物种的剪接位点主要是 GT-AG/GC-AG,使用参数 -d gt,gc 即可。

$ mkdir 1.create_classifiers
$ cd 1.create_classifiers
$ build_classifiers.py -d gt,gc,at -a ag,ac -l create_classifiers.log

程序运行完毕,查看模型的 ROC scores
$ grep roc *.cfg
得到的 donor 的分数一般是 0.9 以上,比 acceptor 的分数要高。

build_classifiers.py 程序首先调用 generate_splice_site_data.py 生成用于 training 的序列,再调用 select_model_parameters.py 生成模型文件。
本步骤的主要结果文件是 classifiers.zip 。
运行 build_classifiers.py 命令一步到位,生成模型文件 classifiers.zip。若需要更加准确的模型文件,则参考软件的使用说明分步选择更好的参数运行。

3.2 过滤 RNA-Seq reads 的比对结果 (Filter Alignment)

使用上一步的模型文件,对 RNA-Seq 的比对结果进行过滤。

$ mkdir ../2.filter_alignments
$ cd ../2.filter_alignments
$ samtools view -h ../tophat.bam > tophat.sam
$ ln -s ../1.create_classifiers/classifiers.zip ./
$ sam_filter.py tophat.sam classifiers.zip -o filtered.sam -v

本步骤得到过滤后的比对结果文件 filtered.sam。

3.3 预测可变剪接 Graph (Predict splice graphs)

$ mkdir ../3.predict_graphs
$ cd ../3.predict_graphs
$ predict_graphs.py ../2.filter_alignments/filtered.sam -v

本步骤得到可变剪接 Graph。结果表现方式为: 默认再当前目录下生成许多的文件夹,每个文件夹以 chromosome 名称的小写字母作为文件名;每个文件夹下有很多 GFF 文件,每个 GFF 文件以 gene model 的大写字母作为文件名前缀; 每个 GFF 的内容即为 Graph 的文本形式,注意其格式不是真正的 GFF 文件格式。

由于 predict_graphs.py 为单线程运行。如果基因组较大,基因模型较多,或 Sam 文件较大,则运行时间极长。可以考虑将文件以染色体为单元进行分割,然后并行运算:

$ sam_collate.py ../2.filter_alignments/filtered.sam
$ for i in `ls *.sam`
do
    echo "predict_graphs.py $i" >> command.predict_graphs.list
done
$ ParaFly -c command.predict_graphs.list -CPU 24

此外, 快速运行的另外一种方法是:将 sam 文件转变成 depths 文件,再用于可变剪接预测。

$ sam_to_depths.py ../2.filter_alignments/filtered.sam -o filtered.depths
$ predict_graphs.py filtered.depths -v

3.3 重新比对 RNA-Seq reads 来对可变剪接 Graphs 进行更新 (Realignment RNA-seq reads)

SpliceGraper 进行可变剪接预测可能得到一些新的 exons, 其中可能包含一些错误的 exon。因此,根据上一步的结果提取 putative transcripts,然后使用 BWA 将 RNA-Seq reads 比对上去,验证 realigned reads 能否完全覆盖新的 exon 及其边界,从而对可变剪接的预测结果进行了更新。

$ mkdir ../4.realignment_pipeline
$ cd ../4.realignment_pipeline
$ realignment_pipeline.py ../3.predict_graphs/ -1 ../reads.1.fatsq -2 ../reads.2.fastq -v

本步骤得到的主要结果为更新后的可变剪接 Graphs,其结果的呈现方式和上一步骤一样。

3.4 可变剪接转录子预测 (Transcript Prediction Pipelines)

SpliceGrapher 提供了 2 种方法进行可变剪接转录子的预测: PSGInfer Pipeline 或 IsoLasso Pipeline。前者的输入文件是 fastq 文件,后者的输入文件是 Sam 文件。这两种方法都需要上一步骤的结果文件夹作为输入。

3.4.1 PSGInfer Pipeline

使用 PSGInfer pipeline 进行可变剪接转录子预测的时候,需要安装 PSGInfer 1.1.3或以上版本。

$ mkdir ../5.Transcript_Prediction
$ psginfer_pipeline.py ../4.realignment_pipeline ../reads.1.fatsq ../reads.2.fastq -d psginfer -l psginfer.log

3.4.2 Isoasso Pipeline

使用 IsoLasso Pipeline 进行可变剪接转录子预测的时候,需要 UCSC 的 gtfToGenePred genePredToBed 程序,以及 IsoLasso 2.6.1或以上版本。

Isolasso 提供了两种算法,使用 CEM 算法则开启 -C 参数。
$ isolasso_pipeline.py ../4.realignment_pipeline ../2.filter_alignments/filtered.sam -t 1.00 -d isolasso_PSG -v
$ isolasso_pipeline.py ../4.realignment_pipeline ../2.filter_alignments/filtered.sam -t 1.00 -Cd isolasso_CEM -v

值得注意的是,软件在chromosome名称的大小写转换上存在bug,导致生成的 bed 格式文件的 chromosome 名称错误,从而使结果不正确。最好了解 Isolasso 软件的使用从而分步运行。
$ generate_putative_sequences.py ../4.realignment_pipeline/.lis -A -f $SG_FASTA_REF -M _forms.gtf -m _forms.map
$ gtfToGenePred _forms.gtf stdout | genePredToBed | sort -k1,1 -k2,2n > _forms.bed
修改 bed 文件中的 chromosome 名称,例如:
$ perl -p -i -e 's/Le01scaffold/LE01Scaffold/' _forms.bed
$ runlasso.py -x _forms.bed --forceref  -o isolasso_PSG ../2.filter_alignments/filtered.sam
$ isolasso_update_graphs.py ../4.realignment_pipeline/.lis _forms.map isolasso_PSG.pred.gtf -t 1.00 -d isolasso_PSG

3.4.3 GTF 转换为 GFF3

以上得到 GTF 的结果文件,一般需要将之转换为 GFF3 格式文件:

$ convert_models.py psginfer/putative_forms.gtf -gff3=spliceGrapher.gff3

3.5 可变剪接图形文件制作

此部分暂时掠过