使用HTSeq进行有参转录组的表达量计算

1. HTSeq简介

HTSeq是使用Python编写的一支用于进行基因Count表达量分析的软件,能根据SAM/BAM比对结果文件和基因结构注释GTF文件得到基因水平的Counts表达量。HTSeq进行Counts计算的原理非常简单易懂,容易上手。

2. HTSeq安装

PYPI下载HTSeq的Python包
$ wget https://pypi.python.org/packages/46/f7/6105848893b1d280692eac4f4f3c08ed7f424cec636aeda66b50bbcf217e/HTSeq-0.7.2.tar.gz
$ tar zxf HTSeq-0.7.2.tar.gz
$ cd HTSeq-0.7.2
$ /opt/sysoft/Python-2.7.11/bin/python setup.py build
$ /opt/sysoft/Python-2.7.11/bin/python setup.py install
$ cd ../ && rm HTSeq-0.7.2 -rf

3. HTSeq使用

3.1 HTSeq的Count模式

HTSeq计算counts数有3种模式,如下图所示(ambiguous表示该read比对到多个gene上;no_feature表示read没有比对到gene上):
HTSeq Count模式

3.2 HTSeq的使用命令

HTseq安装完毕后,在Python软件的bin目录下生成htseq-count命令。
htseq-count运行简单示例:

对于非链特异性真核转录组测序数据
$ /opt/sysoft/Python-2.7.11/bin/htseq-count -f sam -r name -s no -a 10 -t exon -i gene_id -m union hisat2.sam genome.gtf > counts_out.txt
对于链特异性测序真核转录组测序数据
$ /opt/sysoft/Python-2.7.11/bin/htseq-count -f sam -r name -s reverse -a 10 -t exon -i gene_id -m union hisat2.sam genome.gtf > counts_out.txt
对于非链特异性原核生物转录组测序数据
$ /opt/sysoft/Python-2.7.11/bin/htseq-count -f sam -r name -s no -a 10 -t exon -i gene_id -m intersection-strict bowtie2.sam genome.gtf > counts_out.txt

htseq-count命令的常用参数:

-f | --format     default: sam
  设置输入文件的格式,该参数的值可以是sam或bam。
-r | --order     default: name
  设置sam或bam文件的排序方式,该参数的值可以是name或pos。前者表示按read名进行排序,后者表示按比对的参考基因组位置进行排序。若测序数据是双末端测序,当输入sam/bam文件是按pos方式排序的时候,两端reads的比对结果在sam/bam文件中一般不是紧邻的两行,程序会将reads对的第一个比对结果放入内存,直到读取到另一端read的比对结果。因此,选择pos可能会导致程序使用较多的内存,它也适合于未排序的sam/bam文件。而pos排序则表示程序认为双末端测序的reads比对结果在紧邻的两行上,也适合于单端测序的比对结果。很多其它表达量分析软件要求输入的sam/bam文件是按pos排序的,但HTSeq推荐使用name排序,且一般比对软件的默认输出结果也是按name进行排序的。
-s | --stranded     default: yes
  设置是否是链特异性测序。该参数的值可以是yes,no或reverse。no表示非链特异性测序;若是单端测序,yes表示read比对到了基因的正义链上;若是双末端测序,yes表示read1比对到了基因正义链上,read2比对到基因负义链上;reverse表示双末端测序情况下与yes值相反的结果。根据说明文件的理解,一般情况下双末端链特异性测序,该参数的值应该选择reverse(本人暂时没有测试该参数)。
-a | --a     default: 10
  忽略比对质量低于此值的比对结果。在0.5.4版本以前该参数默认值是0。
-t | --type     default: exon
  程序会对该指定的feature(gtf/gff文件第三列)进行表达量计算,而gtf/gff文件中其它的feature都会被忽略。
-i | --idattr     default: gene_id
  设置feature ID是由gtf/gff文件第9列那个标签决定的;若gtf/gff文件多行具有相同的feature ID,则它们来自同一个feature,程序会计算这些features的表达量之和赋给相应的feature ID。
-m | --mode     default: union
  设置表达量计算模式。该参数的值可以有union, intersection-strict and intersection-nonempty。这三种模式的选择请见上面对这3种模式的示意图。从图中可知,对于原核生物,推荐使用intersection-strict模式;对于真核生物,推荐使用union模式。
-o | --samout 
  输出一个sam文件,该sam文件的比对结果中多了一个XF标签,表示该read比对到了某个feature上。
-q | --quiet
  不输出程序运行的状态信息和警告信息。
-h | --help
  输出帮助信息。

3.3 HTSeq使用注意事项

HTSeq的使用有如下注意事项,否则得到的结果是错误的:

1. HTSeq是对有参考基因组的转录组测序数据进行表达量分析的,其输入文件必须有SAM和GTF文件。
2. 一般情况下HTSeq得到的Counts结果会用于下一步不同样品间的基因表达量差异分析,而不是一个样品内部基因的表达量比较。因此,HTSeq设置了-a参数的默认值10,来忽略掉比对到多个位置的reads信息,其结果有利于后续的差异分析。
3. 输入的GTF文件中不能包含可变剪接信息,否则HTSeq会认为每个可变剪接都是单独的基因,导致能比对到多个可变剪接转录本上的reads的计算结果是ambiguous,从而不能计算到基因的count中。即使设置-i参数的值为transcript_id,其结果一样是不准确的,只是得到transcripts的表达量。

3.4 HTSeq的结果

HTSeq将Count结果输出到标准输出,其结果示例如下:

gene00001	0
gene00002	9224
gene00003	880
...
gene12300	1043
gene12301	200
__no_feature	127060
__ambiguous	0
__too_low_aQual	4951
__not_aligned	206135
__alignment_not_unique	0

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

此部分暂时掠过

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

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文件必须为绝对路径。

Genome-guided Trinity for Gene Structure Annotation

使用genome来引导Trinity进行基因结构注释。

RNA-seq的一个主要用途是识别基因组的转录区,重构转录子结构,同时,鉴定转录子的可变剪切。

现在最新的基于genome的转录子预测方法是将RNA-seq的reads使用剪接比对的方法比对到基因组,然后组装比对结果从而得到转录子的结构。(eg. cufflinks, scripture)。我们将这种方法称为:align-reads then assemble-alignments

Trinity可以进行不需要参考基因组的de novo组装,见:Trinity的安装与使用;也能进行有参考基因组支持的组装:即将RNA-Seq比对到genome、RNA-Seq read的de novo组装 和 转录子比对 结合起来。

1. 步骤

1.1 align-reads

使用GSNAP来将reads比对到基因组。将基因组分成各个被reads覆盖的区。

1.2 assemble-reads

对每个区使用Trinity对相应的reads进行组装。

1.3 align-transcripts

使用PASA软件调用GMAP来将Trinity-assembled transcripts比对到genome.

1.4 assemble-transcript_alignments

使用PASA软件来组装上一步骤的比对结果,得出完整的转录子结构,同时,也能解析可变剪接的转录子结构。该步骤和上一步骤其实是在同一个PASA程序中执行得到的。

2. 需要的软件

Trinity
GSNAP & GMAP
PASA

3. 运行

Below, we describe the steps required for running the genome-guided Trinity-based transcript reconstruction pipeline. 适合于真菌物种,其基因密度较大。

3.1 Align RNA-Seq reads to the genome

$ $TRINITY_HOME/util/alignReads.pl --seqType fq --left reads.left.fq --right reads.right.fq --target genome.fasta --aligner gsnap -- -t 8
$ samtools view gsnap_out/gsnap.coordSorted.bam > gsnap.coordSorted.sam

3.2 Assemble the aligned reads using Trinity

$ % $TRINITY_HOME/util/prep_rnaseq_alignments_for_genome_assisted_assembly.pl --SS_lib_type FR --coord_sorted_SAM gsnap.coordSorted.sam -I 1000000
$ find Dir_* -name "*reads" > read_files.list
$ $TRINITY_HOME/util/GG_write_trinity_cmds.pl --reads_list_file read_files.list --paired --SS --jaccard_clip > trinity_GG.cmds
$ $TRINITY_HOME/Inchworm/bin/ParaFly -c trinity_GG.cmds -CPU 6 -failed_cmds trinity_GG.cmds.failed -v
$ find Dir_*  -name "*inity.fasta" -exec cat {} + | $TRINITY_HOME/util/inchworm_accession_incrementer.pl > Trinity_GG.fasta

3.3 Align and assemble the Trinity-reconstructed transcripts using the PASA pipeline

$ cp $PASA_HOME/pasa_conf/pasa.alignAssembly.Template.txt alignAssembly.config
$ perl -p -i -e 's/MYSQLDB=.*/MYSQLDB=sample_mysql_database/' alignAssembly.config
$ $PASA_HOME/scripts/Launch_PASA_pipeline.pl -c alignAssembly.config -C -R -g genome.fasta -t Trinity_GG.fasta --ALIGNERS blat,gmap --transcribed_is_aligned_orient --stringent_alignment_overlap 30.0

GFOLD的安装与使用

一. GFOLD简介

GFOLD是由同济大学生物信息系张勇教授课题组研发的软件。gfold主要用于RNA-seq数据的差异表达基因分析,特别适合于实验没有重复的情况。

二. GFOLD的安装

GFOLD的安装需要安装gsl-1.15

安装gsl-1.15
$ wget ftp://ftp.gnu.org/gnu/gsl/gsl-1.15.tar.gz
$ tar zxf gsl-1.15.tar.gz
$ cd gsl-1.15
$ ./configure
$ make -j 8
$ sudo make install

安装gfold
$ wget https://bitbucket.org/feeldead/gfold/get/e78560195469.zip
$ unzip e78560195469.zip
$ cd feeldead-gfold-e78560195469
$ make
$ cd ../ ; mv feeldead-gfold-e78560195469 gfold

三. GFOLD的使用

使用gfold,发现结果没有cufflinks得出的结果好。该软件相对简单快速。不推荐使用。

Cufllinks的安装与使用

一. 简介

Cufflinks下主要包含cufflinks,cuffmerge,cuffcompare和cuffdiff等几支主要的程序。主要用于基因表达量的计算和差异表达基因的寻找。

二. 安装

Cufflinks下载网页
1. 为了安装Cufflinks,必须有Boost C++ libraries。下载Boost并安装。默认安装在/usr/local。

$ tar jxvf boost_1_53_0.tar.bz2
$ cd boost_1_53_0
$ ./bootstrap.sh
$ sudo ./b2 install

2.安装SAM tools。

下载SAM tools。
$ tar jxvf samtools-0.1.18.tar.bz2
$ cd samtools-0.1.18
$ make
$ sudo su 
# mkdir /usr/local/include/bam
# cp libbam.a /usr/local/lib
# cp *.h /usr/local/include/bam/
# cp samtools /usr/bin/

3. 安装 Eigen libraries。

下载Eigen
$ tar jxvf 3.1.2.tar.bz2
$ cd eigen-eigen-5097c01bcdc4
$ sudo cp -r Eigen/ /usr/local/include/

4. 安装Cufflinks。

$ tar zxvf cufflinks-2.0.2.tar.gz
$ cd cufflinks-2.0.2
$ ./configure --prefix=/path/to/cufflinks/install --with-boost=/usr/local/ --with-eigen=/usr/local/include//Eigen/
$ make
$ make install

5. 可以直接下载Linux x86_64 binary。不需要上述繁琐步骤,解压后的程序直接可用。

三. Cufflinks的使用

1. Cufflinks简介

Cufflinks程序主要根据Tophat的比对结果,依托或不依托于参考基因组的GTF注释文件,计算出(各个gene的)isoform的FPKM值,并给出trascripts.gtf注释结果(组装出转录组)。

2. 使用方法

$ cufflinks [options]* <aligned_reads.(sam/bam)>

一个常用的例子:
$ cufflinks -p 8 -G transcript.gtf --library-type fr-unstranded -o cufflinks_output tophat_out/accepted_hits.bam

3. 普通参数

-h | --help
-o | --output-dir <sting>  default: ./
    设置输出的文件夹名称

-p | --num-threads  default: 1
    用于比对reads的CPU线程数

-G | --GTF <reference_annotation.(gtf/gff)>
    提供一个GFF文件,以此来计算isoform的表达。此时,将不会组装新的transcripts,
程序会忽略和reference transcript不兼容的比对结果

-g | --GTF-guide <reference_annotation.(gtf/gff)>
    提供GFF文件,以此来指导转录子组装(RABT assembly)。此时,输出结果会包含ref
erence transcripts和novel genes and isforms。

-M | --mask-file <mask.(gtf/gff)>
    提供GFF文件。Cufflinks将忽略比对到该GTF文件的transcripts中的reads。该
文件中常常是rRNA的注释,也可以包含线立体和其它希望忽略的transcripts的注释。将这
些不需要的RNA去除后,对计算mRNA的表达量是有利的。

-b | --frag-bias-correct <genome.fa>
    提供一个fasta文件来指导Cufflinks运行新的bias detection and correct
ion algorithm。这样能明显提高转录子丰度计算的精确性。

-u | --multi-read-correct
    让Cufflinks来做initial estimation步骤,从而更精确衡量比对到genome多个
位点的reads。

--library-type  default:fr-unstranded
    处理的reads具有链特异性。比对结果中将会有个XS标签。一般Illumina数据的lib
rary-type为 fr-unstranded。

4. 丰度评估参数

-m | --frag-len-mean default: 200
插入片段的平均长度。不过现在Cufflinks能learns插入片段的平均长度,因此不推荐自主
设置此值。

-s | --frag-len-std-dev default: 80
插入片段长度的标准差。不过现在Cufflinks能learns插入片段的平均长度,因此不推荐自
主设置此值。

-N | --upper-quartile-form
使用75%分为数的值来代替总的值(比对到单一位点的fragments的数值),作normal
ize。这样有利于在低丰度基因和转录子中寻找差异基因。

--total-hits-norm default: TRUE
Cufflinks在计算FPKM时,算入所有的fragments和比对上的reads。和下一个参数
对立。默认激活该参数。

--compatible-hits-norm 
Cufflinks在计算FPKM时,只针对和reference transcripts兼容的fragmen
ts以及比对上的reads。该参数默认不激活,只能在有 --GTF 参数下有效,并且作 RABT
或 ab initio 的时候无效。

5. 组装常用参数

-L | --label  default: CUFF
    Cufflink以GTF格式来报告转录子片段(transfrags),该参数是GTF文件的前缀

--min-frags-per-transfrag <int>  default: 10
    组装出的transfrags被支持的RNA-seq的fragments数少于该值则不被报道。

--min-intron-length <int>  default: 50
    最小的intron大小。

--overlap-radius <int>  default: 50
    Transfrags之间的距离少于该值,则将其连到一起。

6. Cufflinks输出结果

1. transcripts.gtf
该文件包含Cufflinks的组装结果isoforms。前7列为标准的GTF格式,最后一列为attributes。其每一列的意义:

列数   列的名称  例子         描述
1     序列名    chrX        染色体或contig名
2     来源      Cufflinks   产生该文件的程序名
3     类型      exon        记录的类型,一般是transcript或exon
4     起始      1           1-base的值
5     结束      1000        结束位置
6     得分      1000        
7     链        +          Cufflinks猜测isoform来自参考序列的那一条链,
一般是'+','-'或'.'
8     frame    .           Cufflinks不去预测起始或终止密码子框的位置
9     attributes  ...      详见下

每一个GTF记录包含如下attributes:

Attribute      例子       描述
gene_id        CUFF.1    Cufflinks的gene id
transcript_id  CUFF.1.1  Cufflinks的转录子 id
FPKM           101.267   isoform水平上的丰度, Fragments Per Kilobase
 of exon model per Million mapped fragments
frac           0.7647    保留着的一项,忽略即可,以后可能会取消这个
conf_lo        0.07      isoform丰度的95%置信区间的下边界,即 下边界值 =
 FPKM * ( 1.0 - conf_lo )
conf_hi        0.1102    isoform丰度的95%置信区间的上边界,即 上边界值 =
 FPKM * ( 1.0 + conf_hi )
cov            100.765   计算整个transcript上read的覆盖度
full_read_support   yes  当使用 RABT assembly 时,该选项报告所有的intr
ons和exons是否完全被reads所覆盖

2. ispforms.fpkm_tracking
isoforms(可以理解为gene的各个外显子)的fpkm计算结果
3. genes.fpkm_tracking
gene的fpkm计算结果

四. Cuffmerge的使用

1. Cuffmerge简介

Cuffmerge将各个Cufflinks生成的transcripts.gtf文件融合称为一个更加全面的transcripts注释结果文件merged.gtf。以利于用Cuffdiff来分析基因差异表达。

2. 使用方法

$ cuffmerge [options]* <assembly_GTF_list.txt>
输入文件为一个文本文件,是包含着GTF文件路径的list。常用例子:
$ cuffmerge -o ./merged_asm -p 8 assembly_list.txt

3. 使用参数

-h | --help
-o <output_dir> default: ./merged_asm
将结果输出至该文件夹。

-g | --ref-gtf
将该reference GTF一起融合到最终结果中。

-p | --num-threads <int> defautl: 1
使用的CPU线程数

-s | --ref-sequence <seq_dir>/<seq_fastq>
该参数指向基因组DNA序列。如果是一个文件夹,则每个contig则是一个fasta文件;如果是
一个fasta文件,则所有的contigs都需要在里面。Cuffmerge将使用该ref-sequence来
帮助对transfrags分类,并排除repeats。比如transcripts包含一些小写碱基的将归类
到repeats.

4. Cuffmerge输出结果

输出的结果文件默认为 <output_dir>/merged.gtf

五. Cuffcompare的使用

1. Cuffcompare简介

Cuffcompare使用Cufflinks的GTF结果,对GTF结果进行比较。和reference gtf比较寻找novel转录子等。

2. Cuffcompare的使用方法

$ cuffcompare [options]* <cuff1.gtf> [cuff2.gtf] ... [cuffN.gtf]

使用例子:
$ cuffcompare -o cuffcmp cuff1.gtf cuff2.gtf

3. 使用参数

-h
-V
-o <outprefix> default: cuffcmp
输出文件的前缀

-r <reference_mrna.gtf>
参考的GFF文件。用来评估输入的gtf文件中gene models的精确性。每一个输入的gtf的is
oforms将和该参考文件进行比较,并被标注为 overlapping, matching 或 novel。

-R
当有了 -r 参数时,指定该参数时,将忽略参考GFF文件中的一些transcripts。这些tran
scripts不和任何输入的GTF文件overlapped。

-s <seq_dir>/<seq_fastq>
该参数指向基因组DNA序列。如果是一个文件夹,则每个contig则是一个fasta文件;如果是
一个fasta文件,则所有的contigs都需要在里面。小写字母的碱基用来将相应的transcri
pts作为repeats处理。

4. 输出结果

在当前目录下输出3个文件:<coutprefix>.stats,<coutprefix>.combined.gtf 和 <coutprefix>.tracking; 在输入的GTF的同目录下输出<cuff_in>.refmap 和 <cuff_in>.tmap 文件。

六. Cuffdiff的使用

1. Cuffdiff简介

用于寻找转录子表达的显著性差异。

2. Cuffdiff使用方法

$ cuffdiff [options]* <transcripts.gtf> <sample1_1.sam[,...,sample1_M.sam]> <sample2_1.sam[,...,sample2_M.sam]>...[sampleN_1.sam[,...,sampleN_M.sam]]
其中transcripts.gtf是由cufflinks,cuffcompare,cuffmerge所生成的文件,或是由其它程序生成的。

一个常用例子:
$ cuffdiff --lables lable1,lable2 -p 8 --time-series --multi-read-correct --library-type fr-unstranded --poisson-dispersion transcripts.gtf sample1.sam sample2.sam

3. 使用参数

-h | --help
-o | --output-dir <sting> default: ./
输出的文件夹目录。
-L | --lables <lable1,lable2,...,lableN>  default: q1,q2,...qN
给每个sample一个样品名

-p | --num-threads <int> default: 1
使用的CPU线程数

-T | --time-series
让Cuffdiff来按样品顺序来比对样品,而不是对所有的samples都进行两两比对。即第二个
SAM和第一个SAM比;第三个SAM和第二个SAM比;第四个SAM和第三个SAM比...

-N | --upper-quartile-form
使用75%分为数的值来代替总的值(比对到单一位点的fragments的数值),作normalize。
这样有利于在低丰度基因和转录子中寻找差异基因。

--total-hits-norm default: TRUE
Cufflinks在计算FPKM时,算入所有的fragments和比对上的reads。和下一个参数对立。
默认激活该参数。

--compatible-hits-norm
Cufflinks在计算FPKM时,只针对和reference transcripts兼容的fragments以及
比对上的reads。该参数默认不激活,只能在有 --GTF 参数下有效,并且作 RABT 或 ab
 initio 的时候无效。

-b | --frag-bias-correct
提供一个fasta文件来指导Cufflinks运行新的bias detection and correction 
algorithm。这样能明显提高转录子丰度计算的精确性。

-u | --multi-read-correct
让Cufflinks来做initial estimation步骤,从而更精确衡量比对到genome多个位点
的reads。

-c | --min-alignment-count <int>  default: 10
如果比对到某一个位点的fragments数目少于该值,则不做该位点的显著性分析。认为该位点
的表达量没有显著性差异。

-M | --mask-file <mask.(gtf/gff)>
提供GFF文件。Cufflinks将忽略比对到该GTF文件的transcripts中的reads。该文件中
常常是rRNA的注释,也可以包含线立体和其它希望忽略的transcripts的注释。将这些不需
要的RNA去除后,对计算mRNA的表达量是有利的。

-FDR <float> default: 0.05
允许的false discovery rate.

--library-type default:fr-unstranded
处理的reads具有链特异性。比对结果中将会有个XS标签。一般Illumina数据的library-
type为 fr-unstranded。

-m | --frag-len-mean default: 200
插入片段的平均长度。不过现在Cufflinks能learns插入片段的平均长度,因此不推荐自主
设置此值。

-s | --frag-len-std-dev default: 80
插入片段长度的标准差。不过现在Cufflinks能learns插入片段的平均长度,因此不推荐自
主设置此值。

--poisson-dispersion
Use the Poisson fragment dispersion model instead of learning one 
in each condition.

4. Cuffdiff输出

1. FPKM tracking files
2. Count tracking files
3. Read group tracking files
4. Differential expression test
5. Differential splicing tests – splicing.diff
6. Differential coding output – cds.diff
7. Differential promoter use – promoters.diff
8. Read group info – read_groups.info
9. Run info – run.info

七. cufflinks使用中遇到的问题

1. 使用cuffdiff时候,在最新版本下,无重复的RNA-seq样作比较,结果中没有差异表达基因?
在v2.0.1及之后的版本中cuffdiff貌似不支持无重复的RNA-seq数据了。使用之前的版本即可。