高通量测序数据下机后得到了fastq的raw_data,通常测序公司在将数据返还给客户之前会做“clean”处理,即得到clean_data。然而,这些clean_data是否真的“clean”呢?
首先,我们应该做一下质控。如果质控不合格,就需要一些处理,比如去接头、去除量的reads。
(1)去除测序数据中的接头(用到的是fastx_toolkit里面的fastx_clipper工具):
Usage: fastx_clipper [-h] [-a ADAPTER] [-D] [-l N] [-n] [-d N] [-c] [-C] [-o] [-v] [-z] [-i INFILE] [-o OUTFILE] #去掉接头序列
[-a ADAPTER] =接头序列(默认为CCTTAAGG)
[-l N] = 忽略那些碱基数目少于N的reads,默认为5
[-d N] = 保留接头序列后的N个碱基默认 -d 0
[-c] = 放弃那些没有接头的序列.
[-C] = 只保留没有接头的序列.
[-k] = 报告只有接头的序列.
[-n] = 保留有N多序列,默认不保留
[-v] =详细-报告序列编号
[-z] =压缩输出.
[-D] = 输出调试结果.
[-M N] =要求最小能匹配到接头的长度N,如果和接头匹配的长度小于N不修剪
[-i INFILE] = 输入文件
[-o OUTFILE] = 输出文件
Example: fastx_clipper -a AGATCGGAAGAGCACACG -l 25 -d 0 -Q 33 -i SRR306394_1.fastq -o SRR306394_1_trimmed.fastq
(注:-l参数,即忽略那些碱基数目少于N的reads。默认的是5,而对于tophat来说,如果reads长度小于20,就会发出警告。因此,考虑到充分利用这些测序数据同时不损失比对的特异性,这里-l参数,取值为25。-Q参数,对于illumina测序来说,这个参数也是必须要加的,因为illumina测序在给单个碱基作质量打分的时候,加上了33,然后才转成的ASC II, 因此在这里需要加 -Q 33. 否则就会报这样的错误“fastx_clipper: Invalid quality score value (char '#' ord 35 quality value -29) on line 4.”)
(2)去除测序数据中的低质量reads(用到的是fastx_toolkit里面的fastq_quality_filter工具):
Usage:fastq_quality_filter [-h] [-v] [-q N] [-p N] [-z] [-i INFILE] [-o OUTFILE]过滤低质量序列
[-q N] = 最小的需要留下的质量值
[-p N] = 每个reads中最少有百分之多少的碱基需要有-q的质量值
[-z] =压缩输出
[-v] =详细-报告序列编号,如果使用了-o则报告会直接在STDOUT,如果没有则输入到STDERR
Example: fastq_quality_filter -q 20 -p 80 -Q 33 -i SRR306394_1.fastq -o SRR306394_1_filtered.fastq
(注:-q参数,即是指最小需要留下的质量值,这里需要保留Q20(99.99%)以上的,所以对应的是20。-p参数,即是指每个reads中最少有百分之多少的碱基需要有-q的质量值,这里我们要求的是一条reads中最少有80%以上的碱基质量值达到Q20我们才作保留,因此对应的是80. -Q参数,对于illumina测序来说,这个参数也是必须要加的,因为illumina测序在给单个碱基作质量打分的时候,加上了33,然后才转成的ASC II, 因此在这里需要加 -Q 33. 否则就会报这样的错误“fastx_clipper: Invalid quality score value (char '#' ord 35 quality value -29) on line 4.”)
(注:如果在做fastq_quality_filter之后,QC结果仍然可以看到reads的尾部有低质量的reads,此时可以将-p参数调大,比如调至90或者99)。
转载本文请联系原作者获取授权,同时请注明本文来自熊朝亮科学网博客。
链接地址:http://blog.sciencenet.cn/blog-1509670-914439.html