GATK质控概述与常见失败模式
讲GATK质控,说白了就是盯住两个问题:一是你的数据到底能不能信,二是Pipeline到底在哪儿断的。我个人习惯在跑完一套流程后,先不看结果,而是翻一遍QC报告和日志结尾——很多坑在正式产出之前就已经埋下了。
1.1 先读懂QC报告
QC报告不是给你看的装饰品,它直接告诉你当前数据有没有"病"。我印象最深的一次,BQSR之后的碱基质量分布图看起来完美,结果下游却出了大批假阳性——后来发现是因为参考基因组用了旧版本,BQSR算出来的矫正系数全偏了。嗯,从那以后我每到一个新项目组,第一件事就是确认参考基因组的checksum。
基础报告你得看这几个关键字段:
| 报告项目 | 检查什么 | 典型异常信号 |
|---|---|---|
| BaseRecalibrator 特征表 | 碱基质量与协变量分布的偏差 | 某个质量值区间的覆盖度突然跳升 |
| BQSR前后质量分布 | 中位数质量是否提升,低质量尾部是否收窄 | BQSR后反而多了低质量碱基 |
| 覆盖度统计 | 均匀性、深度中位数 | 某条染色体深度为0 |
| 比对率 | 整体比对比例 | 比对率低于90% |
看报告时我习惯用现成的工具快速扫一遍:
# 查看BQSR后的质量汇总
gatk AnalyzeCovariates \
-before recal_data.table \
-after post_recal_data.table \
-plots BQSR_report.pdf
# 快速检查比对情况
samtools flagstat sample.bam | head -20
# 查看覆盖度分布
samtools depth sample.bam | awk '{sum+=$3;cnt++} END{print sum/cnt}'
我的习惯:在正式跑HaplotypeCaller之前,先把BAM按chr21单独过一遍BQSR+Call,全程开着实时日志。这样能在10分钟内发现大部分配置性问题,避免浪费全基因组时间。
1.2 Pipeline崩溃的典型场景
Pipeline崩溃从来不是随机事件,它一定有个明确的触发点。我统计过自己项目里遇到的上百次失败,80%可以归到下面几张表里。
| 失败场景 | 根因 | 症状 | 快速修复 |
|---|---|---|---|
| 内存爆掉 | 线程数开太多 / 堆内存不够 | OutOfMemoryError / 被系统kill | 减少并行,或增加 -Xmx |
| 参考基因组版本不一致 | BAM的@SQ和fasta的contig名称对不上 | MISSING_REF / 索引失败 | 统一参考版本,重建dict |
| BAM索引损坏 | 拷贝中断 / 磁盘满 | FileNotFoundException 或 0长度bai | 重新排序并索引 |
| 样本名冲突 | 同一批次出现两个同名样本 | 重复样本错误 / Merge失败 | 修改ID或分开处理 |
| VCF文件不完整 | 前一步写入时被杀 | TBI索引异常 / 缺失末尾记录 | 重新合并或补跑 |
这里我特别想说一下内存问题——GATK4比GATK3在内存管理上好了很多,但如果你不加限制地用默认线程数(等于CPU核数),很容易把服务器撑爆。我曾经的教训是:在集群上跑,一定要在任务脚本里显式写明-Xmx和并行线程数,别让GATK自己去猜。
1.3 如何快速定位失败原因
别慌,先看最后几行日志。GATK的异常信息其实很直白,只是很多人被海量WARNING淹没了注意力。给你一个标准检查流程:
- 打开日志文件,从末尾往上翻,找到第一个ERROR或FATAL行。
- 记录错误码和异常类名(比如IllegalArgumentException)。
- 检查前一步的输出文件是否完整(用
gzip -t校验gzip,用bcftools index -s检查VCF)。 - 对比参考基因组的md5值,确保BAM的@SQ头完全一致。
- 如果报的是OOM,看系统日志里是否有"Out of memory: Kill process"。
# 快速检查gzip文件完整性
for f in *.bam; do
if ! gzip -t "$f" 2>/dev/null; then
echo "CORRUPT: $f"
fi
done
# 查看BAM头里的参考基因组
samtools view -H sample.bam | grep "^@SQ" | head -3
# 对比参考序列名称
samtools view -H sample.bam | grep "^@SQ" | cut -f2 | sort > bam_ref.txt
grep "^>" reference.fasta | sed 's/>//' | sort > fasta_ref.txt
diff bam_ref.txt fasta_ref.txt
经验之谈:80%的崩溃都能在日志的前5分钟找到线索。别急着重跑,先花10分钟读崩溃日志——这比盲目重跑节省的时间要多得多。
1.4 质控核心逻辑图
下面这张图把GATK标准Germline流程中的质控节点和常见断裂点画了出来。橙色节点标注了最容易出事的步骤,箭头上的红色叉号表示那一步经常崩溃。
1.5 一些实用技巧
小心日志里的WARNING:GATK会打印大量"WARNING: some reads have invalid mapping quality"之类的信息,很多同学看到一堆WARNING就慌了,其实大部分不影响最终结果。真正致命的只有ERROR或FATAL。
我习惯的做法:在Pipeline启动前,把所有中间输出目录清空一遍,避免残留文件干扰。另外,每次改参数后,先在1%的数据上跑通,再投入全量——这个习惯帮我省下至少一周的工期。
最后说一句:质控不是一次性的。你在每步结束后都可以检查产出物的合法性。我常用的是一个三行脚本——检查bam的头、检查vcf的索引、检查日志末尾是否有ERROR——每天早晨看一眼,心里就有底了。