测序数据质控和预处理之fastp

2019-10-05 本文已影响0人吴十三和小可爱的札记

1. 引言

下机的FASTQ数据通常需要进行质控和预处理，以保证下游分析输的准确性。fastp软件仅扫描数据文件一次，就可以完成FASTQC + cutadapt + Trimmomatic 的功能；而且它使用C++开发，利用高效算法，并且支持多线程，加快了处理速度。

2. 特点

对输入和输出数据进行详尽的质量剖析，生成人性化的报告。(quality curves, base contents, KMER, Q20/Q30, GC Ratio, duplication, adapter contents...)
过滤掉低质量，太短和太多N的"bad reads"。
从头部(5')或尾部(3')计算滑动窗内的碱基质量均值，并切除低质量碱基。（类似Trimmomatic，但是非常快）。
头尾剪裁。（上下机序列可能不稳定）
去除接头。（可以不用输入接头序列，算法可以自动识别接头序列并进行剪裁）
不匹配碱基对矫正。（对于双端测序(PE)的数据，软件能够矫正重叠区域的不匹配碱基对。如：重叠区域内，一个碱基质量很高，另一个非常低，fastp会依据质量值进行校正。
修剪尾部(3')的polyG。（修剪ployX尾，可以去掉不想要的ploy。如：mRNA-Seq 中的ployA。）
处理使用了唯一分子标识符（UMI）的数据，并将UMI转换为序列名称。
不仅生成质控和过滤结果的HTML报告（存储在fastp.html，可用-h可指定），还生成了程序可读性非常强的JSON报告（fastp.json，可用-j指定）。
为了适合并行处理，fastp将输出进行分拆。有两种模式可选，a. 指定拆分的个数; b. 指定拆分后每个文件的行数。
fastp支持gzip的输入和输出，同时支持SE和PE数据处理；支持PacBio/Nanopore的长reads数据；也支持交错输入。

3. 安装fastp

Bioconda 安装

# note: the fastp version in bioconda may be not the latest
conda install -c bioconda fastp

Linux 系统安装

# this binary was compiled on CentOS, and tested on CentOS/Ubuntu
wget http://opengene.org/fastp/fastp
chmod a+x ./fastp

源下载安装

# get source (you can also use browser to download from master or releases)
git clone https://github.com/OpenGene/fastp.git
# build
cd fastp
make
# Install
sudo make install

注：如编译时遇到(undefined reference to gzbuffer)，需要更新zlib。

4. 输入和输出

fastp支持双端(PE)数据和单端(SE)数据的输入和输出
- 对于SE数据，使用指令-i 或 --in1输入，-o 或--out1输出。
- 对于PE数据，使用指令-I 或 --in2输入，-O或--out2输出。
- 如果不指定输出文件名，而且输入名以.gz结尾，QC会执行完并压缩。
STDIN输入
- --stdin 指明使用STDIN输入。
- --interleaved_in 指定interleaved paired-end stream.
STDOUT输出 - fastp将结果传到STDOUT，以便传递到bwa和bowtie2。
- --stdout

对于PE数据，输出将与FASTQ交错，因此会有 record1-R1 -> record1-R2 -> record2-R1 -> record2-R2 -> record3-R1 -> record3-R2 ...记录。

store the unpaired reads for PE data
store the reads that fail the filters
process only part of the data
do not overwrite exiting files
split the output to multiple files for parallel processing
merge PE reads

5. 功能介绍

质量过滤 - 质量过滤功能默认开启
- -Q 或 disable_quality_filtering关闭。
- -n 或 --n_base_limit限制碱基N的数量
- -q 或 --qualified_quality_phred指定合格的质量值。（默认为 -q 15，即质量值大于等于Q15的即为合格）
- -u 或 --unqualified_percent_limit指定不合格碱基的最高百分比。（默认40%）
- -e 或--average_qual通过质量百分比过滤，如果质量值低于平均值，就被丢弃。
长度过滤 - 长度过滤默认开启
- -L 或 --disable_length_filtering关闭。
- -l 或 --length_required指定所需最小长度。（-l 30表示低于30个碱基的read会被丢弃）
- --length_limit指定最长read长度。（mall RNA sequencing可用）

结合起来就可以实现常用的discard模式，以保证所有输出的序列都一样长。

低复杂度过滤 - 默认关闭
- -y 或 --low_complexity_filter开启。

复杂度由的定义为reads中与下一个碱基不同的碱基的百分比。

# a 51-bp sequence, with 3 bases that is different from its next base

seq ='AAAATTTTTTTTTTTTTTTTTTTTTGGGGGGGGGGGGGGGGGGGGGGCCCC'

complexity = 3/(51-1) = 6%

-Y 或 --complexity_threshold 指定阈值，默认为30%。

在fastp的报告中，第一个Summary表格很清楚地显示了过滤的统计信息。

接头处理 - 默认开启，自动识别

-A 或 --disable_adapter_trimming关闭。

对于SE数据，可以用a 或 --adapter_sequence参数来输入接头。（fastp通过分析前1M的reads估计接头序列，对于SE数据，估计可能不准确，而且如果接头序列是特制的，也可能导致估计出错。）

对于PE数据，自动检测接头功能是关闭的，--detect_adapter_for_pe开启。亦可以用--adapter_sequence和 --adapter_sequence_r2指定R1, R2接头序列。
--adapter_fasta 还可以用file 的形式指定需要修剪的任何序列，序列长度大＞6bp。

>Illumina TruSeq Adapter Read 1

AGATCGGAAGAGCACACGTCTGAACTCCAGTCA

>Illumina TruSeq Adapter Read 2

AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT

>polyA

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA

修剪顺序：--adapter_sequence | --adapter_sequence_r2 | --adapter_fasta

修剪的接头序列可以在fastp报告中找到.

根据质量修剪每条read - 有三种选择，可以全选

-5或 --cut_front指定从5'开始丢弃低于质量低的，同时多N碱基也会被去除；cut_front_window_size设定滑窗大小； cut_front_mean_quality指定阈值。
-3或 --cut_tail, cut_tail_window_size, cut_tail_mean_quality.
-r或 --cut_right, cut_right_window_size, cut_right_mean_quality.（滑窗从头向尾移动，如果窗内平均质量低于阈值，丢弃滑窗内和右边所有序列。）

*all these three operations will interfere deduplication for SE data, and --cut_front or --cut_right may also interfere deduplication for PE data. *

If --cut_right is enabled, then there is no need to enable --cut_tail, since the former is more aggressive. If --cut_right is enabled together with --cut_front, --cut_front will be performed first before --cut_right to avoid dropping whole reads due to the low quality starting bases.

cut_front will interfere deduplication for both PE/SE data, and --cut_tail will interfere deduplication for SE data.

PE数据的碱基矫正 - 默认关闭

通过overlap检索，矫正后碱基质量和对照同等。

-c or --correction 开启. 可调参数overlap_len_require (default 30), overlap_diff_limit (default 5) and overlap_diff_limit_percent (default 20%).

全局裁剪 - 统一剪裁reads的头部和尾部

用于统一去掉低质量cycle。如：Illumina最后一个cycle通常质量是非常低，通过-t 1 or --trim_tail1=1去掉。

For read1 or SE data, the front/tail trimming settings are given with -f, --trim_front1 and -t, --trim_tail1.

For read2 of PE data, the front/tail trimming settings are given with -F, --trim_front2 and -T, --trim_tail2. （如read2无指定，read2裁剪参数 = read1。）

-b或 --max_len1指定read1长度, -B或 --max_len2指定read2长度。（如read2无指定，read2裁剪参数 = read1。）

1, UMI preprocessing (--umi)

2, global trimming at front (--trim_front)

3, global trimming at tail (--trim_tail)

4, quality pruning at 5' (--cut_front)

5, quality pruning by sliding window (--cut_right)

6, quality pruning at 3' (--cut_tail)

7, trim polyG (--trim_poly_g, enabled by default for NovaSeq/NextSeq data)

8, trim adapter by overlap analysis (enabled by default for PE data)

9, trim adapter by adapter sequence (--adapter_sequence, --adapter_sequence_r2. For PE data, this step is skipped if last step succeeded)

10, trim polyX (--trim_poly_x)

11, trim to max length (---max_len)

polyG裁剪 - 默认对Illumina NextSeq/NovaSeq data开启

对于两色发光法的Illumina设备（NextSeq /NovaSeq），在没有光信号情况下base calling的结果会返回G，所以在序列的尾端可能会出现较多的polyG，需要被去除。

fastp会通过机器ID自动化地识别NextSeq / NovaSeq的数据，然后进行polyG识别和裁剪。

-g or --trim_poly_g 向其他数据开启；-G or --disable_trim_poly_g关闭。

<poly_g_min_len>设定最小长度，默认值10。

polyX尾裁剪 - 默认关闭

-x or --polyX; <poly_x_min_len>, 默认值10。

同时开启时的裁剪顺序：polyG | polyA。

分子标签(UMI)处理 - 消除重复和纠正错误

-U 或 -umi 开启。

--umi_loc指定UMI所在的位置，可以是（index1、 index2、 read1、 read2、 per_index、 per_read ）中的一种；分别表示UMI是在index位置上，还是在插入片段中。如果指定了是在插入序列中，还需要使用 --umi_len 参数来指定UMI所占的碱基长度。

#初始序列

@NS500713:64:HFKJJBGXY:1:11101:1675:1101 1:N:0:TATAGCCT+GACCCCCA

AAAAAAAAGCTACTTGGAGTACCAATAATAAAGTGAGCCCACCTTCCTGGTACCCAGACATTTCAGGAGGTCGGGAAA

+

6AAAAAEEEEE/E/EA/E/AEA6EE//AEE66/AAE//EEE/E//E/AA/EEE/A/AEE/EEA//EEEEEEEE6EEAA



fastp -i R1.fq -o out.R1.fq -U --umi_loc=read1 --umi_len=8



@NS500713:64:HFKJJBGXY:1:11101:1675:1101:AAAAAAAA 1:N:0:TATAGCCT+GACCCCCA

GCTACTTGGAGTACCAATAATAAAGTGAGCCCACCTTCCTGGTACCCAGACATTTCAGGAGGTCGGGAAA

+

EEE/E/EA/E/AEA6EE//AEE66/AAE//EEE/E//E/AA/EEE/A/AEE/EEA//EEEEEEEE6EEAA

输出拆分 - 两种模式

-s or --split指定file数目。

-S or --split_by_lines指定file中的行数。

**过表达序列(Overrepresented sequence analysis) **- 默认关闭

-p or --overrepresentation_analysis开启。( count 长度为10bp, 20bp, 40bp, 100bp or (cycles - 2 )，默认值20兼顾速度与精确性)

在fastp的报告中，不仅提供overrepresented sequence的序列个数和占比，还提供了其在测序cycles中的分布情况。

PE数据reads拼接

对于PE数据，fastp通过-m/--merge选项实现拼接模式。

cheatsheet