3. 数据获取与质量控制

拿到原始测序数据，第一件事不是急着往下跑。我见过太多人从 SRA 或测序中心下载完文件就直接丢给 BWA，结果跑出来一堆垃圾。说到底，GATK 管线下游再花哨的算法，也弥补不了上游的数据质量问题。这一章我们就把 FASTQ 格式和 FastQC 质量评估吃透。

3.1 原始测序数据格式：FASTQ

FASTQ 其实是 FASTA 的升级版。FASTA 只存序列，FASTQ 额外存了每个碱基的质量得分。我早年做转录组时吃过亏：以为格式都一样，结果把 FASTA 当 FASTQ 读，管线直接崩了。嗯，这两者不能混用。

一条标准的 FASTQ 记录包含四行：

@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
CTGTAGCCATCCTTATCGTTCGTGGCTGAAGAGCTGGGCTTGGGACCCCAGCTCTGCGCTCTTCCTTTCTGATCTCT
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

我来拆解一下：

  • 第1行：以 @ 开头，后面是测序仪 ID、流道位置、坐标等信息。说白了就是这条 reads 的身份证。
  • 第2行：碱基序列，A/T/C/G/N。N 代表这个位置测序仪没读出来。
  • 第3行：以 + 开头，有时会重复第1行的 ID，有时直接是 +。
  • 第4行：质量得分，字符个数必须等于第2行。每个字符代表一个碱基的测序置信度。
⚠️ 常见坑：双端测序的 R1 和 R2，两条 reads 的名字必须完全一致（除末尾的 1/2 标记）。我见过有人把不同样本的文件混在一起，结果 bwa 报错说 reads 不配对。检查这条能省很多调试时间。

3.2 质量得分体系：Phred 与 ASCII 编码

为什么质量得分用字符表示？因为历史原因。早期存储空间紧张，用一个字节（字符）存一个质量值比用整数更省。质量得分的计算公式是：

Q = -10 * log10(P)

P：该碱基测错的概率
Q：Phred 质量分数

举个例子：

  • Q30：错误概率 0.001（千分之一），正确率 99.9%
  • Q20：错误概率 0.01（百分之一），正确率 99%

核心指标：我个人习惯用 Q30 占比来判定文库质量。低于 80% 的 Q30 占比，我会谨慎怀疑是文库降解或测序跑偏了。

质量值得分映射成 ASCII 字符时，有两种主流标准：

版本 偏移量 适用场景 最小 ASCII
Sanger / Illumina 1.8+ Phred+33 目前绝大多数据 ! (Q0)
Illumina 1.3–1.7 Phred+64 老数据（2011年前） B (Q2)

你想想看，如果老数据用 Phred+64，新分析软件默认用 Phred+33，那质量分全读偏了。FastQC 会自动检测编码方式，但如果你手动写脚本解析 FASTQ，一定要先确认第一行 @ 后面的版本号。

实用技巧：用 head -n 40 看质量得分行，如果出现大量 ! 或 # 等低 ASCII 字符，通常是 Phred+33；如果最低字符是 B 或 C，那多半是 Phred+64。FastQC 会自动区分，但我建议还是手动扫一眼。

3.3 质量评估工具：FastQC

FastQC 是我用过最顺手的质控工具。它不需要写复杂的参数，跑一个命令就能出一份 HTML 报告。我记得有一次批量处理 48 个样本，全丢给 FastQC 并行跑，半小时不到所有报告就出来了。

3.3.1 安装与运行

FastQC 可以用 conda 安装，非常方便：

conda install -c bioconda fastqc

运行也很简单：

fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o ./qc_reports/ -t 4
  • -o ./qc_reports/：指定输出目录
  • -t 4：用 4 个线程并行处理

输出两个文件：sample_R1_fastqc.html（可视化报告）和 sample_R1_fastqc.zip（原始数据）。我习惯留着 zip，后面批量汇总时会用到 MultiQC。

3.3.2 报告的 11 个核心模块

FastQC 报告里有 11 个模块，但并非每个都需要逐行看。我挑几个最关键的讲：

模块 重点关注 我的经验阈值
Per base sequence quality 每个位置的质量箱线图 前 10 bp 质量偏低常见，末尾低于 Q20 需考虑截断
Per tile sequence quality 流动槽上每个 tile 的质量热图 如果整列 tile 全红，说明该 lane 有污染
Per base sequence content 四个碱基的比例曲线 前 10 bp 有 GC 抖动正常（随机引物），后期应平稳
Adapter Content 接头污染比例 超过 5% 的 reads 含接头，必须切除
Overrepresented sequences 高频出现的异常序列 若匹配已知接头或 rRNA，立即处理
⚠️ 避坑：Per tile sequence quality
有的实验室用老测序仪，偶尔出现单 tile 质量下降。我曾经遇到一个项目，16 个样本里有 3 个的 tile 热图出现一整列橘红色。如果直接跳过不看，下游变异检测会在这个区域引入系统性假阳性。

3.4 从 FASTQ 到 FastQC 的工作流程

我一般把这一阶段的流程画成下图，方便团队理解：

数据获取与质控流程图 原始 FASTQ 文件 FastQC 质量评估 ✅ 合格 ❌ 不合格 质量控制报告（MultiQC） 数据清洗（cutadapt/Trimmomatic） 检查通过后才可进入下游比对 下游分析（BWA/GATK）

3.5 如何快速判读一份 FastQC 报告

报告里每个模块会显示绿色（通过）、黄色（警告）、红色（失败）。但我建议大家按这个顺序快速扫：

  1. 首先看 Per base sequence quality：如果整个图大部分区域低于 Q20，直接判定文库失败。
  2. 再看 Adapter Content：超过 5% 的 reads 含接头，必须切除。
  3. 然后看 Overrepresented sequences：如果占比最高的序列是 rRNA 或线粒体基因组，说明样本富集或建库有问题。
  4. 最后看 Per tile sequence quality：如果热图出现整列红色，建议联系测序中心。

个人经验：我用 FastQC 这么久，发现很多实验室只关注第一个模块就放行了。其实 Adapter Content 和 Overrepresented sequences 才是决定是否要切数据的关键。我团队里有一个质检清单，专门列这 4 项指标必须同时通过才能进入下游。

3.6 批量处理和报告汇总

一个项目动辄几十个样本，一个个看 HTML 报告效率太低。我推荐用 MultiQC 做批量的可视化汇总。MultiQC 会自动扫描目录下的 FastQC 报告，生成一张综合仪表盘。

multiqc ./qc_reports/ -o ./multiqc_report/

它会生成：

  • 所有样本的质量分布对比图
  • GC 含量分布
  • 序列重复率
  • 接头污染的整体情况

这样一眼就能看出哪几个样本是异常值。我记得有一次批量处理 96 个肿瘤 WES 样本，MultiQC 的柱状图里有两个样本的 Q30 明显偏低，单独排查发现是文库浓度估计偏差导致上样量不足。这种问题在小批量里很难发现。

3.7 关于 FASTQ 文件的命名规范

最后说一个容易被忽略的点：文件命名。GATK 管线要求输入文件路径必须清晰可辨。我个人习惯采用这种格式：

项目名_样本ID_库类型_L00X_R1_001.fastq.gz

例如：
Proj001_Sample_Blood_WES_L001_R1_001.fastq.gz
Proj001_Sample_Blood_WES_L001_R2_001.fastq.gz

这样做的好处：

  • 项目、样本、文库、lane、读长方向全在文件名里
  • 后续用脚本批量处理时，可以精准提取元信息
  • 多人协作时不会混淆
⚠️ 警告：
不要在文件名里用空格、中文、特殊符号（如 $ % &）。Linux 系统对这些字符敏感，GATK 也会因为路径格式问题报错。我遇到过有人把样本名写成 "Patient#01"，结果一堆 shell 脚本跑不了。

好了，FASTQ 格式和质量控制部分就这些。下一阶段我们会拿着 FastQC 的评估结果，决定是否需要做数据修剪。记住：质控不是可选项，是任何可复现 GATK 管线的第一道阀。

← 上一章 📖 返回目录 下一章 →