数据获取与质量控制:原始测序数据格式(FASTQ)、质量评估(FastQC)
拿到原始测序数据,第一件事不是急着往下跑。我见过太多人从 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行。每个字符代表一个碱基的测序置信度。
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,立即处理 |
有的实验室用老测序仪,偶尔出现单 tile 质量下降。我曾经遇到一个项目,16 个样本里有 3 个的 tile 热图出现一整列橘红色。如果直接跳过不看,下游变异检测会在这个区域引入系统性假阳性。
3.4 从 FASTQ 到 FastQC 的工作流程
我一般把这一阶段的流程画成下图,方便团队理解:
3.5 如何快速判读一份 FastQC 报告
报告里每个模块会显示绿色(通过)、黄色(警告)、红色(失败)。但我建议大家按这个顺序快速扫:
- 首先看 Per base sequence quality:如果整个图大部分区域低于 Q20,直接判定文库失败。
- 再看 Adapter Content:超过 5% 的 reads 含接头,必须切除。
- 然后看 Overrepresented sequences:如果占比最高的序列是 rRNA 或线粒体基因组,说明样本富集或建库有问题。
- 最后看 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 管线的第一道阀。