4. 测序数据质量控制：FastQC、Trimmomatic 等工具使用

体细胞突变分析的起点，说实话就是原始测序数据。你从测序仪拿到的 FASTQ 文件，里面藏着各种“坑”——接头残留、低质量碱基、甚至还有测序仪自身的系统误差。我最早做这个的时候，觉得直接往下跑比对就行，结果突变检测结果乱七八糟，后来才发现是数据质量根本没做好。

这章我们就聊聊怎么用 FastQC、Trimmomatic 这些工具给原始数据“体检”和“清洗”。这块做扎实了，后面的分析才能站得住脚。

4.1 原始数据长什么样

先看一眼 FASTQ 文件的格式。每条序列四行：

@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

第一行以 @ 开头，是序列标识符。第二行就是碱基序列。第三行是 + 分隔。第四行是质量值，每个字符对应一个碱基的 Phred 分数。

Phred 分数我不多讲公式，你只要记住：分数越高，碱基越可靠。Q30 意味着错误率千分之一，Q20 是百分之一。临床上做突变分析，我通常要求至少 Q30 以上的碱基占比超过 80%。

4.2 FastQC：快速摸清数据底细

FastQC 是第一个上手的工具。它能给你一个 HTML 报告，包含十多项质量指标。我个人习惯在拿到任何一批新数据后，先跑一遍 FastQC，看看有没有明显的异常。

4.2.1 基础用法

fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o ./qc_reports/ -t 4

-o 指定输出目录，-t 是线程数。跑完后每个样本会生成一个 .html 和 .zip。

4.2.2 报告里重点看什么

FastQC 报告有十几个模块，但我每次只盯着这四个看：

模块 红灯代表什么 我的处理方法
Per base sequence quality 某些位置中位数低于 Q20 Trimmomatic 裁掉低质量尾巴
Per sequence GC content 出现双峰或偏离理论分布 怀疑有污染或 PCR 偏差
Adapter Content 接头比例高 必须去除，否则比对会乱
Overrepresented sequences 某条序列出现异常多次 可能是接头二聚体或重复扩增

重点关注：Per base sequence quality 图中如果出现红色警告（背景黄或红），那后面的 Trimmomatic 参数就要设严一点。我以前有一批血样，跑完 FastQC 发现 3' 端在 100bp 后质量直线下降，裁掉后面 20bp 后变异检测准确率明显提升。

4.2.3 批量报告汇总

样本多了一个个看 HTML 太累。我习惯用 MultiQC 把所有报告汇总成一个页面：

multiqc ./qc_reports/ -o ./multiqc_report/

它自动识别 FastQC 输出，生成一个漂亮的总览表。哪个样本有问题，一眼就能看出来。

4.3 Trimmomatic：精细清洗数据

FastQC 告诉你有问题，Trimmomatic 就是来解决问题的。它支持双端数据，能干三件事：去接头、滑窗裁剪、低质量末端切除。

4.3.1 安装与基本命令

Trimmomatic 是 Java 写的，直接下载 jar 包就行。我用得最多的命令是这样的：

java -jar trimmomatic-0.39.jar PE \
  sample_R1.fastq.gz sample_R2.fastq.gz \
  sample_R1_trim.fastq.gz sample_R1_unpaired.fastq.gz \
  sample_R2_trim.fastq.gz sample_R2_unpaired.fastq.gz \
  ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
  LEADING:3 TRAILING:3 \
  SLIDINGWINDOW:4:15 \
  MINLEN:36

参数看着多，我一个个解释：

  • PE 表示双端，SE 用于单端。
  • ILLUMINACLIP：接头文件、允许的最大错配数（2）、比对分值阈值（30）、最优匹配长度（10）。
  • LEADING:3：从 5' 端切掉质量低于 3 的碱基。
  • TRAILING:3：从 3' 端切掉质量低于 3 的碱基。
  • SLIDINGWINDOW:4:15：滑窗大小 4 bp，窗口平均质量低于 15 则切除。
  • MINLEN:36：保留的最短序列长度，太短了比对不到参考基因组。

小技巧：接头的文件可以从 Trimmomatic 的 adapter 目录里选，也可以自己准备。如果你用的是不同建库试剂盒，最好确认一下接头序列。我吃过一次亏——用的试剂盒接头和 TruSeq 不一致，跑了 Trimmomatic 也没去掉，后来用 cutadapt 才解决。

4.3.2 参数怎么定

参数不是一成不变的。我一般分两步走：

  1. 先跑 FastQC 看质量分布。如果碱基质量从第 100bp 开始掉，就把 SLIDINGWINDOW 的阈值调高到 20；如果前 10bp 有低质量，适当调高 LEADING。
  2. 跑一轮 Trimmomatic 后再做 FastQC，对比前后变化。直到所有模块都绿了再往下走。

举个例子：我处理过一个血样，原始数据 R2 端 3' 端质量惨不忍睹。Trimmomatic 我用了 SLIDINGWINDOW:4:20 TRAILING:20，截掉将近 30% 的碱基。看着心疼，但剩下的数据比对率和突变检测一致性都很好。

4.4 其他常用的质控工具

除了 Trimmomatic，还有几个工具我也经常用：

工具 特点 适用场景
cutadapt 聚焦接头去除，灵活支持多类型接头 复杂接头结构、MGI/华大平台数据
fastp 一站式质控，自动检测接头，速度快 全自动批量处理，不需要手动指定接头
BBTools (bbduk) 除接头、去污染、质量过滤，内存高效 超大数据集，需要精细控制每一步

我个人偏爱 Trimmomatic 和 cutadapt 搭配。cutadapt 去接头更灵活，Trimmomatic 做滑窗裁剪和长度过滤效率高。

4.5 质控之后的检查清单

数据清洗完别急着跑下一步。我给自己列了个核对清单：

  • 跑一遍 FastQC 确认所有指标绿色或黄色（没有红色）。
  • 检查总 reads 数损失：如果丢了一多半，要反思参数是不是太严了。
  • 双端数据检查配对数：两个输出文件中的 reads 数应该一致（grep -c "^@" *_trim.fastq.gz）。
  • 用初步比对（如 BWA-MEM）看一眼比对率是否达到预期（一般健康样本 ≥ 95%）。

注意：不要以为质控跑完就没问题了。我曾经遇到一批冷冻样本，Trimmomatic 后的 FastQC 全绿，结果比对时发现 GC 含量异常。后来发现是样本降解严重，大量片段长度低于 50bp，MINLEN 设成 36 根本拦不住。所以质控后还要结合比对结果综合判断。

4.6 本章知识体系总览

为了让你对整个质控流程有个整体认识，我画了张图。从 FASTQ 输入，到 FastQC 诊断，再到 Trimmomatic 清洗，最后用 FastQC 验证，形成一个闭环。

测序数据质控流程总览 原始 FASTQ FastQC 质量诊断 接头残留 → 去除接头 低质量尾部 → 滑窗裁剪 过短序列 → 长度过滤 Trimmomatic 清洗 清洗后 FASTQ 再次 FastQC 验证 参数迭代 根据 FastQC 结果调整参数

这张图的核心就是：不依赖一次性的质控。FastQC 发现问题 → Trimmomatic 清洗 → 再 FastQC 确认。循环迭代直到数据干净。我建议你在项目里把这个流程写成脚本，每次换参数改改变量就行。

4.7 避坑总结

最后分享几个我用真金白银换来的教训：

  • 双端数据别用错模式：PE 命令要同时给 R1 和 R2，不要只处理一条。
  • 接头文件版本要对：不同版本的 TruSeq 接头序列可能不同。我习惯自己从试剂盒说明书里提取序列，手动建一个 FASTA。
  • 保存中间文件：Trimmomatic 输出四个文件——两个配对的清理后文件，两个未配对的。未配对的下游可能也有用，别直接删。
  • 不要过度修剪：如果一条 read 被裁掉 90% 以上，还不如直接丢弃。一般来说，我保留 reads 损失在 30% 以内。
← 上一章 📖 返回目录 下一章 →