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