📚 课程目录第2章 GATK工具套件介绍

2. GATK工具套件介绍

说到GATK，做生物信息的人应该都绕不开它。这玩意儿，说白了就是基因组分析领域的工业级标准。特别是做重测序、找变异的朋友，GATK几乎是必经之路。

我记得2013年刚接手第一个全外显子组项目时，用的还是GATK3的版本。那时候是MapReduce架构，跑起来又慢又麻烦，动不动就内存炸了。后来GATK4在2018年发布，架构来了个大翻新，简直像换了个引擎。我个人非常推崇GATK4——它更轻、更快，而且新加的Spark模式让大基因组分析变得舒服多了。

2.1 GATK4核心架构——它为什么变得更好了？

一个工具好不好用，架构决定一半。GATK4从GATK3的“Queue”架构切换到了统一的“Tool”架构。什么意思呢？

以前：各个工具散落，每次分析要组合很多分散的模块
现在：每个分析任务打包成一个独立工具（Tool），用起来很直接

GATK4取消了Python写的Queue调度层，所有工具变成了标准的Java命令行。开发团队也加入了对Spark的原生支持。

嗯，这里要注意：Spark模式不是必须的。小项目用单机模式就够，大项目（比如1000+样本）才值得上Spark。

我的习惯：小团队做20-50个样本，就用GATK4的常规模式。样本量破200了，我会考虑用Spark或者走云计算。

GATK4启动的方式很简单。下载好GATK的jar包或压缩包后，直接调主程序就行：

# GATK4 最基本的使用方式
gatk --list  # 列出所有可用工具

gatk HaplotypeCaller \
    -R reference.fasta \
    -I sample.bam \
    -O sample.g.vcf.gz

看到没？一条命令就是一个独立分析动作。比起GATK3那套复杂的workflow配置，清爽多了。

2.2 核心工具模块——你常用哪些？

GATK4工具虽多，但高频使用的其实就那么20来个。我这里按功能把它们分成4大模块。

模块分类	核心工具	主要用途
数据预处理	MarkDuplicates, BaseRecalibrator, ApplyBQSR, PrintReads, SortSam	BAM文件标重复、碱基质量分数校正、排序
变异发现	HaplotypeCaller, Mutect2, GenotypeGVCFs, CombineGVCFs	单样本变异检测、联合基因分型、体细胞突变检测
变异过滤	VariantFiltration, SelectVariants, HardFilter, CNNScoreVariants	硬过滤与基于模型的过滤、子集提取
功能注释	VariantAnnotator, Funcotator, VariantsToTable	变异功能注释、生成分析报表

其中HaplotypeCaller是变异发现的核心引擎。我做联合分型项目时，大部分时间都跟它打交道。不过它的计算资源消耗也不小，一个全基因组样本大概需要8-16GB内存。

我在项目中遇到过一种情况：样本bam文件排序顺序不对，导致HaplotypeCaller报错“SAM headers don't match”。后来养成了习惯——跑GATK前，一定先用samtools做一次排序校验。

经验之谈：我见过不少新手直接拿未排序的BAM跑HaplotypeCaller，结果一跑就是好几个小时，最后报错终止。白费力气不说，还浪费计算时间。排序确认跑起来很快，千万别省这步。

来看看GATK4的整体模块关系。我画的这张图能帮你一见看清各模块的定位：

2.3 GATK最佳实践工作流——按步骤来，别跳步

Broad Institute提出了两套官方最佳实践（Best Practices）：

种系变异发现流程（Germline）：适用于正常组织、群体变异分析
体细胞变异发现流程（Somatic）：适用于肿瘤-正常配对分析

咱们的HapotypeCaller联合分型项目走的是种系流程，具体步骤如下：

第一步：BAM预处理

# 标记重复
gatk MarkDuplicates \
    -I sample.bam \
    -O sample_md.bam \
    -M sample_metrics.txt

# 碱基质量分数校正（BQSR）
gatk BaseRecalibrator \
    -R reference.fasta \
    -I sample_md.bam \
    --known-sites dbsnp.vcf.gz \
    --known-sites mills.vcf.gz \
    -O recal.table

gatk ApplyBQSR \
    -R reference.fasta \
    -I sample_md.bam \
    --bqsr-recal-file recal.table \
    -O sample_recal.bam

第二步：单样本GVCF产出

gatk HaplotypeCaller \
    -R reference.fasta \
    -I sample_recal.bam \
    -O sample.g.vcf.gz \
    -ERC GVCF

第三步：多样本联合分型

# 先合并所有样本的GVCF
gatk GenomicsDBImport \
    -R reference.fasta \
    -V sample1.g.vcf.gz \
    -V sample2.g.vcf.gz \
    --genomicsdb-workspace-path my_database \
    -L intervals.list

# 然后联合分型
gatk GenotypeGVCFs \
    -R reference.fasta \
    -V gendb://my_database \
    -O cohort.vcf.gz

第四步：VCF过滤

gatk VariantFiltration \
    -R reference.fasta \
    -V cohort.vcf.gz \
    -O cohort_filtered.vcf.gz \
    --filter-name "QD_fail" --filter-expression "QD < 2.0" \
    --filter-name "FS_fail" --filter-expression "FS > 60.0" \
    --filter-name "MQ_fail" --filter-expression "MQ < 40.0"

一个省时间的建议：在第二步和第三步之间，如果样本量很大（比如超过100个），我会考虑用slice方式分染色体并行。每条染色体建一个GenomicsDB，最后合并。跑起来快很多。

2.4 实战注意事项——我自己踩过的坑

参考基因组必须一致——所有样本必须使用同一个版本的参考基因组，否则联合分型会出严重错误
GVCF模式不可省略——很多人贪快直接输出VCF，但联合分型需要GVCF模式的完整中间信息
数据量大时关注磁盘IO——GenomicsDBImport是IO密集型操作，SSD比HDD能快数倍

我曾经犯过一个错误：有一次做200个样本的联合分型，忘了加-L bed文件指定目标区域。结果GenomicsDBImport跑了一整天还没跑完——它在全基因组上建数据库，数据量爆炸了。从那以后，每次跑多样本都会确认区间文件配对了没。

GATK4的功能足够强大，但也需要你理解它的设计哲学。说白了，它追求的是可重复性、可扩展性和工业级稳定性。只要流程搭对了，跑起来很稳当。