fastp-fastq-preprocessing

fastp — FASTQ 质量控制与接头剪切工具

概述

fastp 可通过一次多线程处理，为 Illumina FASTQ 文件剪切接头、过滤低质量 reads 并生成 QC 报告。它会自动根据双端 reads 的重叠区域检测接头序列，无需手动指定接头。fastp 会校正双端重叠区域中的错配，按质量分数和长度过滤 reads，去除 polyX 尾（RNA-seq 中为 polyA），并生成交互式 HTML 和机器可读的 JSON QC 报告。fastp 的速度比 Trim Galore 和 Trimmomatic 快 3–10 倍，结果相当或更好，因此已成为 WGS、RNA-seq 和 ChIP-seq 流程中比对前的标准预处理步骤。

适用场景

• 在任何 NGS 流程（RNA-seq、WGS、WES、ChIP-seq、ATAC-seq）中，于比对前剪切 Illumina 接头和低质量碱基
• 在流程第一步生成每个样本的 QC 报告（HTML + JSON），供后续 MultiQC 汇总
• 处理双端 reads，且倾向于根据重叠区域自动检测接头，而非手动指定接头
• 去除来自 3′ 端富集方案（Smart-seq、QuantSeq）的 RNA-seq reads 中的 polyA 尾
• 按 UMI 或索引拆分 FASTQ 文件，用于解复用工作流
• 当剪切时还需要 FastQC 提供的详细逐碱基质量报告时，可使用 Trim Galore 作为替代方案
• 当需要精细控制滑动窗口剪切步骤时，可使用 Trimmomatic 作为替代方案

前置条件

• 软件：fastp（conda 或预编译二进制文件）
• 输入：原始 Illumina FASTQ 文件（单端或双端，.fastq 或 .fastq.gz）

安装前检查：该工具可能已在当前环境中可用（例如，在 pixi / conda 环境中）。先运行 command -v fastp；如果返回路径，则跳过以下安装命令。在 pixi 项目中运行时，请使用 pixi run fastp，而不是直接运行 fastp。

# Install with conda
conda install -c bioconda fastp

# Or download pre-compiled binary (Linux)
wget https://github.com/OpenGene/fastp/releases/download/v0.24.0/fastp
chmod +x fastp
./fastp --version
# fastp 0.24.0

# Verify
fastp --version

准备阶段访谈

编写任何分析代码之前，先与用户确认以下事项。

decisions:
  - id: D1
    param: -i / -I
    kind: required
    source: data
    ask: "Were these reads sequenced from single-end or paired-ends of each fragment?"
    default: "read off the file inventory - paired when R1/R2 mates are present"

  - id: D2
    param: library_prep
    kind: required
    source: user
    ask: "Which library prep produced these reads - 5'-, 3'-, small RNA, lnRNA or amplicon?"
    default: null

  - id: D3
    param: --trim_poly_x
    kind: optional
    source: user
    depends_on: [D2]
    ask: "Trim polyA/polyT tails left by the capture step?"
    default: "on for 3'-enriched RNA-seq, off otherwise"

  - id: D4
    param: --adapter_sequence / --detect_adapter_for_pe
    kind: required
    source: user
    depends_on: [D1]
    ask: "Auto-detect the adapter, or supply the kit's sequence explicitly?"
    default: "auto-detect"

  - id: D5
    param: --length_required
    kind: required
    source: user
    depends_on: [D2]
    ask: "Reads shorter than this after trimming are discarded - how short is too short to map?"
    default: "15 (raise to ~36 for standard RNA-seq; keep low for small RNA)"

  - id: D6
    param: --qualified_quality_phred
    kind: optional
    source: user
    ask: "How confident must a base call be before it counts as good quality?"
    default: "15 (20 = 1% error)"

  - id: D7
    param: --correction
    kind: optional
    source: user
    depends_on: [D1]
    ask: "Correct mismatched bases where the two mates overlap?"
    default: "off"
    skip_if: "single-end input - there is no overlap to compare"

  - id: D8
    param: --low_complexity_filter
    kind: optional_conditional
    source: data
    ask: "Drop low-complexity reads (homopolymer and short-repeat runs)?"
    default: "off"

  - id: D9
    param: --thread, --split
    kind: never_ask
    source: data
    reason: "Parallelism and output sharding; affect runtime and file layout, not the reads kept"
    default: "8 threads, no split"

直接询问 D2，而不是自行推测，因为样本制备方式会影响后续两个答案：它决定 polyA 剪切是否正确（D3；3′ 富集文库末端带有捕获尾，该尾部不属于生物序列），也决定长度下限应如何设置（D5——适用于 RNA-seq 的 36 bp 最小长度，会丢弃小 RNA 测序中整条 read）。D4 和 D7 取决于 D1，因为两者都涉及如何处理一对 mate reads。

快速开始

# Paired-end adapter trimming with QC report
fastp \
    -i sample_R1.fastq.gz \
    -I sample_R2.fastq.gz \
    -o sample_R1.trimmed.fastq.gz \
    -O sample_R2.trimmed.fastq.gz \
    -h sample_qc.html \
    -j sample_qc.json \
    --thread 8

echo "Trimmed reads in: sample_R1.trimmed.fastq.gz"

工作流

步骤 1：单端接头剪切

在单端 FASTQ 上运行 fastp，自动检测接头。

# Single-end with auto adapter detection
fastp \
    -i sample.fastq.gz \
    -o sample.trimmed.fastq.gz \
    -h sample_qc.html \
    -j sample_qc.json \
    --thread 8 \
    --qualified_quality_phred 20 \
    --length_required 36

echo "Input reads:   $(zcat sample.fastq.gz | wc -l | awk '{print $1/4}')"
echo "Output reads:  $(zcat sample.trimmed.fastq.gz | wc -l | awk '{print $1/4}')"

步骤 2：双端接头剪切

处理双端 FASTQ 文件，根据重叠区域检测接头并进行校正。

# Paired-end with overlap-based adapter auto-detection
fastp \
    -i sample_R1.fastq.gz \
    -I sample_R2.fastq.gz \
    -o sample_R1.trimmed.fastq.gz \
    -O sample_R2.trimmed.fastq.gz \
    -h sample_qc.html \
    -j sample_qc.json \
    --thread 8 \
    --correction \
    --detect_adapter_for_pe \
    --qualified_quality_phred 20 \
    --length_required 36

# Specify adapters explicitly (if auto-detection fails)
# fastp -i R1.fq.gz -I R2.fq.gz \
#   --adapter_sequence AGATCGGAAGAGCACACGTCTGAACTCCAGTCA \
#   --adapter_sequence_r2 AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT \
#   -o R1.out.fq.gz -O R2.out.fq.gz

步骤 3：质量过滤与 read 长度剪切

设置质量和长度阈值，以进行更严格或更宽松的过滤。

# Strict quality filtering (e.g., for variant calling)
fastp \
    -i sample_R1.fastq.gz \
    -I sample_R2.fastq.gz \
    -o sample_R1.filtered.fastq.gz \
    -O sample_R2.filtered.fastq.gz \
    -h sample_qc.html \
    -j sample_qc.json \
    --thread 8 \
    --qualified_quality_phred 25 \
    --unqualified_percent_limit 20 \
    --length_required 50 \
    --max_len1 150 \
    --max_len2 150 \
    --low_complexity_filter \
    --complexity_threshold 30

echo "Filtering complete. Check sample_qc.html for pass/fail rates."

步骤 4：去除 RNA-seq polyA 尾

在比对前，去除 3′ 富集 RNA-seq 方案中的 polyA 尾。

# Remove polyA tails (QuantSeq 3′ mRNA-seq)
fastp \
    -i quantseq_R1.fastq.gz \
    -o quantseq_R1.trimmed.fastq.gz \
    -h quantseq_qc.html \
    -j quantseq_qc.json \
    --thread 8 \
    --trim_poly_x \
    --poly_x_min_len 10 \
    --qualified_quality_phred 20 \
    --length_required 25

# For Smart-seq2 paired-end with polyA
fastp \
    -i smartseq_R1.fastq.gz \
    -I smartseq_R2.fastq.gz \
    -o smartseq_R1.trimmed.fastq.gz \
    -O smartseq_R2.trimmed.fastq.gz \
    --trim_poly_x --poly_x_min_len 10 \
    --thread 8 \
    -h smartseq_qc.html -j smartseq_qc.json

步骤 5：解析 QC 报告 JSON 以监控流程

从 fastp 的 JSON 输出中提取关键 QC 指标，用于自动化质量门控。

import json
from pathlib import Path

def parse_fastp_json(json_path: str) -> dict:
    with open(json_path) as f:
        data = json.load(f)
    
    before = data["summary"]["before_filtering"]
    after = data["summary"]["after_filtering"]
    
    return {
        "total_reads_in":  before["total_reads"],
        "total_reads_out": after["total_reads"],
        "pct_passed":      after["total_reads"] / before["total_reads"] * 100,
        "q30_rate_before": before["q30_rate"] * 100,
        "q30_rate_after":  after["q30_rate"] * 100,
        "mean_len_before": before["read1_mean_length"],
        "mean_len_after":  after["read1_mean_length"],
        "adapter_trimmed": data["filtering_result"]["adapter_trimmed"],
    }

metrics = parse_fastp_json("sample_qc.json")
for key, val in metrics.items():
    print(f"{key:25s}: {val:.1f}" if isinstance(val, float) else f"{key:25s}: {val:,}")

# Quality gate: fail if < 70% reads pass filter
if metrics["pct_passed"] < 70:
    print("WARNING: Low pass rate — check raw data quality")

步骤 6：批量预处理流程

依次处理多个样本，并汇总每个样本的 QC 结果。

#!/bin/bash
# Batch paired-end preprocessing for multiple samples
SAMPLES=(ctrl_1 ctrl_2 treat_1 treat_2)
DATA="data"
OUT="trimmed"
QC="qc/fastp"
THREADS=8

mkdir -p "$OUT" "$QC"

for sample in "${SAMPLES[@]}"; do
    echo "=== Processing $sample ==="
    fastp \
        -i "$DATA/${sample}_R1.fastq.gz" \
        -I "$DATA/${sample}_R2.fastq.gz" \
        -o "$OUT/${sample}_R1.fastq.gz" \
        -O "$OUT/${sample}_R2.fastq.gz" \
        -h "$QC/${sample}.html" \
        -j "$QC/${sample}.json" \
        --thread $THREADS \
        --correction \
        --detect_adapter_for_pe \
        --qualified_quality_phred 20 \
        --length_required 36 \
        2>&1 | grep -E "Read[12]|Filtering|Adapter|passed"
done

# Aggregate QC metrics
python3 - << 'EOF'
import json, pandas as pd
from pathlib import Path

rows = []
for jf in sorted(Path("qc/fastp").glob("*.json")):
    with open(jf) as f: data = json.load(f)
    after = data["summary"]["after_filtering"]
    before = data["summary"]["before_filtering"]
    rows.append({
        "sample": jf.stem,
        "reads_in": before["total_reads"],
        "reads_out": after["total_reads"],
        "pct_passed": round(after["total_reads"]/before["total_reads"]*100, 1),
        "q30_after": round(after["q30_rate"]*100, 1),
    })
df = pd.DataFrame(rows)
print(df.to_string(index=False))
df.to_csv("fastp_summary.tsv", sep="\t", index=False)
EOF

# Run MultiQC to aggregate all fastp JSON reports
multiqc qc/fastp/ -o qc/ -n fastp_multiqc_report

关键参数

参数	默认值	范围/选项	作用
-i / -I	必需	文件路径	输入 FASTQ（双端测序时为 R1 和 R2）
-o / -O	必需	文件路径	输出修剪后的 FASTQ（R1 和 R2）
-h / -j	—	文件路径	HTML 和 JSON QC 报告输出路径
--thread	3	1–16	CPU 线程数；8 是较好的平衡
--qualified_quality_phred	15	0–40	最低碱基质量（Phred）；20 = 1% 错误率
--length_required	15	1–1000	修剪后的最短读长；丢弃更短的 reads
--correction	关闭	标志	校正 PE 重叠区域中的错配
--detect_adapter_for_pe	关闭	标志	为 PE 数据启用基于重叠的接头自动检测
--adapter_sequence	自动	字符串	显式指定 R1 接头；覆盖自动检测
--trim_poly_x	关闭	标志	修剪 polyX（polyA/polyT）尾部；用于 3′ 富集 RNA-seq
--low_complexity_filter	关闭	标志	过滤低复杂度 reads（默认复杂度低于 30%）
--split	关闭	整数	将输出按方向拆分为 N 个文件（用于并行处理）

常用配方

配方 1：将 fastp 集成到 Snakemake 工作流中

# Snakefile — fastp trimming rule
configfile: "config.yaml"
SAMPLES = config["samples"]

rule fastp_pe:
    input:
        r1 = "data/{sample}_R1.fastq.gz",
        r2 = "data/{sample}_R2.fastq.gz"
    output:
        r1 = "trimmed/{sample}_R1.fastq.gz",
        r2 = "trimmed/{sample}_R2.fastq.gz",
        html = "qc/{sample}_fastp.html",
        json = "qc/{sample}_fastp.json"
    threads: 8
    shell:
        """
        fastp -i {input.r1} -I {input.r2} \
              -o {output.r1} -O {output.r2} \
              -h {output.html} -j {output.json} \
              --thread {threads} \
              --correction --detect_adapter_for_pe \
              --qualified_quality_phred 20 \
              --length_required 36
        """

配方 2：使用 Python 汇总 fastp JSON 报告

import json
import pandas as pd
from pathlib import Path

qc_dir = Path("qc/fastp")
records = []

for jf in sorted(qc_dir.glob("*.json")):
    with open(jf) as f:
        d = json.load(f)
    b = d["summary"]["before_filtering"]
    a = d["summary"]["after_filtering"]
    records.append({
        "sample": jf.stem.replace("_fastp", ""),
        "reads_in_M": b["total_reads"] / 1e6,
        "reads_out_M": a["total_reads"] / 1e6,
        "pct_passed": a["total_reads"] / b["total_reads"] * 100,
        "q30_pct": a["q30_rate"] * 100,
        "mean_len_bp": a["read1_mean_length"],
        "adapter_pct": d["filtering_result"]["adapter_trimmed"] / b["total_reads"] * 100,
    })

df = pd.DataFrame(records).round(2)
print(df.to_string(index=False))

# Flag low-quality samples
low_q = df[df["pct_passed"] < 80]
if not low_q.empty:
    print(f"\nSamples with < 80% reads passing: {list(low_q['sample'])}")

预期输出

输出	格式	描述
*_R1.trimmed.fastq.gz	FASTQ.gz	修剪后的 R1 reads（已移除接头和低质量碱基）
*_R2.trimmed.fastq.gz	FASTQ.gz	修剪后的 R2 reads（仅限双端测序）
*.html	HTML	交互式 QC 报告，包含逐碱基质量、GC 含量和接头图表
*.json	JSON	机器可读的 QC 指标，用于自动化处理和 MultiQC 解析
fastp.log	文本	stderr 摘要，包含通过/失败的 read 数量和过滤统计信息

故障排除

问题	原因	解决方案
SE 模式下未检测到接头	SE reads 需要显式指定接头或使用 --adapter_sequence	--detect_adapter_for_pe 仅用于 PE；为 SE 指定接头：--adapter_sequence AGATCGGAAGAGC
接头含量过高（> 50%）	插入片段较短（small RNA、miRNA）或文库制备效果不佳	检查文库方案；使用 --overlap_len_require 10 调整重叠灵敏度
过滤掉的 reads 过多（通过率 < 60%）	质量阈值过严或测序运行质量较差	将 --qualified_quality_phred 放宽至 15；将 --length_required 降至 25
JSON 输出缺少字段	fastp 版本过旧	升级：conda update fastp，或从 GitHub 下载最新二进制文件
MultiQC 未解析 fastp JSON	JSON 文件不在扫描目录中	运行 multiqc qc/ 而不是 multiqc .；使用 ls qc/*.json 确认 JSON 文件存在
输出 FASTQ 为空	所有 reads 均被过滤（输入错误或阈值过严）	使用 zcat sample.fq.gz | head -8 检查输入 FASTQ；先在不启用 --low_complexity_filter 的情况下运行
大文件处理速度慢	线程数较少	将 --thread 增至 8–12；确保输入文件位于高速存储设备（SSD）上
未移除 polyA	未设置 --trim_poly_x	对 3′ 富集方案添加 --trim_poly_x --poly_x_min_len 10

参考资料

• fastp GitHub：OpenGene/fastp — 源代码、更新日志和使用指南
• Chen S 等（2018），“fastp：超高速一体化 FASTQ 预处理器” — Bioinformatics 34(17):i884-i890。DOI:10.1093/bioinformatics/bty560
• fastp 文档 — 完整参数参考和配方
• MultiQC fastp 模块 — 汇总多个样本的 fastp JSON 报告