第三章 原始数据质量控制：Trimmomatic和FastQC的参数调整与过滤策略

拿到测序数据后，第一件事绝对不是直接跑物种鉴定——你得先看看数据质量怎么样。说白了，垃圾进垃圾出，这道理放在哪都一样。我最早做宏基因组项目时，就因为这一步没认真做，后面折腾了整整一周才发现是接头污染导致的假阳性。嗯，从那以后我就养成了“先质检再处理”的死规矩。

核心思想：FastQC 是照妖镜，Trimmomatic 是手术刀。先用 FastQC 看到底哪里脏，再拿 Trimmomatic 精确修剪。
原始数据质量控制流程 原始 FASTQ FastQC 评估 判断？ 直接分析 Trimmomatic 再次 FastQC

3.1 FastQC 指标解读：哪些指标值得你重点关注？

跑 FastQC 很简单，一行命令的事：

fastqc sample_R1.fastq sample_R2.fastq -o ./qc_report -t 4

出来会得到一个 HTML 报告。我个人习惯首先看三个最要命的指标：

  • Per base sequence quality —— 每个碱基位置的质量分布。如果整个读段尾部大面积飘红（质量低于 Q20），就必须考虑切尾部。
  • Adapter Content —— 接头含量。宏基因组双端测序经常因为插入片段太短导致两头残留接头，这个不干掉，后面 mapping 时会引入大量非真实比对。
  • Overrepresented sequences —— 过度代表序列。有的是 rRNA 污染，有的是真实的保守区域。我遇到过一种情况：全是同一条重复序列，后来发现是文库构建时扩增偏倚。这个得仔细排除。
小技巧：FastQC 的“Per tile sequence quality”可以帮你发现 flowcell 上某个特定 tile 的质量异常。如果某个 tile 整块不行，你得考虑要不要舍弃这个泳道的数据。我做过一个项目，就是因为忽略了 tile 异常，导致拼接结果里出现大量嵌合体。

3.2 Trimmomatic 参数配置：怎么切才不伤数据？

Trimmomatic 的参数看着多，其实大部分场景只用调整 5 个核心参数。下面这张表整理了我的常用取值和推荐理由：

参数 作用 我常用的值 使用场景
ILLUMINACLIP去除接头和引物2:30:10:8:true双端通用，接头库用 TruSeq3-PE-2.fa
LEADING从 5' 端切除低质量碱基3几乎必开，质量<3 的碱基直接扔掉
TRAILING从 3' 端切除低质量碱基3同上，但建议结合滑动窗口一起用
SLIDINGWINDOW滑动窗口平均质量4:15核心过滤，窗口4bp，平均质量<15 则切除
MINLEN最短保留长度50过短读段对物种分类贡献极低，直接舍弃

具体到命令行，长这样（我用的是双端模式）：

trimmomatic PE \
  sample_R1.fastq.gz sample_R2.fastq.gz \
  sample_R1_trimmed.fastq.gz sample_R1_unpaired.fastq.gz \
  sample_R2_trimmed.fastq.gz sample_R2_unpaired.fastq.gz \
  ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10:8:true \
  LEADING:3 \
  TRAILING:3 \
  SLIDINGWINDOW:4:15 \
  MINLEN:50 \
  -threads 8

这里面的 8:true 值得多说一句。很多人不知道这个参数是干嘛的——它控制了在允许的错配数下，只匹配接头的最短长度。我一般设到 8 个碱基，既能识别残留接头，又不会把真正的插入序列误剪掉。

3.3 参数调整的避坑指南

避坑 1：不要盲目用全局默认值。 我曾经接手一个土壤宏基因组项目，同事直接用了网上抄的 Trimmomatic 参数，SLIDINGWINDOW 设成了 4:20。结果呢？读段长度从 150 bp 被削到 60 bp，而且大量本来 Q30 的碱基也被无情切掉。后来我改成 4:15，数据使用率从 35% 直接升到 72%。你想想看，窗口平均质量 20 对于低质量样本来说太苛刻了，宏基因组本身 GC 含量波动大，部分区域 Q 值偏低是正常的。

避坑 2：ILLUMINACLIP 参数中的 palindromeClipThreshold 和 simpleClipThreshold 怎么调？ 双端数据中接头常常以“回文”形式出现（两条 read 的 3' 端互相配对）。我一般把回文匹配阈值设 30（2:30:10:8:true），简单匹配阈值设 10。如果数据中接头含量特别高（比如超过 5%），我会把简单匹配阈值降到 8，让算法更积极地剪掉短接头。

避坑 3：注意 paired-end 模式下未配对输出的处理。 双端过滤后，可能有一端被切掉而另一端保留，这时 Trimmomatic 会把配对失败的读段分别写到 *_unpaired.fastq.gz 里。我强烈建议你把这些 unpaired 保留下来——它们在后续能用于 k-mer 分析（比如 MetaPhlAn 只要求单端也行），别直接删掉。

警告：不要一次地把所有低质量碱基都砍光。我见过有人用 LEADING:20 + TRAILING:20 + SLIDINGWINDOW:1:30，结果一万条 read 只剩下 800 条能用。质量控制的目的是去除被污染的部分，而不是追求100%高质量——我们做的是丰度检测，不是单碱基测序。

3.4 如何根据 FastQC 报告决定参数？实战案例

假设你拿到一份 FastQC 报告，发现如下异常：

  • Per base quality 在 100-150 bp 位置掉到 Q20 以下
  • Adapter Content 在 120 bp 处开始抬头
  • 没有 rRNA 或过度代表序列

这种情况最典型：插入片段约 300 bp，但读长 150 bp，两端覆盖不到读段中间，导致 3' 端质量下降且接头未完全测通。我的调整方案是：

trimmomatic PE \
  ... \
  LEADING:3 \
  TRAILING:3 \
  SLIDINGWINDOW:4:15 \
  MINLEN:90 \
  # 不对！MINLEN 设得太高会丢掉很多全长片段
  # 改成 MINLEN:75

实际操作中我会先把 MINLEN 设为 36（默认最小值），让数据先跑一轮，再用 FastQC 检查修剪后的质量。如果 3' 端仍然有残留接头，就把 SLIDINGWINDOW 的窗口缩小到 4:13，或者干脆加一个 CROP 参数把 3' 端固定切掉 20 bp。我的经验是：能调滑动窗口解决的，就别动固定切除。 因为固定切除对长度变化敏感，而宏基因组样本的插入片段长度本就分布广。

3.5 多线程与输出文件管理

Trimmomatic 支持多线程，但注意它的 -threads 参数只在线程内部并行处理压缩和解压缩。如果你在 HPC 上跑，建议把线程数设成核心数减 1，给系统留点余地。另外，输出文件命名最好带上版本号或参数摘要，方便后期回溯。我习惯这么写：

# 例：sample_R1_TRIMMED_L3_T3_W4x15_MIN75.fastq.gz
# 这样看到文件名就知道参数：LEADING 3, TRAILING 3, SLIDINGWINDOW 4:15, MINLEN 75

3.6 综合参数配置模板（适用于 150 bp PE 宏基因组）

下面是我个人最常用的一套模板，经过 20+ 个项目的验证，通用性较好：

trimmomatic PE -threads 8 -phred33 \
  in_R1.fastq.gz in_R2.fastq.gz \
  out_R1_paired.fastq.gz out_R1_unpaired.fastq.gz \
  out_R2_paired.fastq.gz out_R2_unpaired.fastq.gz \
  ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10:8:true \
  LEADING:3 \
  TRAILING:3 \
  SLIDINGWINDOW:4:15 \
  MINLEN:50

如果是低质量样本（比如环境水样、冻存组织），我会把 SLIDINGWINDOW 降到 4:12，同时把 MINLEN 提到 60，因为短读段在这种样本中几乎都是引物二聚体或噪声。

我的黄金法则：修剪前先跑一次 FastQC，修剪后再跑一次 FastQC。两次对比，一看每碱基质量箱线图有没有整体抬升，二看接头含量是不是趋近于零。如果这两点都做到了，数据质量就算过关了。
← 上一章 📖 返回目录 下一章 →