Archives

NGS数据质量过滤:Usearch fastq-filter、seqtk trimfq和seqtk_filter

本文介绍NGS数据(这里我们针对FASTQ序列数据)质量控制, 质量控制常见的内容有:

  1. 过滤低复杂度比例高的序列;
  2. 包含有’N’s的序列;
  3. 低质量区域,一般指5’端和3’端;
  4. phix序列;
  5. Expected errors 期望误差(这里按照这个字面解读);

其它几种比较常见, 这篇文章主要介绍期望误差值的方式。

期望误差可以参考:Expected errors predicted by Phred (Q) scores,简而言之就是将Phred值使用错误概率值(P_error)描述,对所有碱基的的P_error求和,该值我们可以称之为”E”值, 可以使用该值对序列进行过滤,文章: Error filtering, pair assembly and error correction for next-generation sequencing reads 描述的和详细, 并推荐使用E值阈值为 1.0

P_error 、Q值 和 ASCII 表映射关系如下

P_error

Robert Edgar 在其网站上也描述了为什么使用平均Q值过滤不好, 参考:Average Q score

下面介绍几个小工具,如何使用E值模式进行序列质量控制。

一、fastq-filter

fastq-filter为 Usearch 实现的序列质量控制工具,命令行参数介绍:cmd_fastq_filter

支持的命令行选项:

-fastqout           filename     FASTQ 输出文件;
-fastaout           filename     FASTA 输出文件;
-fastqout_discarded filename     过滤掉的序列 FASTQ 输出文件;
-fastaout_discarded filename     过滤掉的序列 FASTA 输出文件;
-relabel            prefix       重新设定序列输出标识符前缀,其后添加 .1, .2, .3 ...
-fastq_eeout                     序列标识符添加"ee=xx;" 信息。
                                 其值为切除低质量后计算得到的值;
-sample             string       序列标识符添加"sample=string" 信息;
-fastq_truncqual    N            从前往后,确定第一个碱基质量值小于N的位置,舍弃该碱基及其后所有碱基;
-fastq_maxee        E            过滤掉所有碱基EE值总和大于E的序列;
                                 其值为切除低质量后计算得到的值;
-fastq_trunclen     L            保留前N碱基,过滤掉长度小于L的序列;
-fastq_minlen       L            过滤掉长度小于L的序列;
-fastq_stripleft    N            删除左边N个碱基;
-fastq_maxee_rate   E            过滤掉满足条件:E > 总EE值/序列长度的序列;
-fastq_maxns        k            过滤掉序列包含多于k个"Ns"的序列;
-fastq_ascii        Q            64或者33指定Illumina质量值类型;

以上所有选项遵循这样原则: 先处理序列(主要是截断), 然后再计算特征值(比如EE值和序列长度)去判断过滤条件;

实例:

usearch -fastq_filter reads.fastq -fastq_trunclen 150 -fastq_maxee 0.5 -fastaout reads.fasta
usearch -fastq_filter reads.fastq -fastq_minlen 100 -fastq_truncqual 15 -fastqout filtered_reads.fastq

该工具对扩增子测序数据比较友好,比如我们测序数据质量不是很好,可以舍弃R2端序列,并使用`-fastq_truncqual` 截取质量值大于15的碱基,舍弃掉长度小于100的序列(`-fastq_minlen`), 如果5’端序列质量不好,或者检测到有bias, 可使用 -fastq_stripleft 删除左侧N个碱基。

二、seqtk trimfq

SEQTK 错误率概念,但和fastq-filter模式还不一样,使用trimfq进行质量控制。

命令行接口:

$ seqtk trimfq

Usage:   seqtk trimfq [options] <in.fq>

Options: -q FLOAT    error rate threshold (disabled by -b/-e) [0.05]
         -l INT      maximally trim down to INT bp (disabled by -b/-e) [30]
         -b INT      trim INT bp from left (non-zero to disable -q/-l) [0]
         -e INT      trim INT bp from right (non-zero to disable -q/-l) [0]
         -L INT      retain at most INT bp from the 5'-end (non-zero to disable -q/-l) [0]

trimfq执行模式有两种:

  1. 使用期望错误进行切除末端序列,保留最小长度(默认最小长度30);
  2. 直接切除左侧和右侧指定碱基数,保留最小长度为从5′->3’方向截取;

trimfq的error rate threshold执行模式为为寻找一个大于该阈值和小于该阈值的区间(该区间质量值总和最高),如error rate threshold不能保证有最小保留长度,则切换成最小长度为窗口的区间,寻找质量值总和最高的区间。

三、 DADA2 filterAndTrim

DADA2 提供了 DADA2 R函数也采用了maxEE模式对序列进行过滤;

filterAndTrim(fwd, filt, rev = NULL, filt.rev = NULL, compress = TRUE,
truncQ = 2, truncLen = 0, trimLeft = 0, maxLen = Inf, minLen = 20,
maxN = 0, minQ = 0, maxEE = Inf, rm.phix = TRUE, primer.fwd = NULL,
matchIDs = FALSE, id.sep = "\\s", id.field = NULL,
multithread = FALSE, n = 1e+05, OMP = TRUE, verbose = FALSE)

实例:

filt <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(240,160),
              maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
              compress=TRUE, multithread=TRUE) 

支持双端序列,可使用参数比较多:

文件信息:

fwd、filt、rev、filt.rev 为序列路径,filt和filt.rev为输出文件,会自动创建文件;
compress 支持压缩文件格式

切除操作:

truncQ  :Truncate reads at the first instance of a quality score less than or equal to truncQ。
truncLen:截取前指定个数碱基;
trimLeft:切除左侧指定长度,如提供truncLen,最终长度为truncLen-trimLeft

过滤参数:

minQ   : 碱基中如存在小于minQ值的序列被过滤;
maxLen : 最大序列长度阈值;
minLen : 最小序列长度阈值,所有trimming 操作后的有效长度;
maxEE  : 最大期望误差值;EE = sum(10^(-Q/10))
rm.phix: 过滤phix序列;

原则:过滤操作总是发生在切除操作后;

四、 Seqtk_filter

重新造的小轮子,解决几个小问题:

  1. 支持双端序列模式;
  2. 切除左侧低于指定质量值序列;
  3. 切除右侧低于指定质量值序列;
  4. 期望误差值参数过滤;EE = sum(10^(-Q/10))
  5. 期望误差率参数过滤;EE = sum(10^(-Q/10)), EE/length

命令行接口:

$ seqtk_filter
Usage: seqtk_fqchk [options] <in1.fq> <in2.fq> <prefix>
Options:
  -q INT     offset value: 33 for Sanger quality, 64 for Illumina Phred+64 quality, default: 33
  -E FLOAT   Discard reads >E expected errors. Calculated after any truncation options.
  -R FLOAT   Discard reads >E expected errors per base.Calculated after any truncation options.
  -L INT     Discard sequences with < L letters.
  -b INT     trim base quality <Q  from left  default:[0]
  -e INT     trim base quality <Q  from right default:[0]
  -v         print version number

实例:

seqtk_filter -L100 -E 1.5   -R 0.01  -e2 -b2   OS2_1.R1.fq.gz OS2_1.R2.fq.gz  OS

该程序现阶段功能实现数据质量控制功能。

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归 上海逻捷信息科技有限公司 所有。

Last Update: 2017-09-25 11:56 AM

Comments are closed.