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

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. 输入和输出

  1. fastp支持双端(PE)数据和单端(SE)数据的输入和输出

    • 对于SE数据,使用指令-i 或 --in1输入,-o 或--out1输出。

    • 对于PE数据,使用指令-I 或 --in2输入,-O或--out2输出。

    • 如果不指定输出文件名,而且输入名以.gz结尾,QC会执行完并压缩。

  2. STDIN输入

    • --stdin 指明使用STDIN输入。

    • --interleaved_in 指定interleaved paired-end stream.

  3. STDOUT输出 - fastp将结果传到STDOUT,以便传递到bwabowtie2

    • --stdout

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

  1. store the unpaired reads for PE data

  2. store the reads that fail the filters

  3. process only part of the data

  4. do not overwrite exiting files

  5. split the output to multiple files for parallel processing

  6. merge PE reads

5. 功能介绍

  1. 质量过滤 - 质量过滤功能默认开启

    • -Q 或 disable_quality_filtering关闭。

    • -n 或 --n_base_limit限制碱基N的数量

    • -q 或 --qualified_quality_phred指定合格的质量值。(默认为 -q 15,即质量值大于等于Q15的即为合格)

    • -u 或 --unqualified_percent_limit指定不合格碱基的最高百分比。(默认40%)

    • -e 或--average_qual通过质量百分比过滤,如果质量值低于平均值,就被丢弃。

  2. 长度过滤 - 长度过滤默认开启

    • -L 或 --disable_length_filtering关闭。

    • -l 或 --length_required指定所需最小长度。(-l 30表示低于30个碱基的read会被丢弃)

    • --length_limit指定最长read长度。(mall RNA sequencing可用)

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

  1. 低复杂度过滤 - 默认关闭

    • -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表格很清楚地显示了过滤的统计信息。

  1. 接头处理 - 默认开启,自动识别
  • -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报告中找到.

  1. 根据质量修剪每条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.

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

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

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

  1. 全局裁剪 - 统一剪裁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)

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

  1. polyX尾裁剪 - 默认关闭

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

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

  1. 分子标签(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
  1. 输出拆分 - 两种模式

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

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

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

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

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

  1. PE数据reads拼接

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

cheatsheet

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,185评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,445评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,684评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,564评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,681评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,874评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,025评论 3 408
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,761评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,217评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,545评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,694评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,351评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,988评论 3 315
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,778评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,007评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,427评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,580评论 2 349

推荐阅读更多精彩内容