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 ....