Quake的安装和使用

一. Quake简介

Quake是由CBCB(Center for Bioinformatics and Computational Biology)开发的运用于修正序列错误的软件。Quake采用k-mer的错误修正方式,特别适合于Illumina测序的short reads数据,将reads中的错误碱基进行修正,同时,必须要满足碱基的基因组覆盖度要>15X。其文章2010年发表在Genome Biology上。

二. Quake的安装

1. 下载并安装Boost
2. 下载Quake并安装

$ wget http://www.cbcb.umd.edu/software/quake/downloads/quake-0.3.4.tar.gz
$ tar zxf quake-0.3.4.tar.gz
$ cd Quake/src
$ make

3. 安装JELLYISH,并将jellyfish链接到quake

$ ln -s ..../jellyfish ../bin/jellyfish
这一步需要安装jellyfish,然后将其软链接到Quake的bin目录下。或者修改Quake的bin
目录下的quake.py脚本,将jellyfish所在的目录进行修正。

4. 安装R及R的软件包VGAM

$ R
> install.packages("VGAM")
> q(save="no")

5. 使用Quake的bin目录下的quake.py来运行程序

三. quake的使用参数

1. 主要参数:

-r READSF
    Fastq文件名
-f READS_LISTF
    一个文件,其中包含fastq文件名,每行一个文件名;如果是双末端reads,则是一行两
个文件名
-k K
    使用的k-mer的长度。如果基因组大小为G,则k-mer长度选择为: k ~= log(200G)
/log(4)
-p PROC  default: 4
    使用的CPU线程数
-q QUALITY_SCALE
    使用的碱基质量格式,一般是64或33.如果不给出,则软件会自行猜测

2. 计算k-mer参数:

--no_jelly
    使用quake自带的一个简单的程序来进行k-mer计数,而不是使用jellyfish。该程序
是单线程运行,速度相对慢,能满足像微生物基因组这样的小基因组测序,但是不适合于大的基
因组。
--no_count  default: false
    不进行k-mer计数,而其计数结果已经存于在目标文件[readsfile].qcts和[reads
file].cts中了。
--int  default:false
    对kmers以整数的方式进行计数,而不使用碱基质量值。默认情况下是利用到了碱基质量,
计数的称为qmer。
--hash_size=HASH_SIZE
    Jellyfish的参数,用来设置Hash的大小。如果不设置,Quake则会使用k来估计出一
个值来。

3. 覆盖度模型参数:

--no_cut  default: false
    Quake使用k-mer计数来画直方图,通过调用R脚本cov_model_qmer.r来画直方图,并
决定出Coverage cutoff的阈值。默认情况下是最优化的模型,从而得出一个cutoff的Co
verage值,此值将输出到文件cutoff.txt中。
--ratio=RATIO
    确定Coverage cutoff值的时候,该值处对应的qmer错误的可能性是正确的可能性的
倍数(默认是200倍),即正确率不足0.5%。该值越小,则阈值越松。

4. reads的修正参数:

-l MIN_READ
    输出的长度>=此值的修正后的reads。很多reads由于修正或修剪后会较短
--headers
    输出的fastq文件使用的头和原始文件一致,不包含修正信息
--log
    输出一个log文件,里面包含所有的修正日志,包括"碱基质量 位置 新的碱基 旧的碱基"。

四. Quake的结果

使用Quake对双末端测序文件 A.fastq 和 B.fastq 进行reads修正,则产生如下结果文件:

A.cor.fastq  B.cor.fastq
    修正过后的reads文件               
A.cor_single.fastq  B.err_single.fastq
    A文件中修正过后的reads结果;而B文件中的序列则是error reads
A.err_single.fastq  A.cor_single.fastq
    B文件中修正过后的reads结果;而A文件中的序列则是error reads
A.err.fastq  B.err.fastq
    A和B文件中成对的reads都是error reads

五. Quake的correct程序

Quake的bin文件中有一支C++的程序,用于reads的修正。其实Quake.py在运行过程会调用该程序,但是当出现基因组覆盖度较低,而cutoff的值出现异常结果,比如小于1时,可以考虑使用该程序来重新对reads进行修正,取cutoff值为1.0。其使用方法为:

$ correct -f [fastq list file] -k [k-mer size] -c [cutoff] -m [counts file] -p 24

jellyfish的安装和使用

一. JELLYFISH简介

JELLYFISHCBCB(Center for Bioinformatics and Computational Biology)的Guillaume Marçais 和 Carl Kingsford 研发的一款计数 DNA 的 k-mers 的软件。该软件运用 Hash 表来存储数据,同时能多线程运行,速度快,内存消耗小。该软件只能运行在64位的Linux系统下。其文章于2011年发表在杂志 Bioinformatics 上。

二. JELLYFISH安装

$ wget http://www.cbcb.umd.edu/software/jellyfish/jellyfish-1.1.10.tar.gz
$ tar zxvf jellyfish-1.1.10.tar.gz
$ mkdir jellyfish
$ cd jellyfish-1.1.10
$ ./configure --prefix=$PWD/../jellyfish
如果安装在当前目录中,会报错。
$ make -j 8
$ make install

三. JELLYFISH的使用

1. jellyfish的使用方法

jellyfish的功能有:kmer计数;融合二进制的Hash结果;统计Hash结果;通过Hash结果来画直方图;将Hash结果输出成文本格式;查询指定k-mer的数目。

$ jellyfish count [-o prefix] [-m merlength] [-t threads] [-s hashsize] [--both-strands] fasta [fasta ...]
$ jellyfish merge hash1 hash2 ...
$ jellyfish dump hash
$ jellyfish stats hash
$ jellyfish histo [-h high] [-l low] [-i increment] hash
$ jellyfish query hash
$ jellyfish cite

2. k-mer的计数

使用 count 的命令来执行计数功能,使用例子:

$ jellyfish count -m 16 -s 100M -t 24 -o mer_counts -c 7 input.fastq
使用fastq文件在默认参数上和fasta文件没有区别。生成的hash结果为二进制文件。

常用参数:

-m | --mer-len=<num>
    使用的k-mer的长度。如果基因组大小为G,则k-mer长度选择为: k ~= log(200G)
/log(4)。
-s | --size=<num>
    Hash 的大小。最好设置的值大于总的独特的(distinct)k-mer数,这样生成的文件只
有一个。若该值不够大,则会生成多个hash文件,以数字区分文件名。如果基因组大小为G,每
个reads有一个错误,总共有n条reads,则该值可以设置为『(G + n)/0.8』。该值识别 
M 和 G。
-t | --threads=<num>  default: 1
    使用的CPU线程数
-o | --output=<string>  default: mer_counts
    输出的结果文件前缀
-c | --counter-len=<num>  default:7
    k-mer的计数结果所占的比特数,默认支持的最大数字是2^7=128。对于基因组测序覆盖
度为N,则要使设置的该值要大于N。该值越大,消耗内存越大。
-out-counter-len=<num>  default:4
    输出的二进制hash文件中的计数结果所占的字节数,一个字节是8比特。则默认支持的最大
数字是2^32=4.3G
-C | --both-strand  default: false
    对正义链和反义链都进行计数
-q | --quake  default: false
    quake兼容模式
--quality-start=<num>  default: 64
    起始碱基质量的ASCII值,默认为PHRED64
--min-quality=<num>  default: 0
    支持的最小的碱基质量值,低于此值的碱基将由N代替
-L | --lower-count=<num>
    不输出数目低于此值的k-mer
-U | --upper-count=<num>
    不输出数目高于此值的k-mer

3. 融合二进制的输出结果

上一步的输出结果为二进制文件,可能输出了多个hash文件,因此需要将这些hash文件合并成一个文件,此时用到 merge 命令。使用方法:

$ jellyfish merge -o mer_counts_merged.jf hash1 hash2 ...

常用参数:

-o | --output=<string>  default: mer_counts_merged.jf
    输出的结果文件
--out-counter-len=<num>  default: 4
输出的二进制hash文件中的计数结果所占的字节数,一个字节是8比特。则默认支持的最大数字
是2^32=4.3G

4. 对hash结果进行统计

k-mer的结果以hash的二进制文件结果给出,需要统计出k-mer总数,特异的k-mer数目,只出现过一次的kmer数,出现了最多的k-mer的数目等信息。以stats命令来运行。使用方法:

$ jellyfish stats hash
示例结果为:
Unique:    32355544    #只出现过一次的k-mer的数目
Distinct:  88414020    #特异性的k-mer数目,包含上一个的数据
Total:     432232807   #总的k-mer数目
Max_count: 85348       #同一个k-mer出现的最多的数目

常用参数:

-L | --lower-count=<num>
    不统计数目低于此值的k-mer
-U | --upper-count=<num>
    不统计数目高于此值的k-mer

5. 通过Hash结果来画直方图

对k-mer的计数结果有个直观的认识,则需要统计出现了x(x=1,2,3…)次的kmer的数目y,以x,y为横纵坐标画出直方图。使用 histo 命令能给出 x 和 y 对应的值,将结果默认输出到标准输出。其使用方法为

$ jellyfish histo -l 1 -h 1000 hash

常用参数:

-l | --low=<num>  default: 1
    最低的 x 轴的值。同时结果会将低于此值的所有的k-mer的数目作为 (x-1) 的值。因
此该值为 2 和 1 的结果是一致的。
-h | --high=<num>  default: 10000
    最高的 x 轴的值。同时结果会将高于此值的所有的k-mer的数目的和作为 (x+1) 的值。
-i | --increment=<num>  default: 1
    x 轴取值是每隔该数值取值
-t | --threads=<num>  default: 1
    使用的CPU线程数
-f | --full  default: false
    全部的直方图

6. 将二进制Hash结果转换成文本文件

由于count命令生成的结果为二进制的,如有需要,则可以转换成可读文本文件。使用 dump 命令,使用方法:

$ jellyfish dump -c -t -U 1000 hash

常用参数:

-c | --colum  default: false
    生成结果为2列,第一列为k-mer序列,第二列为对应的数目。默认情况下是是fasta格
式,fasta的头为k-mer的数目,fasta的序列为k-mer的序列。
-t | --tab  default: false
    当 -c 参数存在时,以tab来进行分隔两行。默认是以空格来分开的。
-L | --lower-count=<num>
    不输出小于该值的k-mer
-U | --upper-count=<num>
    不输出高于该值的k-mer
-o | --output=<file>
    输出文件的路径和名称

7. 查询指定的k-mer出现的次数

如果需要从Hash结果中查询指定的k-mer出现的次数,则要是用 query 命令。从标准输入读取k-mer的序列,从标准输出得到k-mer对应的数目。使用方法

$ jellyfish query hash

常用参数:

-C | --both-strands  default: false
    同时查询k-mer序列的正负链
-i | --input=<file>
    输入的文件
-o | --output=<file>
    输出的文件

四. 思考

1. 对Illumina paired-end测序结果进行jellyfish分析

由于paired-end序列有一定的顺序,需要将第2个文件的序列进行反向重复后,在和第一个文件的序列合到一起进行分析。可以使用Trinity中附带的软件fastool来将fastq文件转换成fasta文件,以及反向重复的转换。

$ $Trinity_Home/trinity-plugins/fastool/fastool --illumina-trinity --to-fasta reads_1.fastaq > reads_1.fasta
$ $Trinity_Home/trinity-plugins/fastool/fastool --rev --illumina-trinity --to-fasta reads_2.fastaq > reads_2.fasta
$ cat reads_1.fasta reads_2.fasta > both.fasta
$ jellyfish count ....

ALLPATHS-LG的使用

一、ALLPATH简介

ALLPATHS-LG是一个基因组组装软件,适合于组装short reads数据,由Computational Research and Development group at the Broad Institute开发。ALLPATHS-LG是现在行业内公认进行基因组De novo组装效果最好的软件。

二. 基础注意事项

1. 不能只使用一个library数据进行组装;
2. 必须有一个"overlapping"的片段文库的paired-reads数据。比如,reads长度~
100bp,插入片段库长度~180bp;
3. 必须有jumping library数据;
4. 基因组组装需要100x或以上基因组覆盖度的碱基,这个覆盖度是指raw reads数据(在
error correction和filtering之前)的覆盖度;
5. 可以使用PacBio数据;
6. 不能使用454数据和Torrent数据。主要是这两者测序太贵,如果什么时候价格降低,有
需求的话,会写出相应的代码来满足要求;
7. 官方提供了测试用数据;
8. 不支持在整个计算机集群上进行运算;
9. 需要消耗的内存峰值大约是1.7bytes每个碱基,即输入10G的碱基数据量,大约需要17
G内存;
10. 对于试探性的参数,比如K,原则上可以调整。但是我们不会自行调整,并也不推荐。AL
LPATHS-LG不像其它De novo一样,Kmer大小的参数K和read大小之间没有直接的联系,
ALLPATHS-LG会在运行过程中运用一系列的K值。

三. ALLPATHS-LG使用方法

1. 基础的使用方法和命令

使用RunAllPathsLG这个命令来运行。虽然有很多参数,但是在没有指导的情况下不要随意使用,使用默认设置即可。其使用方法为:

$ RunAllPathsLG arg1=value1 arg2=value2 ...

参数主要是设置程序辨别的一些目录,在程序的运行过程,会输入相应目录中的数据,将结果输入到指定的目录。一个简单的命令使用例子:

#!/bin/sh

# ALLPATHS-LG needs 100 MB of stack space.  In 'csh' run 'limit stacksize 100000'.
ulimit -s 100000
# ALLPATHS-LG命令的写法与一般的linux参数写法不是很一样。采用 ‘参数=值’ 的方法,并使之成每行一个参数,使用'\'来连接各个参数,这样看起来直观易懂。初始接触的人可能会不适应。

RunAllPathsLG \
 PRE=$PWD\
 REFERENCE_NAME=species.genome\
 DATA_SUBDIR=data\
 RUN=run\
 SUBDIR=test\
 EVALUATION=STANDARD\
 TARGETS=standard\
 OVERWRITE=True\
 MAXPAR=8
 | tee -a assemble.out

2. 详细的参数说明

必须的参数
PRE (String)
    程序运行的根目录,所有的其它目录全在该目录下
REFERENCE_NAME (String)
    参考基因组目录名称,位于PRE目录下。如果有一个参考基因组,可将参考基因组放到该
目录中;若没有,则创建该文件夹用于基因组组装
DATA_SUBDIR (String)
    DATA子目录名称,位于REFERENCE_NAME目录下。程序从该目录中读取数据。
RUN (String)
    运行目录名称,位于DATA_SUBDIR下。程序将生成的中间文件和结果文件存储于该目录
。比如组装结果是一个名为ASSEMBLES的目录,位于该目录下。

部分可选参数:
SUBDIR (String) default: test
    子目录名,在REF/DATA/RUN/ASSEMBLIES目录下创建的存放基因组组装结果的目录
名。
K (int) default: 96
    核心Kmer大小,只有K=96能很好地运行。
EVALUATION (String: {NONE,BASIC,STANDARD,FULL,CHEAT})default:BASIC
    给定一个参考基因组,pipeline能在基因组组装的不同阶段对组装过程和结果进行评估。
    BASIC:基础评估,不需要参考基因组;
    STANDARD:使用参考基因组来运行评估模块;
    FULL:在某些组装模块下打开in-place评估,不会影响组装结果;
    CHEAT:稍微使用参考基因组指导组装,产生更详细的分析,能对组装结果产生小的(好方
向的)改变。
REFERENCE_FASTA (String) default: REF/genome.fasta
    评估中使用的参考基因组。
MAXPAR (int) default: 1 
    有些模块的运行是独立的,不相互依赖,能同时运行。该参数设定能同时运行的模块的最
大数目。由于pipeline中的绝大部分模块都能多线程运行,因此将该值设定大于1,效果不明
显。
THREADS (String) default: max 
    有些模块能多线程程运行,默认使用最大线程数运行。
OVERWRITE (Bool) default: False 
    是否覆盖存在的文件。可以设置该选项为True,在每次运行程序的时候设定RUN参数为
一个新的目录名,则比较好。
TARGETS (vec) default: standard 
    pipeline会生成一系列的文件,不同的文件的生成需要call不同的模块。如果某文件
已经存在了并且是最新的,则跳过相应的模块的运行。本参数指定生成哪些拟定的目标文件(p
seudo targets)。若目标文件没有相应的模块能生成,则会得到报错。
    none:没有拟定的目标文件,仅仅生成指定的目标文件;
    standard:生成组装文件和选定的评估文件;
    full_eval:生成组装文件和额外的评估文件。
TARGETS_REF (String) 
    在ref_dir目录中生成的目标文件。
    多个目标文件的书写方法为: TARGETS_REF="{target1,target2,target3}" 。
TARGETS_DATA (String) 
    在data目录中生成的目标文件。
TARGETS_RUN (String) 
    在run目录中生成的目标文件。
TARGETS_SUBDIR (String)
    在subdir中生成的目标文件。 
FORCE_TARGETS (Bool) default: False
    生成目标文件,即使文件已经存在并且看起来是很新的。

3. 输入文件与目录的准备

两个文库:插入片段长度为180bp和3000bp,illumina测序文件结果为fastq格式。以此为例来准备ALLPATHS-LG运行所需的文件和目录。

(1) 准备 in_groups.csv 和 in_libs.csv 文件。

这两个文件内容由逗号隔开,in_groups.csv文件内容如下:

group_name, library_name, file_name
firest, Illumina_180bp, seq/species_500bp_read?.fastq
second, Illumina_3000bp, seq/species_3000bp_read?.fastq

in_groups.csv文件的解释:

group_name:数据独特的代号,每一份数据有一个代号;
library_name:数据所属文库的名字,体现出该;
filename:数据文件所存放位置。可以为相对位置,文件名可以包含'*'和'?'(但是扩展名
中不能有该符号,因为要根据扩展名识别文件类型),从而代表paired数据。支持的文件类型有
'.bam','fasta','fa','fastq','fq','fastq.gz'和'fq.gz'。

in_libs.csv文件内容如下:

library_name, project_name, organism_name, type, paired, frag_size, frag_stddev, insert_size, insert_stddev, read_orientation, genomic_start, genomic_end
Illumina_180bp, species, species.genome, fragment, 1, 180, 10, , , inward, 0, 0
Illumina_3000bp, species, species.genome, jumping, 1, , , 3000, 500, outward, 0, 0

in_libs.csv文件的解释:

library_name:和in_groups.csv中的相匹配;
project_name:project的名字;
organism_name:测序物种的名字;
type:仅仅只是一个信息;
paired:0:Unpaired reads;1:paired reads;
frag_size:小片段文库插入片段长度的均值;
frag_stddev:小片段文库的插入片段长度估算的标准偏差;
insert_size:大片段文库插入片段长度的均值;
insert_stddev:大片段文库插入片段长度估算的标准偏差;
read_orientation:reads的方向,小片段文库为inward,大片段文库为outward;
genomic_start:reads从该位置开始,读入数据,如果不为0,之前的碱基都被剪掉;
genomic_end:reads从该位置开始,停止读入数据,如果不为0,之后的碱基都被剪掉。

(2) 使用PrepareAllPathsInputs.pl来对数据进行转换

ALLPATHS-LG接受的输入数据要求如下:

1. ALLPATHS-LG的输入数据支持小片段文库(fragment library)、大片段文库(jum
ping library)和超大片段文库(long jumping library)。并且前两种文库至少各有
一个才能进行基因组组装。超大片段文库是只插入片段>20kb的文库,其测序方向和小片段文
库一致,为inward。

2. ALLPATHS-LG的输入数据放置在//文件夹下,包含3种文件:碱基文件,质量文件和配
对信息文件
   frag_reads_orig.fastb
   frag_reads_orig.qualb 
   frag_reads_orig.pairs 

   jump_reads_orig.fastb 
   jump_reads_orig.qualb 
   jump_reads_orig.pairs

以下是可选的超大插入片段文库对应的数据文件(非必须):

  long_jump_reads_orig.fastb 
  long_jump_reads_orig.qualb 
  long_jump_reads_orig.pairs

使用PrepareAllPathsInputs.pl来将fastq等格式的测序结果转换成ALLPATHS-LG可接受的文件。以下是该程序的参数:

DATA_DIR
    将转换后的数据文件放到此文件夹下。
PICARD_TOOLS_DIR
    若输入数据为bam格式,则需要用到Picard软件,该参数Picard的路径
IN_GROUPS_CSV
    输入的in_groups.csv文件名
IN_LIBS_CSV
    输入的in_libs.csv文件名
INCLUDE_NON_PF_READS default: 1
    1:包含non-PF reads;0:仅仅只包含PF reads.
PHRED_64 default: 0
    0:碱基质量是ASCII的33到126,一般情况下Illumina数据的最低碱基质量是'B';
1:碱基质量的ASCII码是从64到126,一般情况下Illumina数据的最低碱基质量是'#'。
PLOIDY
    生成ploidy文件。该文件就包含一个数字 1 或者 2 。1表示基因组为单倍体型,2表
示双倍体型。
HOSTS
    列出平行forking的host主机(这些主机必须要能无密码直接ssh连上)。比如“2,3.
host2,4.host3"表示使用本地机器的2个CPU线程,host2机器的3个CPU线程和host3机
器的4个CPU线程。

以下是不常用的参数,主要用来选择转换的数据量的大小。当测序数据量太多,而只想使用其
中一部分数据的时候,可以用到
FRAG_FRAC
    使用小片段库reads的比例。比如 30% 或 0.3 。如果设定了此值,则不能同时设定
FRAG_COVERAGE。
JUMP_FRAC
    使用大片段库reads的比例。比如 20% 或 0.2 。如果设定了此值,则不能同时设定
JUMP_COVERAGE。
LONG_JUMP_FRAC
    使用超大片段库reads的比例。 比如 90% 或 0.9 。如果设定了此值,则不能同时
设定LONG_JUMP_COVERAGE。
GENOME_SIZE
    估计的基因组大小,用来计算对应覆盖度所对应的reads数
FRAG_COVERAGE
    所期望的小片度库的覆盖度,比如 45. 要求GENOME_SIZE有设定
JUMP_COVERAGE
    所期望的大片度库的覆盖度,比如 45. 要求GENOME_SIZE有设定
LONG_JUMP_COVERAGE
    所期望的超大片度库的覆盖度,比如 1. 要求GENOME_SIZE有设定

一个 PrepareAllPathsInputs.pl 的例子如下:

#!/bin/sh

# ALLPATHS-LG needs 100 MB of stack space.  In 'csh' run 'limit stacksize 100000'.
ulimit -s 100000

mkdir -p species.genome/data

# NOTE: The option GENOME_SIZE is OPTIONAL. 
#       It is useful when combined with FRAG_COVERAGE and JUMP_COVERAGE 
#       to downsample data sets.
#       By itself it enables the computation of coverage in the data sets 
#       reported in the last table at the end of the preparation step. 

# NOTE: If your data is in BAM format you must specify the path to your 
#       picard tools bin directory with the option: 
#
#       PICARD_TOOLS_DIR=/your/picard/tools/bin

PrepareAllPathsInputs.pl\
 DATA_DIR=$PWD/species.genome/data\
 PLOIDY=1\
 IN_GROUPS_CSV=in_groups.csv\
 IN_LIBS_CSV=in_libs.csv\
 OVERWRITE=True\
 | tee prepare.out

(3) 使用Fasta2Fastb来转换得到参考基因组的输入数据

将参考基因组的fasta文件改名为genome.fasta并放到PRE/REFRENCE_NAME/文件夹下。使用Fasta2Fastb来将fasta文件转换为其二进制文件,在PRE/REFRENCE_NAME/目录下生成genome.fastb文件。

$ Fasta2Fastb IN=PRE/REFRENCE_NAME/genome.fasta

4. ALLPATHS Cache的使用

PrepareAllPathsInputs.pl脚本实际上是包含两个步骤:
1. 将测序的原始文本数据(fastq等文件)转换成二进制数据(fastb,qualb文件),并将各个group测序数据的二进制数据放置到缓存文件夹:<DATA>/read_cache。
2. 将部分或者全部的二进制数据结合并生成所需要的运用于基因组组装的数据文件。因此,使用ALLPATHS Cache的优点:a. 当加入新的测序数据或进行不同数据量的基因组组装评估时,需要再次进行基因组组装,这时可以仅仅只转换新的测序数据,节约了重新进行数据转换的时间。只需要将二进制数据根据需要来生成ALLPATHS-LG接受的文件; b. 当测序数据的测序质量有phred64和phred33兼有的时候,则需要使用Cache来分别转换数据。

(1). 使用CacheLibs.pl将文库信息读入到Cache

使用方法:

CacheLibs.pl\
  CACHE_DIR=<CACHE_DIR>\
  IN_LIBS_CSV=in_libs.csv\
  ACTION=Add

使用参数:

CACHE_DIR
    缓存文件夹的绝对路径名
IN_LIBS_CSV
    输入的in_libs.csv文件名
ACTION  deault: List
    Add,List 或者 Remove,从Cache中添加,列出或移除文库信息。

(2). 使用CacheGroups.pl将文本数据文件转换成二进制数据文件到Cache

使用方法:

CacheGroups.pl\
  CACHE_DIR=<CACHE_DIR>\
  PICARD_TOOLS_DIR=/opt/picard/bin\
  IN_GROUPS_CVS=in_groups.csv\
  TMP_DIR=/tmp\
  HOSTS='2,3.host2,4.host3'\
  ACTION=Add

使用参数:

CACHE_DIR
    缓存文件夹的绝对路径名
PICARD_TOOLS_DIR
    若输入数据为bam格式,则需要用到Picard软件,该参数为Picard的路径
IN_GROUPS_CSV
    输入的in_groups.csv文件名
PHRED_64  default: 0
    输入的fastq的碱基质量格式,True: 碱基质量格式为PHRED64; False: 碱基质量
格式为PHRED33.
OVERWRITE  default: 0
    是否覆盖已存在的数据文件
TMP_DIR
    临时文件夹的路径,如果数据量够大,则必须要够大
INCLUDE_NON_PF_READS default: 1
    1:包含non-PF reads;0:仅仅只包含PF reads.
HOSTS
    列出平行forking的host主机(这些主机必须要能无密码直接ssh连上)。比如“2,3.
host2,4.host3"表示使用本地机器的2个CPU线程,host2机器的3个CPU线程和host3机
器的4个CPU线程。
ACTION  default: List
    Add,List 或者 Remove,从Cache中添加,列出或移除group的数据信息。

(3). 使用CacheToAllPathsInputs.pl来生成ALLPATHS的输入数据。
使用方法:

CacheToAllPathsInputs.pl\
  CACHE_DIR=<CACHE_DIR>\
  DATA_DIR=DATA_DIR\
  GROUPS="{group1,group2}"\
  FRAG_FRAC=100%\
  JUMP_FRAC=100%

常用参数:

CACHE_DIR
    缓存文件夹的绝对路径名
DATA_DIR
    DATA文件夹就绝对路径名
GROUPS
    需要转换的到ALLPATHS的输入数据的groups名。格式为"{group1,group2,group
...}"
IN_GROUPS_CSV
    in_groups.csv文件名,和 GROUPS 参数二选一
FRAG_FRAC
    使用小片段库reads的比例。比如 30% 或 0.3 。如果设定了此值,则不能同时设定
FRAG_COVERAGE。
JUMP_FRAC
    使用大片段库reads的比例。比如 20% 或 0.2 。如果设定了此值,则不能同时设定
JUMP_COVERAGE。
LONG_JUMP_FRAC
    使用超大片段库reads的比例。 比如 90% 或 0.9 。如果设定了此值,则不能同时
设定LONG_JUMP_COVERAGE。
FRACTIONS
    同时设置上述3个参数。比如"{0.5,30%,100%}"
GENOME_SIZE
    估计的基因组大小,用来计算对应覆盖度所对应的reads数
FRAG_COVERAGE
    所期望的小片度库的覆盖度,比如 45. 要求GENOME_SIZE有设定
JUMP_COVERAGE
    所期望的大片度库的覆盖度,比如 45. 要求GENOME_SIZE有设定
LONG_JUMP_COVERAGE
    所期望的超大片度库的覆盖度,比如 1. 要求GENOME_SIZE有设定
COVERAGES
    同时设置上述3个参数,比如"{45,50,2}"
LONG_READ_MIN_LEN  default: 500
    设置被称为long unpaired read的阈值(适用于PacBio reads)
PLOIDY  
    生成ploidy文件。该文件就包含一个数字 1 或者 2 。1表示基因组为单倍体型,2表
示双倍体型。如果没有该参数,则不生成ploidy文件

四. 思考题

对一个简单的真菌物种,运用Illumina Hiseq2000平台,构建了180bp,500bp,3000bp长度的3个DNA片段文库,分别进行了测序。获得了相应的fastq数据文件,如何使用ALLPATHS-LG来进行基因组的De novo组装?

1. 安装ALLPATHS-LG到Unix系统服务器上

2. 准备测序的数据文件

测序的数据文件如下,并放置到当前工作目录下的 seq 文件夹中。

180.reads1.fastq    180.reads2.fastq
500.reads1.fastq    500.reads2.fastq
3000.reads1.fastq   3000.reads2.fastq

3. in_libs.csv 和 in_groups.csv 文件的准备
两个文件都放置到当前工作目录下。 in_libs.csv 内容如下:

library_name, project_name, organism_name, type, paired, frag_size, frag_stddev, insert_size, insert_stddev, read_orientation, genomic_start, genomic_end
Illumina_180bp, species, species.genome, fragment, 1, 180, 20, , , inward, 0, 0
Illumina_500bp, species, species.genome, fragment, 1, 500, 50, , , inward, 0, 0
Illumina_3000bp, species, species.genome, jumping, 1, , , 3000, 500, outward, 0, 0

in_groups.csv 内容如下:

group_name, library_name, file_name
180, Illumina_180bp, ./seq/180.reads?.fastq
500, Illumina_500bp, ./seq/500.reads?.fastq
3000, Illumina_3000bp, ./seq/3000.reads?.fastq

4. 使用PrepareAllPathsInputs.pl来对数据进行转换

将以下内容写入prepare.sh并运行

#!/bin/sh
ulimit -s 100000
mkdir -p species.genome/data
PrepareAllPathsInputs.pl\
 DATA_DIR=$PWD/species.genome/data\
 PLOIDY=1\
 IN_GROUPS_CSV=in_groups.csv\
 IN_LIBS_CSV=in_libs.csv\
 OVERWRITE=True\
 | tee prepare.out

5. 运行ALLPATHS-LG主程序进行基因组的De novo组装

将以下内容写入到assemble.sh中,并运行

#!/bin/sh
ulimit -s 100000
RunAllPathsLG\
 PRE=$PWD\
 REFERENCE_NAME=species.genome\
 DATA_SUBDIR=data\
 RUN=run\
 SUBDIR=test\
 OVERWRITE=True\
 MAXPAR=8\
 | tee -a assemble.out

6. 结果文件

基因组组装结果文件位于./species.genome/data/run/ASSEMBLIES/test/final.assembly.fasta文件。同一个目录下有一个文件final.assembly.efasta也是组装结果。
efasta格式即“enhanced fasta”。
efasta和fasta的区别:

fasta: ATGTCNTGTCG
efasta:ATGTC{A,T}GTCG

efasta使结果更明确,但是在数据处理的时候,不易兼容。

Trinity的安装与使用

一、 Trinity简介

Trinity,是由 the Broad Institute 开发的转录组de novo组装软件,由三个独立的软件模块组成:Inchworm,Chrysalis和Butterfly。三个软件依次来处理大规模的RNA-seq的reads数据。Trinity的简要工作流程为:
Inchworm: 将RNA-seq的原始reads数据组装成Unique序列;
Chrysalis: 将上一步生成的contigs聚类,然后对每个类构建Bruijn图;
Butterfly: 处理这些Bruijn图,依据图中reads和成对的reads来寻找路径,从而得到具有可变剪接的全长转录子,同时将旁系同源基因的转录子分开。
Trinity发表在 Nature Biotechnology

二、 Trinity的安装

1. 下载Trinity
2. 安装Trinity。

$ tar zxvf trinityrnaseq_r2013-02-25.tgz
$ cd trinityrnaseq_r2013-02-25
$ make
仅需要在安装目录下进行make即可。该命令编译了由C++编写的Inchworm和Chrysalis,而
使用Java编写的Butterfly则不需要编译,可以直接使用。

三、Trinity的使用

1. 直接运行安装目录下的程序Trinity.pl来使用该软件,不带参数则给出使用帮助。其典型用法为:

Trinity.pl --seqType fq --JM 50G --left reads_1.fq  --right reads_2.fq --CPU 6

2. Trinity参数

必须的参数:
--seqType     reads的类型:(cfa, cfq, fa, or fq)
--JM          jellyfish使用多少G内存用来进行k-mer的计算,包含‘G’这个字符
--left        左边的reads的文件名
--rigth       右边的reads的文件名
--single      不成对的reads的文件名

可选参数:

Misc:
--SS_lib_type        reads的方向。成对的reads: RF or FR; 不成对的reads
: F or R。在数据具有链特异性的时候,设置此参数,则正义和反义转录子能得到区分。默认
情况下,不设置此参数,reads被当作非链特异性处理。FR: 匹配时,read1在5'端上游, 
和前导链一致, read2在3'下游, 和前导链反向互补. 或者read2在上游, read1在下游反
向互补; RF: read1在5'端上游, 和前导链反向互补, read2在3'端下游, 和前导链一致;
--output             输出结果文件夹。默认情况下生成trinity_out_dir文件夹并
将输出结果保存到此文件夹中。
--CPU                使用的CPU线程数,默认为2
--min_contig_length  报告出的最短的contig长度。默认为200
--jaccard_clip       如果两个转录子之间有UTR区重叠,则这两个转录子很有可能在
de novo组装的时候被拼接成一条序列,称为融合转录子(Fusion Transcript)。如果有
fastq格式的paired reads,并尽可能减少此类组装错误,则选用此参数。值得说明的是:
1. 适合于基因在基因组比较稠密,转录子经常在UTR区域重叠的物种,比如真菌基因组。而对
于脊椎动物和植物,则不推荐使用此参数; 2. 要求fastq格式的paired reads文件(文件
中reads名分别以/1和/2结尾,以利于软件识别),同时还需要安装bowtie软件用于reads
的比对; 3. 单独使用具有链特异性的RNA-seq数据的时候,能极大地减少UTR重叠区很小的
融合转录子; 4. 此选项耗费运算,若没必要,则不用此参数。
--prep               仅仅准备一些文件(利于I/O)并在kmer计算前停止程序运行
--no_cleanup         保留所有的中间输入文件
--full_cleanup       仅保留Trinity fasta文件,并重命名成${output_dir}.
Trinity.fasta
--cite               显示Trinity文献引证和一些参与的软件工具
--version            报告Trinity版本并推出

Inchworm 和 K-mer 计算相关选项:
--min_kmer_cov      使用Inchworm来计算K-mer数量时候,设置的Kmer的最小值。
默认为1
--inchworm_cpu      Inchworm使用的CPU线程数,默认为6和--CPU设置的值中的
小值。

Chrysalis相关选项:
--max_reads_per_graph   在一个Bruijn图中锚定的最大的reads数目,默认为200
000
--no_run_chrysalis      运行Inchworm完毕,在运行chrysalis之前停止运行
Trinity
--no_run_quantifygraph  在平行化运算quantifygrahp前停止运行Trinity

Butterfly相关选项:
--bfly_opts                    Butterfly额外的参数
--max_number_of_paths_per_node 从node A -> B,最多允许多少条路径。默认
为10
--group_pairs_distance         最大插入片读长度,默认为500
--path_reinforcement_distance  延长转录子路径时候,reads间最小的重叠碱基
数。默认PE:75; SE:25
--no_triplet_lock              不锁定triplet-supported nodes
--bflyHeapSpaceMax             运行Butterfly时java最大的堆积空间,默认
为20G
--bflyHeapSpaceInit            java初始的堆积空间,默认为1G
--bflyGCThreads                java进行无用信息的整理时使用的线程数,默
认由java来决定
--bflyCPU                      运行Butterfly时使用的CPU线程数,默认为2
--bflyCalculateCPU             计算Butterfly所运行的CPU线程数,由公式
 80% * max_memory / maxbflyHeapSpaceMax 得到
--no_run_butterfly             在Chrysalis运行完毕后,停止运行Butterfly

Grid-computing选项:
--grid_computing_module   选定Perl模块,在/Users/bhaas/SVN/trinityr
naseq/trunk/PerlLibAdaptors/。

3. 适合于illumina测序数据的真菌物种转录组组装的Trinity命令为:

Trinity.pl --seqType fq --JM 50G --left reads_1.fq  --right reads_2.fq --SS_lib_type FR --output transcriptome_tissue --CPU 24 --jaccard_clip --inchworm_cpu 24 --group_pairs_distance 500 --bflyCPU 24

4. Trinity生成的结果文件
运行程序结束后,转录组结果为trinity_out_dir/Trinity.fasta。可以使用软件所带的一支程序分析转录组统计信息。

$ $TRINITY_HOME/util/TrinityStats.pl trinity_out_dir/Trinity.fasta
Total trinity transcripts:	30706
Total trinity components:	26628
Contig N50: 554

三. Trinity运行原理与过程

1. 检测java的可运行性,因为buttfly会用到

2. 运行jellyfish,使用其dump命令得到jellyfish.kmers.fa文件

3. Inchworm(Linear contig construction from k-mers)

assembles the RNA-seq data into the unique sequences of transcripts, often generating full-length transcripts for a dominant isoform, but then reports just the unique portions of alternatively spliced transcripts.

4. Chrysalis

clusters the Inchworm contigs into clusters and constructs complete de Bruijn graphs for each cluster. Each cluster represents the full transcriptonal complexity for a given gene (or sets of genes that share sequences in common). Chrysalis then partitions the full read set among these disjoint graphs.

5. Butterfly

then processes the individual graphs in parallel, tracing the paths that reads and pairs of reads take within the graph, ultimately reporting full-length transcripts for alternatively spliced isoforms, and teasing apart transcripts that corresponds to paralogous genes.

四. 注意事项

3.1 Trinity分步运行

当数据量比较大的时候,trinity运行的时间会很长,同时,内存不够等情况出现的时候有可能程序运行崩溃。最好是分步运行。下一步会接着前一步进行下去。

Stage 1: generate the kmer-catalog and run Inchworm: –no_run_chrysalis

Stage 2: Chrysalis clustering of inchworm contigs and mapping reads: –no_run_quantifygraph

Stage 3: Chrysalis deBruijn graph construction: –no_run_butterfly

Stage 4: Run butterfly, generate final Trinity.fasta file. (exclude –no_ options)

3.2 计算资源

Ideally, you will have access to a large-memory server, ideally having ~1G of RAM per 1M reads to be assembled (but often, much less memory may be required).

The assembly from start to finish can take anywhere from ~1/2 hour to 1 hour per million reads (your mileage may vary). 个人记录了一次,使用dell服务器,64GB RAM,24 threads : 53M 的reads,运行了16.5h(平均3.2M/h),内存使用峰值为43G.

ALLPATHS-LG的安装

ALLPATHS-LG的安装步骤

1. 下载GMP并安装。可以跳过1,2,3步。

下载的压缩包为lz格式,必须使用lzip解压缩。
$ sudo yum install lzip
$ lzip -d gmp-5.1.1.tar.lz
$ tar xvf gmp-5.1.1.tar
$ cd gmp-5.1..1
$ ./configure
$ make
$ sudo make install

2. 下载MPFR并安装。
3. 下载MPC并安装
4. 下载gcc并安装。

直接在源目录中configure并make易发生错误;
若有错误"fatal error: gnu/stubs-32.h 没有那个文件或目录",则需安装glibc-
devel.i686(centos6.3)或glibc-devel.i386(centos5.8)。

$ sudo yum install glibc-devel.i686
$ tar -jxvf gcc-4.7.2.tar.bz2
$ cd gcc-4.7.2
$ ./contrib/download_prerequisites
下载上述的三个所需的库,并创造链接。因此不需要一个个下载并安装那么麻烦。

$ cd ../
$ mkdir gcc-build
$ cd gcc-build
$ ../gcc-4.7.2/configure --prefix=/opt/gcc-4.7.2
推荐不安装在默认路径,不然不宜卸载。现在gcc取消了make uninstall命令。默认路径
是/usr/local;系统自带的gcc路径为/usr。

$ make -j 8
多线程运行,数字依机器配置而定,不然默认情况下极度耗费时间。

$ sudo make install
$ sudo make clean (optional)

$ vim /etc/ld.so.conf.d/gcc-4.7.2.x86_64.conf
在文本文件中添加以下两行
/opt/gcc-4.7.2/lib
/opt/gcc-4.72/lib64
$ sudo ldconfig
将库文件加入到访问路径,并更新库。

5. 下载ALLPATHS-LG并安装。

新版的ALLPATHS-LG需求gcc 4.7.0及以上版本,安装gcc 4.7.0及以上版本又需要gmp,
mpfr和mpc具有较高的版本。

$ tar zxvf LATEST_VERSION.tar.gz
$ cd allpathslg-?????
$ ./configure --prefix=/home/chenlianfu/programs/ALLPATHS-LG/ CXX
=/opt/gcc-4.7.2/bin/c++ CXXPP=/opt/gcc-4.7.2/bin/cpp
$ make
$ make install

InterProScan的三种使用方法

Interproscan,通过蛋白质结构域和功能位点数据库预测蛋白质功能。是EBI开发的一个集成了蛋白质家族、结构域和功能位点的非冗余数据库。Interproscan整合了一些使用最普及的一些数据库,并应用于功能未知的蛋白进行Interpro注释和GO注释。
以下介绍3中interpro注释的方法:

一、网页版的Interpro注释

打开InterProScan的官网地址:http://www.ebi.ac.uk/Tools/pfa/iprscan/。将序列粘贴到输入框中进行Interpro注释。

优点:使用网页版,方便快捷;不消耗本地计算资源。
缺点:输入必须为蛋白质序列;一次只能比对条蛋白质序列。

二、使用EBI提供的perl程序进行远程比对

程序下载网页地址:http://www.ebi.ac.uk/Tools/webservices/services/pfa/iprscan_rest
其实,除了perl程序,Python和Ruby也各有一支程序。分别是:
iprscan_lwp.pl ; iprscan_urllib2.py ; iprscan_net_http.rb

优点:不消耗本地计算资源;可以自己编写脚本来大批量比对本地的protein序列到EBI的
interpro服务器,获得得interpro注释。
缺点:比对结果为xml格式,普通科研人员玩不来。

三、本地化的InterProScan注释

3.1 本地化的InterProScan安装与配置

3.1.1 从ftp://ftp.ebi.ac.uk/pub/databases/interpro/iprscan下载以下5个文件:

RELEASE/latest/iprscan_v4.8.tar.gz
BIN/4.x/iprscan_bin4.x_[PLATFORM].tar.gz 
DATA/iprscan_DATA_[LATESTDATAVERSION].tar.gz
DATA/iprscan_PTHR_DATA_[LATESTDATAVERSION].tar.gz
DATA/iprscan_MATCH_DATA_[LATESTDATAVERSION].tar.gz

3.1.2 将5个文件解压到一个文件夹中,然后运行其中的文件Config.pl,来对InterProScan进行配置。
3.1.3 配置的过程中,若选择进行本地web配置,则修改本地www服务的配置文件,以能进行本地化网页版的运行。

3.2 本地化InterProScan的使用。

3.2.1 命令行运行iprscan的方法:

$ bin/iprscan -cli -iprlookup -goterms -format xml -i test.fasta -o test.out

3.2.2 iprscan的参数说明:

-cli         设定程序在unix命令下运行,如果不设此参数,程序会被当作CGI程序运行。
-iprlookup   结果里显示相应的interpro注释信息。
-goterms     结果里显示相应的GO注释信息,但前面要加上-iprlookup参数。
-format      输出结果的格式,有raw, xml, txt, html(default), ebixml
(EBI header on top of xml) gff。
-appl  数据库和扫描方法的的选择。无此参数表示默认选择全部数据库(配置Interpro
scan时候设置的数据库)。选择多个数据库则需该参数多次。本地运行Interproscan相
比官网运行,有两个数据库不能选择为:tmhmm 和 signalp。这两个数据库的选择需要
commercial license。
-i           输入文件,InterProScan支持输入蛋白质序列和核酸序列,如果输入核
酸系列,程序会将其翻译成蛋白质序列,你可以指定翻译用到的密码表,用下面的-trtable
参数,序列格式可以是raw,Fasta或者EMBL。
-o           结果输出文件,如果不选择此参数,结果将输出到标准输出,输出格式可以
用下面的-format参数设定。
-trtable   选择核酸翻译蛋白质的密码表,同时可以设定-trlen 参数来控制核酸翻译
的转录子长度。
-nocrc       不对输入蛋白质序列进行crc64匹配。不加此参数,则默认是会对蛋白质
序列开启了crc64匹配。Interpro数据库(memember database)已经包含了大量序列
搜索的结果,就是如果你的蛋白序列已经包含在interpro的数据库里面,iprscan会直接
给出搜索结果,无需进行本地运算。interpro数据库不包括tmhmm, coil和signalp,
所以crc64匹配不到这3个数据库。一般情况下,由于commercial licese而无法使用
tmhmm和signalp数据库,coil数据库也不会使用。故不使用此参数,加快程序运行速度。
-seqtype     输入的序列类型,蛋白质序列(-seqtype p)(defult)或者(-se
qtype n)。
-email       设定一个Email地址,程序运行完毕向信箱发送邮件通知分析完毕。
-verbose     程序运行过程中显示运行的状态。
-help        显示帮助信息。

3.2.3 InterProScan其它附带的重要程序

meter.pl     reports the progress of a job.有百分之几的chunk已经运行完毕。
converter.pl 将raw的格式转换成其它的格式,比如html,xml,txt等格式。
iterator.pl  运用于逐条去注释序列

3.2.4 多线程运行
hmmpfam, hmmscan 和 hmmsearch 能多线程运行。
经过测试 PIR superfamily 和 SUPERFAMILY 这两个数据库的应用比较费时,可以设置conf文件夹中的配置文件hmmpir.conf和superfamily.conf,将其中的cpu_opt的值设置高一些。默认是1。
3.2.5 逐条序列地运行
InterProScan不能立马给出结果文件并相继把结果放入到结果文件中。可以采用逐条比对的方法来得到注释结果。可以采用如下的方法来随时终结掉程序,并拿到部分结果。

$ bin/iterator.pl -i test.fa -o test.out -c "bin/iprscan -cli -i %infile -iprlookup -goterms -format xml"

3.2.6 优缺点

优点:使用本地化的数据库,在断网和计算机资源充足的情况下,能加快注释速度;本地化网
页版能同时比对多条序列;本地化能对DNA序列进行interpro注释。
缺点:本地化安装InterProScan比较复杂耗时;需要不时更新本地数据库;本地化运行耗
费计算资源大;

SAM格式简介

查看文献请点击:The Sequence Alignment/Map format and SAMtools。也可以google搜索“The SAM Format Specification”查询SAM Format的详细文档。
SAMtools作为一个通用工具,使用各种方法来,以SAM格式来对reads的比对结果进行预处理。比如:索引排序,查找变异和比对查看等。
SAM格式由两部分组成:头部区和比对区,都以tab分列。

头部区:以‘@’开始,体现了比对的一些总体信息。比对的SAM格式版本,比对的参考序列,
比对使用的软件等。
比对区: 比对结果,每一个比对结果是一行,有11个主列和1个可选列。
@HD VN:1.0 SO:unsorted  
头部区第一行:VN是格式版本;SO表示比对排序的类型,有unkown(default),unsorted,
queryname和coordinate几种。samtools软件在进行排序后不能自动更新bam文件的SO值。
picard却可以。
@SQ SN:A.auricula_all_contig_1 LN:9401
参考序列名。这些参考序列决定了比对结果sort的顺序。SN是参考序列名;LN是参考序列
长度;
@RG ID:sample01
Read Group. 1个sample的测序结果为1个Read Group;该sample可以有多个library
的测序结果。
@PG ID:bowtie2 PN:bowtie2 VN:2.0.0-beta7
比对所使用的软件。
比对区11个列和可选列的解释
1  QNAME  比对的序列名
2  FLAG   Bwise FLAG(表明比对类型:pairing,strand,mate strand等)
3  RNAME  比对上的参考序列名
4  POS    1-Based的比对上的最左边的定位
5  MAPQ   比对质量
6  CIGAR  Extended CIGAR string (操作符:MIDNSHP) 比对结果信息:匹配碱基数,可变剪接等。
7  MRNM   相匹配的另外一条序列,比对上的参考序列名
8  MPOS   1-Based leftmost Mate POsition
9  ISIZE  插入片段长度
10 SEQ    和参考序列在同一个琏上的比对序列(若比对结果在负意链上,则序列是其反向重复序列)
11 QUAL   比对序列的质量(ASCII-33=Phred base quality)
12 可选的行,以TAG:TYPE:VALUE的形式提供额外的信息
比对区第2列FLAG说明:
Bitwise(按位计算)
Bit(0x表示十六进制) 十进制数值 Description
0x1    1   比对时有多个种子序列比对上

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数据了。使用之前的版本即可。