第三章 输入数据准备：FASTQ 质量控制、BAM/SAM 格式规范、参考基因组索引构建

大家好，我是老张。这一章咱们聊聊 GATK 管线的“地基”——输入数据的准备工作。别看这些都是基础活儿，但我敢说，至少一半的 GATK 报错或结果异常，根源都在这一步没做好。

我个人习惯是：在拿到测序中心提供的 FASTQ 后，先花半天时间做一轮质量控制，然后把格式搞清楚，最后才去建索引。这三个环节环环相扣，缺一不可。

一、FASTQ 质量控制：堵住垃圾数据进管线

FASTQ 文件是测序仪直接吐出来的原始数据。它里面藏着什么问题呢？

  • 低质量碱基：测序质量跟位置有关，末尾经常掉链子。
  • 接头污染：如果片段长度比测序读长短，接头序列就会被读出来。
  • 重复序列：PCR 扩增带来的重复，会干扰变异检测的统计模型。
  • N 过多：未知碱基太多，下游比对直接歇菜。

我常用的工具是 FastQC（看报告）加 fastp（做清洗）。为什么不推荐老旧的 Trimmomatic？因为 fastp 单线程就能处理压缩文件，而且自动识别接头格式。我在项目中遇到过两次，对方给的 FASTQ 是 Illumina NextSeq 的双端数据，接头藏在 Phred+33 质量值里，Trimmomatic 需要手动指定 adapter 文件，而 fastp 一跑就能自动检测出来。

我的日常处理流程：
  1. 用 FastQC 生成 HTML 报告，重点关注「Per base sequence quality」和「Adapter Content」。
  2. 用 fastp 一键清洗，跑完再跑一次 FastQC 确认干净。
  3. 关键参数：--qualified_quality_phred 20 和 --length_required 50。

举个例子，我处理人全基因组 30x 数据时的命令：

fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
      -o sample_R1.clean.fastq.gz -O sample_R2.clean.fastq.gz \
      --detect_adapter_for_pe \
      --qualified_quality_phred 20 \
      --unqualified_percent_limit 40 \
      --n_base_limit 5 \
      --length_required 50 \
      --cut_front --cut_front_mean_quality 20 \
      --cut_tail --cut_tail_mean_quality 20 \
      -j fastp_report.json -h fastp_report.html

嗯，这里要注意：--cut_front 和 --cut_tail 是滑动窗口切除低质量碱基，适用于读段开头质量不稳（比如 WGBS 数据），如果你做的是标准 WES 可以不开。为什么呢？因为捕获试剂盒的酶切步骤不会破坏 5' 端质量，切多了反而损失覆盖度。

我曾经踩过的坑： 跑 GATK HaplotypeCaller 之前没有去除 PCR 重复，结果在纯合区域检出了 20% 的杂合位点。后来发现是 DuplicateReadSet 标记（MarkDuplicates）没做，那些 PCR 重复被当作独立证据用于组装。所以，一定要在比对 BAM 之后、变异检测之前做 MarkDuplicates，或者用 fastp 自带的 --dedup 参数在原始 FASTQ 阶段先标记一遍。

二、BAM/SAM 格式规范：GATK 的“身份证”要求

BAM/SAM 是比对后的标准格式。GATK 对它们有非常严格的要求，很多报错都跟格式不满足规范有关。说白了，GATK 不信任任何来源不明的 BAM，它有一套自己的校验逻辑。

我记得第一次跑 GATK 时，直接拿 BWA 生成的原生 BAM 输进去，结果 PrintReads 阶段报错说「SAM file has no header， or header has no @SQ lines」。当时我还以为是 GATK 版本问题，折腾了半天才发现 BWA 默认不输出头信息。从那以后，我对 BAM 头的完整性特别敏感。

2.1 必须满足的硬性规范

字段/条件 要求 检查工具
@HD 头 必须包含 VN:1.x 和 SO:coordinate（排序） samtools head -1
@SQ 头 每一条参考序列都要有 SN 和 LN，且顺序必须与参考索引一致 samtools view -H
@RG 头 必须有 ID、SM、LB、PL，而且是唯一的 GATK 自带 ValidateSamFile
比对标志位 必须符合 SAM 规范，0x4 标记 unmapped 的 read 不得出现在 paired BAM 的 mate 位置上 samtools flagstat
CIGAR 操作 不允许包含 I/D 混合的非法组合（比如 2I3D） GATK SplitNCigarReads

最常踩的雷：@RG 缺失或重复。 GATK 在 BaseRecalibrator 阶段需要 ReadGroup 信息来区分不同文库，如果缺失会直接报 MISSING_READ_GROUP 错误。而且同一个 BAM 里不能有两个相同 ID 的 @RG，否则工具会混乱。

添加 @RG 的正确姿势是在 BWA 比对时就加上参数：

bwa mem -t 8 -R '@RG\tID:Sample1\tSM:Sample1\tLB:Lib1\tPL:ILLUMINA' \
  reference.fa R1.fastq.gz R2.fastq.gz | samtools sort -o sorted.bam

如果已经比对完了怎么办？可以用 samtools addreplacerg 修改，但要注意同时更新头信息和所有记录的 RG tag，否则 GATK 会认为 BAM 被篡改过。

2.2 我推荐的 BAM 预处理流程

  1. 排序： samtools sort -@ 8 -m 4G，按坐标排序。
  2. 标记重复： GATK MarkDuplicates 或 Picard 的相同工具，输出一个包含 .metrics.txt 的文件。
  3. 建立索引： samtools index sorted.bam，生成 .bai 文件。
  4. 校验： gatk ValidateSamFile -I sorted.bam -M SUMMARY

为什么我坚持要求验证呢？因为 GATK 的工具链 —— 尤其是 HaplotypeCaller —— 在读入 BAM 时会对每条 read 的 CIGAR 做合法性检查，非法 CIGAR 会导致读取中断。我在一次 RNA-seq 数据分析中就遇到过，STAR 比对产生的 BAM 在剪接位点用了 N 操作，但 GATK 默认用 SplitNCigarReads 处理 RNA-seq 数据之前必须先做一次 NCIGAR 的拆分，否则直接报错。

三、参考基因组索引构建：一步错，全盘错

参考基因组索引是 GATK 的“地图”。地图画错了，在哪条染色体上找变异都会南辕北辙。GATK 需要三种索引：FASTA 本身的 .fai、.dict（字典文件），以及 bwa 索引（如果是 BWA 比对）。

我的习惯是用一套脚本同时生成所有索引：

# 假设参考基因组为 Homo_sapiens_assembly38.fasta （来自 GATK Resource Bundle）
# 1. 生成 .fai
samtools faidx Homo_sapiens_assembly38.fasta

# 2. 生成 .dict
gatk CreateSequenceDictionary -R Homo_sapiens_assembly38.fasta \
      -O Homo_sapiens_assembly38.dict

# 3. 生成 BWA 索引（如果需要比对）
bwa index -a bwtsw Homo_sapiens_assembly38.fasta

我的建议： 永远从 GATK Resource Bundle 下载预构建好的参考基因组（如 hg38 的 uction 版），而不是自己从 UCSC/NCBI 下载 raw FASTA 再手动构建。为什么呢？因为 Resource Bundle 里的 FASTA 已经去掉了 ALT contigs 和 decoy 序列，而且序列名符合 GATK 约定（比如 chr1 vs 1）。我自己以前从 NCBI 下载 GRCh38.p13，对比时发现序列名是 NC_000001.11，结果 GATK 不认识，只好重新索引。

3.1 索引之间的交叉验证

建好索引后，一定要做两件事：

  • 核对 .fai 和 .dict 的序列数是否一致——用 diff <(cut -f1 .fai) <(grep "^@SQ" .dict | cut -f2 | cut -d: -f2)。
  • 确认 BAM 的 @SQ 头与参考序列完全一致——包括大小写、顺序。

我记得有一回在提交 SLURM 作业时，集群管理员把参考基因组版本混了：WES 数据用的是 hg38 比对，但 BaseRecalibrator 用的是 hg19 的索引。结果跑出来的 BQSR 表格全是负值，校准后 BAM 质量反倒降低。检查半天才发现是参考基因组不匹配。所以，索引文件一定要跟 FASTA 文件放在同一个目录下，并且随 FASTA 一起版本控制。

知识结构总览

输入数据准备三阶段 FASTQ 质量控制 FastQC + fastp BAM/SAM 格式规范 @RG、CIGAR、排序 参考基因组索引 .fai + .dict + bwa 三阶段联合验证 FASTQ 清洗 → BAM 格式校验 → 索引一致性检查 确保 GATK 管线零异常输入 进入 GATK 变异检测流程

一句话总结本章要点：

  • FASTQ 要经过质量清洗和接头去除，再用 FastQC 复验。
  • BAM 必须包含完整 @RG 头、正确排序、去重复，并通过 ValidateSamFile。
  • 参考基因组索引三件套（.fai + .dict + bwa 索引）缺一不可，且必须与 BAM 的 @SQ 头完全一致。

做好了这三件事，你后面调 GATK 参数时就可以把注意力放在真正的生物学问题上了，而不是被各种 “IllegalArgumentException: SAMFileReader does not have any index” 这种低级报错牵着鼻子走。嗯，这一节就到这儿，咱们下节见。

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