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

此部分暂时掠过

GPU blast 的安装与使用

1. 安装 Nvidia 显卡驱动

安装 Nvidia 显卡驱动,得要禁用 CentOS 默认的 nouveau 驱动,然后安装 Nvidia 驱动

下载显卡驱动
# wget  http://cn.download.nvidia.com/XFree86/Linux-x86_64/340.65/NVIDIA-Linux-x86_64-340.65.run

禁用 Nvidia 显卡驱动
# echo "blacklist nouveau" >> /etc/modprobe.d/blacklist.conf

备份并重建系统引导的 image 文件
# mv /boot/initramfs-$(uname -r).img /boot/initramfs-$(uname -r).img.bak
# dracut -v /boot/initramfs-$(uname -r).img $(uname -r)

重启后在 init3 模式下安装显卡驱动
# init 3
# sh NVIDIA-Linux-x86_64-340.65.run

2. 检测显卡性能参数

推荐使用 CUDA-Z 来检测显卡性能参数。正常使用 CUDA-Z 需要依赖一些库文件和 GCC 。

$ wget http://downloads.sourceforge.net/project/cuda-z/cuda-z/0.9/CUDA-Z-0.9.231.run
$ mv CUDA-Z-0.9.231.run ~/bin/CUDA-Z

运行 CUDA-Z 前安装如下一些系统软件
$ sudo yum install libXrender*
$ sudo yum install libXrender*i686*
$ sudo yum install libXext*
$ sudo yum install libXext*i686*
$ sudo yum install libz*
$ sudo yum install libz*i686*
将 GCC 的库文件路径设置正确
$ export C_INCLUDE_PATH=/opt/gcc-4.7.2/include
$ export LD_LIBRARY_PATH=/opt/gcc-4.7.2/lib64:/opt/gcc-4.7.2/lib

运行 CUDA-Z 图形化界面,界面和 windows 下的 CPU-Z 一致。
$ CUDA-Z

显卡一些高级常识,例如:
我的笔记本电脑 GTX570M 显卡: 该 GPU 有 7个 SM,每个 SM 有 48 个 SP,每个 SP 有 32 个 warp; SP 称为流处理器,该值越大,则显卡性能越好; warp 是最小的硬件执行单位,则该显卡最大的线程数是 7*48*32=10752。
SM 之间可以理解是物理隔绝的,不同的 SM 的计算是独立的。将数据用显卡进行计算时,进行并行计算时,我个人理解成这样: 把运行的程序做成很多小份进行并行计算;所有的并行计算的总体称为 Grid,并行运算的最小单元为 thread;gpu-blast在默认设置中,将 64 个 threads 合并为 1 个 block,因此 1 个 block 的运行需要 2 个 SP 进行运算;gpu-blast在默认设置中,运行一个程序,将之分割成 521 个 blocks, 这 512 个 blocks 则需要 1024 个 SP 进行运算,可以达到最快的运算速度。所以,如果 GPU 的 SP 总数超过 1024, 则不能完全发挥显卡的性能。
对于 GTX570M 显卡,我设置 gpu-blast 的参数 -gpu_threads 32 -gpu_blocks 336 我觉得是比较好的。这样,将 blast 程序分割为 336 个 blocks 进行并行计算;将所有的 block 分配给所有的 SP ,刚好每个 SP 都运行 32 线程。但实际上,使用 gpu-blast 默认的参数也更好地发挥显卡性能,我个人实验结果是,-gpu_threads 128 -gpu_blocks 512 的结果更加不错。

3. 安装 CUDA Toolkits

CUDA 的下载网址: https://developer.nvidia.com/cuda-downloads

$ wget http://developer.download.nvidia.com/compute/cuda/6_5/rel/installers/cuda_6.5.14_linux_64.run
$ sudo init 3
# sh cuda_6.5.14_linux_64.run

4. 安装 GPU-Blast

GPU-Blast官网: http://archimedes.cheme.cmu.edu/?q=gpublast
GPU-Blast 提供的最高 blast 版本是 2.2.28 版本,需要 ncbi-blast 同版本源码文件支持。

$ mkdir /opt/biosoft/gpu_blast
$ cd /opt/biosoft/gpu_blast

从源码包安装 ncbi-blast 2.2.28 版本(此步骤可选)
$ wget ftp://ftp.ncbi.nih.gov/blast/executables/blast+/2.2.28/ncbi-blast-2.2.28+-src.tar.gz
$ tar zxf ncbi-blast-2.2.28+-src.tar.gz
$ cd ncbi-blast-2.2.28+-src/c++
$ ./configure && make -j 8 

安装 GPU-Blast 前需要将 CUDA Toolkits 的库文件与头文件加入到能识别的路径
$ echo 'export LD_LIBRARY_PATH=/usr/local/cuda/lib:/usr/local/cuda/lib64:$LD_LIBRARY_PATH' >> ~/.bashrc
$ echo 'export C_INCLUDE_PATH=/usr/local/cuda/include/:$C_INCLUDE_PATH' >> ~/.bashrc
$ source ~/.bashrc

下载并安装 GPU-Blast
$ wget http://eudoxus.cheme.cmu.edu/gpublast/gpu-blast-1.1_ncbi-blast-2.2.28.tar.gz 
$ tar zxf gpu-blast-1.1_ncbi-blast-2.2.28.tar.gz
$ perl -p -i -e 's/LATEST/2.2.28/' install
$ perl -p -i -e 's/make /make -j /' install
$ ./install
Do you want to install GPU-BLAST on an existing installation of "blastp" [yes/no]
如果输入 no 则会自动下载 ncbi-blast-2.2.28+-src.tar.gz,并进行安装,需要做上述修改,否则下载地址不正确。同时,默认下 make 编译是单线程运行,加入 -j 参数来并行运算,从而加快编译速度。
如果选择 yes,则需要输入 blastp 的路径,该路径是 blastp 的编译后所在的路径,而不是安装路径。输入 /opt/biosoft/gpu_blast/ncbi-blast-2.2.28+-src/c++/GCC447-Debug64/bin/ , 按 Enter 进行安装。

安装完毕后,将程序文件所在路径加入环境变量 PATH
$ ln -s ncbi-blast-2.2.28+-src/c++/GCC447-ReleaseMT64/bin/ bin
$ echo 'PATH=$PATH:/opt/biosoft/gpu_blast/bin' >> ~/.bashrc
$ source ~/.bashrc

安装完毕后,程序文件位于

4. 运行 GPU-Blast

以 Swissprot 数据库为例:

首先,以 fasta 文件构建 blast 数据库,需要加入 -sort_volumes 参数
$ makeblastdb -dbtype prot -in uniprot_sprot.fasta -out uniprot_sprot -sort_volumes -max_file_sz 500MB

然后,在普通 blast 数据库基础上构建 gpu 数据库。使用 -method 2 参数生成 gpu 数据库文件
$ blastp -db uniprot_sprot -query test.fa -gpu t -method 2

再使用 -gpu t 运用 gpu_blast,否则是单纯用 CPU 进行计算
$ blastp -db uniprot_sprot -query test.fa -gpu t

使用 gpu_blast 确实能加快速度,在软件的示例中,在 -num_threads 2 下,GPU 有 2.5 倍的加速。
我的使用经验: 对 40 条蛋白序列使用 blastp 比对到 Swissprot 数据库, 设置 -num_threads 8, gpu_blast 耗时 1m41.053s, cpu_blast 耗时 0m52.264s; 对 8 条蛋白质序列使用 blastp 比对到 Swissprot 数据库, 设置 -num_threads 2,gpu_blast 耗时 0m7.892s, cpu_blast 耗时 0m12.717s;对同样 8 条蛋白质序列使用 blastp 比对到 Swissprot 数据库, 设置 -num_threads 8,gpu_blast 耗时 0m5.288s, cpu_blast 耗时 0m5.836s;
个人使用感觉: 虽然有 GPU 加速,但是依然需要耗费大量 CPU,而不是只运用了 GPU 进行计算; 如果 CPU 本身运行速度较快的时候,使用 GPU 加速的比率则较低了。

对转录组数据进行质量控制分析

1. 转录组数据质量控制分析的内容

主要包含几个方面:

1. 测序数据的质量控制,推荐使用 Trimmomatic 对低质量数据进行截除;
2. 是否是链特异性测序,以及链特异性测序的种类;
3. 转录组数据对参考基因组的匹配情况,包括分析匹配到基因内或基因间区的数据量;
4. 转录组数据中 rRNA 数据含量;
5. 转录子从 5' 到 3' 端的覆盖度分析;
6. 转录组测序数据量的饱和度分析。

本文主要讲解后面 4 个内容的分析方法。

2. 使用 RNA-SeQC

下载 RNA-SeQC 及其说明文档

$ wget http://www.broadinstitute.org/cancer/cga/tools/rnaseqc/RNA-SeQC_v1.1.8.jar
$ wget http://www.broadinstitute.org/cancer/cga/sites/default/files/data/tools/rnaseqc/RNA-SeQC_Help_v1.1.2.pdf

程序运行需求的输入文件:

genome.fasta        基因组的序列文件,同时需要有相应的 .fai 和 .dict 文件。
genome.gtf          基因组注释文件,该软件仅支持 GTF 格式。标准的 gtf 文件必须有 start_codon, stop_codon 和 CDS 内容,但本软件要求 gtf 文件必须还含有 exon 内容。
accepted_hits.bam   转录组数据比对到基因组的结果,例如 tophat 的结果。此 bam 文件必须含有 RG 参数。
rRNA.fasa           rRNA 序列,用于计算 rRNA 测序数据含量。此文件不是必须的。
rRNA_intervals.list rRNA 的结构信息文件,每行以 "chromosome_id:start-end" 信息表示 rRNA 结构。此文件不是必须的。但和上一个文件必须要二选一。

准备输入文件

得到基因组的 .fai 和 .dict 文件
$ samtools faidx genome.fasta
$ java -jar /opt/biosoft/picard-tools-1.95/CreateSequenceDictionary.jar R=genome.fasta O=genome.dict

若运行 tophat 没有输入 RG 参数,则运行下面命令得到 RG 参数
$ java -jar /opt/biosoft/picard-tools-1.95/AddOrReplaceReadGroups.jar I=accepted_hits.bam O=accepted_hits.RG.bam ID=A1 LB=A1 PL=ILLUMINA SM=A1 PU=run

若仅有 GFF3 文件,没有标准的 GTF 文件,需要自己编写程序进行转换
$ gff3ToGtf.pl genome.fasa genome.gff3 > genome.gtf
$ perl -p -i -e 's/^\s*$//' genome.gtf
若不去除 genome.gtf 文件中的空行,则会出现错误提示: No attributes were found in the GTF file on line 。

运行 RNA-SeQC 命令

对一个数据运行 RNA-SeQC,并结果输出到制定文件夹
$ java -jar RNA-SeQC_v1.1.8.jar -o A1 -r genome.fasta -s "A1|accepted_hits.RG.bam|A1_note" -t genome.gtf -BWArRNA rRNA.fasta

对多个数据运行 RNA-SeQC
$ echo -e "A1\tA1.bam\tA1_note
A2\tA2.bam\tA2_note
B1\tB1.bam\tB1_note
B2\tB2.bam\tB2_note" > samples.txt
$ java -jar RNA-SeQC_v1.1.8.jar -o rna-seqc_out -r genome.fasta -s samples.txt -t genome.gtf -BWArRNA rRNA.fasta

注意,多个数据运行 RNA-SeQC,输入的 samples.txt 文件内容为:
第一行信息为: "Sample ID\tBam File\tNotes", 这行头部会被忽略,仅起注释作用。
后面多行则为相关的信息,用3列表示,用tab分割,注意Bam文件必须为绝对路径。

使用 HISAT 将转录组数据 mapping 到基因组序列

1. HISAT 简介

HISAT:hical Indexing for Spliced Alignment of Transcripts。该软件和 tophat 开发于同一个单位,它相当于 tophat 的升级版,比对速度比 tophat2 快50倍。

2. HISAT 的下载和安装

$ wget http://www.ccb.jhu.edu/software/hisat/downloads/hisat-0.1.5-beta-Linux_x86_64.zip
$ unzip hisat-0.1.5-beta-Linux_x86_64.zip -d /opt/biosoft/
$ echo 'PATH=$PATH:/opt/biosoft/hisat-0.1.5-beta/' >> ~/.bashrc
$ source ~/.bahsrc

2. HISAT 使用

2.1 构建索引文件
$ hisat-build genome.fasta genome

2.2 进行比对 产用的命令与参数示例:

$ hisat -x genome -u 1000000 -p 24 -I 0 -X 500 --rna-strandness RF -1 reads.1.fastq -2 reads.2.fastq -U single.fastq -S result.sam 

hisat 与 bowtie2 的参数基本一致,常用的参数:

-x    输入索引数据库
-1    输入 reads1
-2    输入 reads2
-U    输入单段序列
-u    仅对前多少个reads进行比对
-I    最小插入片段长度
-X    最大插入片段长度
-S    输出sam文件
--rna-strandness 链特异性测序,一般为 RF
-k    一个reads比对到多个地方,报告几个结果,默认为 5 。

Centos 安装 HP Deskjet 1010 驱动

Linux 下安装 HP Deskjet 1010 驱动需要安装 HPLIP, 可以参考: HPLIP Manual Install.

$ sudo yum -y -d 10 -e 1 install cups cups-devel gcc-c++ ghostscript libjpeg-devel glibc-headers libtool libusb-devel make python python-devel PyXML openssl-devel net-snmp-devel policycoreutils-gui PyQt4 PyQt4-devel dbus-python notify-python sane-backends sane-backends-devel sane-frontends xsane python-imaging python-imaging-devel *crypto*
$ sudo ln -s /lib/libcrypt-2.12.so /lib/libcrypt.so

$ wget http://downloads.sourceforge.net/project/hplip/hplip/3.15.2/hplip-3.15.2.tar.gz
$ tar zxf hplip-3.15.2.tar.gz
$ cd hplip-3.15.2
$ ./configure --with-hpppddir=/usr/share/cups/model/HP --libdir=/usr/lib64 --prefix=/usr --enable-qt4 --disable-libusb01_build --enable-doc-build --enable-cups-ppd-install --disable-foomatic-drv-install --disable-foomatic-ppd-install --disable-hpijs-install --disable-udev_sysfs_rules --disable-policykit --disable-cups-drv-install --enable-hpcups-install --enable-network-build --enable-dbus-build --enable-scan-build --enable-fax-build
$ make -j 4 && sudo make install

$ sudo hp-setup
设置打印机

使用 r8s 估算物种分歧时间

1. r8s 简介

r8s用于通过系统发育树估计分歧时间和分子演化速率。
软件运行需要:

一个带有枝长的有根树;
比对的序列长度。

2. r8s 下载与安装

$ wget http://loco.biosci.arizona.edu/r8s/r8s.dist.tgz
$ tar zxf r8s.dist.tgz -C /opt/biosoft/
$ mv /opt/biosoft/dist /opt/biosoft/r8s
$ echo 'PATH=$PATH:/opt/biosoft/r8s/' >> ~/.bashrx
$ source ~/.bashrc

2. r8s 的使用

软件解压缩后,其中带有一个 mannual 的 PDF 文件。

2.1 r8s 的使用方法

直接输入命令 r8s 则会进入该软件的命令行界面,推荐编写了 r8s 的脚本文件后,直接运行。运行 r8s 的方式如下:

$ r8s -b -f r8s_in.txt > r8s_out.txt

-b    运行 r8s 后推出软件的命令行界面
-f    输入的 r8s 脚本文件,该文件包含 r8s 的命令行

r8s_in.txt 的一个示例如下:

#nexus
begin trees;
tree tree_1 = [&R] ((gluc:0.46,cneg:0.54):4.3,(scer:1.07,tree:0.68):0.52);
end;

begin r8s;
blformat lengths=persite nsites=300000 ulrametric=no;
MRCA asc tree scer;
MRCA bas gluc cneg;
fixage taxon=asc age=520;
constrain taxon=bas min_age=350 max_age=410;
#divtime method=PL crossv=yes cvstart=0 cvinc=1 cvnum=18;
set smoothing=100;
divtime method=PL algorithm=TN;
showage;
describe plot=cladogram;
describe plot=chrono_description;
end;

按照上述示例中,第一部分是输入进化树数据;第二部分是运行 r8s 的命令。

2.2 r8s 命令

按照上述示例,需要依次输入上面的 r8s 命令:

blformat

length      设置树的枝长单位。若枝长单位是位点替换率,则其值为 persize,则 枝长 * 序列长度 = 替换数;若枝长单位是替换数,则该参数值为 total。默认参数是 total。
nsites      设置多序列比对的序列长度。
ultrametric 是否是超度量树,一般进化树选 no。默认参数是 no。

MRCA

该命令用来设置内部节点名。示例中设置了 tree 和 scer 的 most recent common ancestor 的节点名为 asc。

fixage

该命令用来设置 MRCA 指定的节点的分歧时间,使用该指定的分歧时间作为校正,来预测其它节点的分歧时间。
r8s 需要至少有一个内部节点进行 fixage。

constrain

该命令用来约束 MRCA 指定的节点的分歧时间,设置该节点允许的最大或最小分歧时间。

divtime

该命令用来进行分歧时间和替换速率计算。总共有 4 种计算方法(LF | NPRS | PL)和 3 种数学算法(Powell | TN | QNEWT)。 一般常用且较优,是使用 PL 和 TN。
但是使用 PL 方法,则需要设置参数 smoothing 的值。 通过设置多个 smoothing 的值来进行一些计算,选择最优的值即可。一般情况下,该值位于 1~1000 能得到一个最佳(最低)的分值。

divtime method=PL crossv=yes cvstart=0 cvinc=1 cvnum=18;
上述命令,是设置 smoothing 的值从 1, 10, 100, 1000 ... 1e17, 来计算,最后得到最佳的 smoothing 值。

若使用 fixage 对 2 个节点的分歧时间进行了固定,则可以运行命令:
divtime method=PL crossv=yes fossilfixed=yes cvstart=0 cvinc=1 cvnum=18;

若使用 fixage 对 1 个节点进行分歧时间固定,同时使用 constrain 对 2 个节点进行了约束,则可以运行命令:
divtime method=PL crossv=yes fossilconstrained=yes cvstart=0 cvinc=1 cvnum=18;

得到最优的 smoothing 值后,使用 set 进行设置,然后运行 divtime 进行分歧时间和替换速率的计算。

showage

使用该命令能得到各个节点的分歧时间和替换速率。

describe

使用该命令得到进化树的图和文字结果。 其 plot 参数如下:

cladogram    得到分支树的图,图上有各个节点的编号,和 showage 的结果结合观察。
phylogram    得到进化树的图,枝长表示替换数。
chronogram   得到超度量树的图,枝长表示时间。
ratogram     得到树图,枝长表示替换速率。

phylo_description     得到树的ASCII文字结果,枝长表示替换数。
chrono_description    得到树的ASCII文字结果,枝长表示时间。
rato_description      得到树的ASCII文字结果,枝长表示替换速率。

node_info    得到节点的信息表格

使用 CAFE 进行基因家族扩张分析

1. CAFE 简介

CAFE (Computational Analysis of gene Family Evolution) 用于进行基因家族的扩张收缩分析。

2. CAFE 下载和安装

$ wget http://heanet.dl.sourceforge.net/project/cafehahnlab/cafe.linux.x86_64
$ wget http://downloads.sourceforge.net/project/cafehahnlab/CAFEv3.1_Manual.pdf
$ wget http://downloads.sourceforge.net/project/cafehahnlab/CAFEv3.1_Manual.doc
$ mv cafe.linux.x86_64 ~/bin/cafe
$ cafe

3. CAFE 的简单使用

CAFE 需要的输入:

1. 基因家族在各个物种中的数目。该文件内容有多行,以 tab 分割,第一行内容必须如下:
Description    ID    species1    species2    ...
上面 species1 和 species2 为物种名简称。后面的行为相应的值。 Description 的值可以空着; ID 为基因家族的名称;其它的为数字,表示该基因家族在相应物种中的数目。
该文件可以由 OrthoMCL 的结果统计得到。

2. 超度量树,枝长表示分歧时间。可以使用 BEAST 或 r8s 等生成。树中的物种名必须和第一个文件中的物种名对应。

3. 输入参数  λ 的值。该值表示在物种进化过程中,每个单位时间内基因获得与丢失的概率。可以由软件自己进行计算。

直接输入命令 cafe 则进入了软件的命令行界面,也可以将命令写入到 cafe 脚本中,直接运行。一个简单的示例如下:

#!~/bin/cafe
version
date
load -i orthoMCL2cafe.tab -t 16 -l log.txt -p 0.05
tree (((chimp:6,human:6):81,(mouse:17,rat:17):70):6,dog:93)
lambda -s -t (((1,1)1,(2,2)2)2,2)
report out

以下简单介绍如上命令

version
    显示软件的版本

date
    显示当前日期

load
    -i 输入的数据文件
    -t 设置程序运行的线程数,默认为 8
    -l 设置输出的日志文件,默认标准输出
    -p 设置 p_value 的阈值,默认为 0.01

tree 输入超度量树

lambda
    -l 设置 λ 的值
    -s 设置软件自动寻找最优的 λ 值
    -t 默认下所有分支的 λ 值是相同的,若需要不同的分支有不同的 λ 值,则用该参数进行设置。该参数的值和 tree 命令中的树的内容一致,只是去除了分歧时间,并将物种名换成了表示 λ 值的编号。其中,相同的编号表示有相同的 λ 值。例如,该参数的值为 (((1,1)1,(2,2)2)2,2) ,它表示 chimp,human 和 紧邻的分支有相同的 λ 值,其它分支有另外相同的 λ 值。

report out
    设置输出文件的前缀为 out

4. CAFE 的输出结果

CAFE 的输出结果为 out.cafe,该文件内容包含如下几部分:

4.1 Tree

第 1 行为输入的树的信息。

4.2 Lambda

第 2 行为 λ 值。

4.3 节点 ID

第 3 行为节点的 ID。同样是树的文字内容,不过给每个节点进行了编号,有利于后面的数值的对应。

4.4 每个基因家族中扩张或收缩的基因数目

第 4 行给出了一系列的 节点ID 对。CAFE 对这些 ID 对进行了基因家族扩张的统计。
第 6 行的值为平均每个基因家族中扩张的基因数目,负数表示基因家族收缩。

4.5 基因家族数目

第 7 行给出发生了扩张的基因家族数目;
第 8 行给出没有发生改变的基因家族数目;
第 9 行给出发生了收缩的基因家族数目;

4.6 每个基因家族的具体扩张与收缩情况

第 10 行之后是每个基因家族具体的扩张情况。其内容分为 4 列:

第 1 列: 基因家族 ID
第 2 列: 树的信息,其中每个物种名后面多个一个下划线和数字,该数字表示基因家族的数目。特别是每个节点的基因家族数目都计算出来了,从而知道在某一个分化过程中基因家族的扩张情况。
第 3 列: 该基因家族总体上的扩张情况的 p_value,不同物种中的基因数目差异越大,则 p 值越小。
第 4 列: 计算出了每个分枝的基因家族扩张的显著性。

根据 CAFE 的结果,可以自行编写程序提取信息。