samtools 和 bamtools 安装

1. samtools 安装

SAMtools 新的 1.0 版本改进比较大。 SAMtools 包含 3 部分:

Samtools : 对 SAM/BAM/CRAM 格式文件进行读、写、索引和查看操作。
BCFtools : 对 BCF2/VCF/gVCF 格式文件进行读写操作; 对 SNP 和 short indel 进行 calling/filtering/summarising
HTSlib : 用于高通量数据读写操作的 C library

新版的 samtools 将 Samtools 和 BCFtool 分离开,成为 2 个独立的软件包。 这 2 个都运用了 HTSlib,并且其软件包中都包含了 HTSlib 。
新版不的 Samtools 的官方网址更改为了: http://www.htslib.org/

1.1 安装 Samtools

$ wget http://downloads.sourceforge.net/project/samtools/samtools/1.1/samtools-1.1.tar.bz2
$ tar jxf samtools-1.1.tar.bz2 -C /opt/biosoft/
$ cd /opt/biosoft/samtools-1.1/
$ make -j 8
$ echo 'PATH=$PATH:/opt/biosoft/samtools-1.1/' >> ~/.bashrc

1.2 安装 BCFtools

$ wget http://sourceforge.net/projects/samtools/files/samtools/1.1/bcftools-1.1.tar.bz2
$ tar jxf bcftools-1.1.tar.bz2 -C /opt/biosoft/
$ cd /opt/biosoft/bcftools-1.1
$ make -j 8
$ echo 'PATH=$PATH:/opt/biosoft/bcftools-1.1' >> ~/.bashrc

2. Bamtools

Bamtools 是 Samtools 的 c++ 版本的 Samtools, 其官网: https://github.com/pezmaster31/bamtools

$ wget https://github.com/pezmaster31/bamtools/archive/master.zip -O bamtools-master.zip
$ unzip bamtools-master.zip -d /opt/biosoft/
$ cd /opt/biosoft/bamtools-master/
$ mkdir build
$ cd build/
$ cmake ..
$ make -j 8
$ echo 'PATH=$PATH:/opt/biosoft/bamtools-master/bin/' >> ~/.bashrc 
$ source ~/.bashrc

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

使用 ssh 代理

1. ssh 代理的好处

若用户有一台服务器,该服务器网络速度快,或是可以正常访问国外网站的国外主机。则可以以此服务器为跳板,通过 ssh 代理来快速访问网络内容,或翻墙访问国外网页。

2. ssh 代理设置

此方法是在使用 linux 系统下的设置,在本地机器运行命令:

$ ssh -qTfnN -D 4241 chenlianfu@122.205.95.116

表示可以通过 IP 为 122.205.95.116 的服务器为跳板,访问网络。

3. firefox 设置

设置代理服务器为 localhost, 端口 4241, socks5 类型,即可使用代理了。

3. sshd 设置保持 ssh 持续连接

修改 /etc/ssh/sshd_config 文件内容:

ClientAliveInterval 60

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 个基因,则需要打断,认为该区域是多个基因。

Apache 和 Tomcat 整合

1. Apache 和 Tomcat 的联系

Apache 主要用于支持静态网页, 而 Tomcat 主要用于支持动态网页。 比如运行 WortPress 使用 apahce, 运行 WebApollo 这样的 java 网页则使用 Tomcat。Apache 是 web 服务器, Tomcat 是应用(java)服务器。 Tomcat 可以认为是 apache 的扩展,但是可以独立于 apache 运行。
由于学校只开放了 80 端口,而 Tomcat 运行在 8080 端口。为了能让 WebApollo 在 80 端口被访问,因此需要将 java 请求通过 apache 转发给 tomcat 处理。

2. 安装 Apache 和 Tomcat

整合 Apache 和 Tomcat 需要安装 jakarta-tomcat-connectors

使用 yum 安装 Apache
$ sudo yum install httpd httpd-devel

安装  jakarta-tomcat-connectors
$ wget http://mirror.bjtu.edu.cn/apache/tomcat/tomcat-connectors/jk/tomcat-connectors-1.2.40-src.tar.gz
$ tar zxf tomcat-connectors-1.2.40-src.tar.gz 
$ cd tomcat-connectors-1.2.40-src/native/
$ ./configure --with-apxs=/usr/sbin/apxs
$ make

安装 Tomcat,直接下载二进制包解压缩即可使用

3. 配置 JK 以整合 Apache 和 Tomcat

首先,在 /etc/http/conf.d/ 目录下建立文件 mod_jk.conf,内容如下:

# 载入 jk 模块
LoadModule jk_module modules/mod_jk.so
# 指出mod_jk模块工作所需要的工作文件workers.properties的位置
JkWorkersFile /etc/httpd/conf/workers.properties
JkMountFile   /etc/httpd/conf/uriworkermap.properties
JkLogFile  /etc/httpd/logs/mod_jk.log
JkLogLevel info
JkLogStampFormat "[%a %b %d %H:%M:%S %Y]"
JkOptions +ForwardKeySize +ForwardURICompat -ForwardDirectories
JkRequestLogFormat "%w %V %T"
# 将所有servlet 和jsp请求通过ajp13的协议送给Tomcat,让Tomcat来处理
JkMount /servlet/* worker1
JkMount /*.jsp worker1

将以下内容写入 /etc/http/conf/workers.properties :

# Defining a worker named worker1 and of type ajp13
# 创建名为 worker1 的 worker
worker.list=worker1

# Set properties for worker1
# work1 使用ajp13协议与tomcat 进程通讯
worker.worker1.type=ajp13
# 监听ajp13请求的的tomcat worker主机地址
worker.worker1.host=localhost
# 默认情况下tomcat在ajp13协议中使用的端口为8009
worker.worker1.port=8009
worker.worker1.lbfactor=50
worker.worker1.cachesize=10
worker.worker1.cache_timeout=600
worker.worker1.socket_keepalive=1
worker.worker1.socket_timeout=300

将以下内容写入 /etc/httpd/conf/uriworkermap.properties :

/*=worker1
!/*.jpg=worker1
!/*.gif=worker1
!/*.png=worker1
!/*.bmp=worker1
!/*.html=worker1
!/*.htm=worker1
!/*.swf=worker1
!/*.css= worker1
!/*.js= worker1

上述内容表示: ‘/’表示网页中的相对路径; /*=worker1 表示将‘/’目录下所有的请求都交给 worker1 (workers.properties 文件中指定的)来处理;模块会首先在‘/’目录过滤掉所有图片、flash、js文件、css文件和静态网页,将剩下的文件类型自己全部处理。

4. 配置 Apache 和 Tomcat 展示指定的目录

对 Tomcat 中的 webApp 目录进行展示。
首选修改 Apache 的 /etc/httpd/conf.d/httpd.conf 文件

Alias /tomcat "/opt/biosoft/WebApollo-2014-04-03/apache-tomcat-7.0.57/webapps"
<Directory "/opt/biosoft/WebApollo-2014-04-03/apache-tomcat-7.0.57/webapps">
    Options Includes FollowSymLinks 
    AllowOverride None
    Order deny,allow
    Allow from all
</Directory>

再修改 /opt/biosoft/WebApollo-2014-04-03/apache-tomcat-7.0.57/conf/server.xml 文件

<Host name="localhost"  appBase="webapps"
      unpackWARs="true" autoDeploy="true"
      errorReportValveClass="org.bbop.apollo.web.ErrorReportValve">
      <Context path="/tomcat" docBase="/opt/biosoft/WebApollo-2014-04-03/apache-tomcat-7.0.57/webapps/" debug="0" reloadable="true" crossContext="true"/>
      ...
</Host>

以上修改都表示将目录 /opt/biosoft/WebApollo-2014-04-03/apache-tomcat-7.0.57/webapps/ 展示到网页中的 /tomcat 目录下。

5. 使用代理进行整合

此外,还有两种非常简易的方法进行整合。
默认情况下,Tomcat 对外开启了 8080 端口,以 HTTP/1.1 协议运行; 同时,也开启了 8009 端口,以 AJP/1.3 协议运行。这两种协议可以同时生效,也可以按需要仅保留一个即可。
因此,只需要使用 Apache 的 mod_proxy.so, mod_proxy_http.so, mod_proxy_ajp.so 等模块,即可使用代理的方法来整合 Apache 和 Tomcat。这些模块默认情况下都存在并是开启着的。
因此,配置如下设置即可:

ProxyPass /tomcat http://www.chenlianfu.com:8080
或
ProxyPass /tomcat ajp://www.chenlianfu.com:8009

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/

Perl 模块使用心得

使用 Math::Random 获取随机数

当需要获取一组符合正太分布的随机数:

use Math::Random qw(random_normal);
my @aa = random_normal(10,5000,250);
my $bb = random_normal(1,5000,250);

以上结果得到@aa得到均值为5000,标准差为250的10个数据。
$bb得到一个均值为5000,标准差为250的数据。

Excel::Writer::XLSX

输出 Excel 结果文件。

Spreadsheet::Read

读取 Excel 文件。

CentOS 语音生成 TTS

TTS, 即 test to speech ,表示将文本转为为阅读的语音。

1. 中文 TTS

使用 ekho 是个不错的选择,其官网:http://www.eguidedog.net/ekho_cn.php
ekho 的下载与安装

$ wget http://downloads.sourceforge.net/e-guidedog/ekho-5.7.tar.xz
$ tar Jxf ekho-5.7.tar.xz
$ cd ekho-5.7
$ sudo yum install libsndfile*
$ ./configure && make && make install

上述安装的 ekho 紧能阅读中文

$ ekho "这次公务员涨工资,博士后也跟着涨吗?现在一个月才2000多,是2006年的标准,有点少啊。"

2. 英文 TTS

espeak 和 festival 都能进行英文 TTS

$ echo "It's such a beautiful day! Why are you in front of the computer?" | festival --tts
$ echo "It's such a beautiful day! Why are you in front of the computer?" | espeak

使用 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 可变剪接图形文件制作

此部分暂时掠过