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

此部分暂时掠过