gatk-variant-calling

GATK — 胚系变异检测流程

概览

GATK（Genome Analysis Toolkit）实现了 GATK 最佳实践工作流，用于从 Illumina WGS 和 WES 数据中检测 SNP 和 indel。该流程对每个样本运行 HaplotypeCaller（生成 GVCF 文件），使用 GenomicsDBImport 汇总 GVCF，通过 GenotypeGVCFs 进行联合分型，并使用 VQSR（Variant Quality Score Recalibration，变异质量评分重校准）或硬过滤筛选变异。GATK 需要以经 BWA-MEM2 比对、标记重复并完成碱基质量评分重校准（BQSR）的 BAM 文件作为输入。它集成了 Picard 工具、samtools 和 bcftools，用于预处理和后处理。GATK4 工作流是 NIH/ENCODE 在研究和临床基因组学中进行胚系变异检测的标准。

适用场景

• 从 WGS 或 WES 样本中检测胚系 SNP 和 indel，用于群体遗传学或临床变异分析
• 在多个样本间运行联合分型，用于队列规模研究（家系、病例对照）
• 在运行 HaplotypeCaller 前进行碱基质量评分重校准（BQSR），以提高变异检测准确性
• 生成 GVCF 文件以便扩展队列：添加新样本时无需重新处理已有样本
• 生成变异检出集，以便后续使用 Ensembl VEP、ANNOVAR 或 SnpEff 进行注释
• 如果需要速度更快且准确性相当的深度学习方法，请改用 DeepVariant（Google）
• 如果需要快速检测变异且不进行基于组装的局部重比对，请改用 bcftools call

前置条件

• 软件：GATK4、Java 17+、samtools、BWA-MEM2
• 参考文件：基因组 FASTA + 已知变异 VCF（dbSNP、1000G、Mills indel）
• 输入：已标记重复并排序、带有 @RG 读组标头的 BAM（来自 BWA-MEM2）

安装前请检查：当前环境中可能已经有此工具（例如，在 pixi / conda 环境中）。先运行 command -v gatk；如果返回路径，请跳过下方的安装命令。在 pixi 项目中运行时，请使用 pixi run gatk 调用工具，而不是直接运行 gatk。

# Install GATK4
wget https://github.com/broadinstitute/gatk/releases/download/4.6.0.0/gatk-4.6.0.0.zip
unzip gatk-4.6.0.0.zip
export GATK="$PWD/gatk-4.6.0.0/gatk"

# Or with conda
conda install -c bioconda gatk4

# Verify
gatk --version
# GATK v4.6.0.0

# Download GATK resource bundle files (GRCh38)
# From gs://gcp-public-data--broad-references/hg38/v0/ (requires gsutil or Broad FTP)

预先访谈

编写任何分析代码前，先与用户确认以下事项。

decisions:
  - id: D1
    param: referenceBundle
    kind: derived
    source: upstream
    ask: "Which reference FASTA were the BAMs aligned to?"
    default: "carried from the alignment stage; must match exactly"

  - id: D2
    param: callingIntervals
    kind: required
    source: user
    ask: "Call across the whole genome, or restrict to capture targets or regions of interest?"
    default: null

  - id: D3
    param: cohortMode
    kind: required
    source: user
    ask: "Call each sample on its own, or emit GVCFs now and genotype the cohort jointly afterwards?"
    default: "joint genotyping when more than one sample is being analysed together"

  - id: D4
    param: filteringStrategy
    kind: required
    source: user
    depends_on: [D3]
    ask: "Filter calls with fixed hard thresholds, or train a variant-recalibration model on the cohort?"
    default: "hard filters; recalibration needs a cohort large enough to train on"

  - id: D5
    param: knownSites
    kind: optional
    source: user
    ask: "Annotate calls against a known-variant resource such as dbSNP?"
    default: "none"

  - id: D6
    param: callConfidence
    kind: optional
    source: user
    ask: "How confident must a genotype be before the variant is emitted at all?"
    default: 30

  - id: D7
    param: threadsAndHeap
    kind: never_ask
    source: data
    reason: "Pair-HMM threads and JVM heap size affect runtime and memory, not the calls"
    default: "4 threads, heap sized to the genome"

D4 取决于 D3，因为变异重校准是一种队列方法——对于单个样本，它没有可用于训练的数据，因此只有在回答队列相关问题后，过滤选择才有意义。

快速开始

GENOME="GRCh38.fa"
DBSNP="dbsnp_146.hg38.vcf.gz"

# Run HaplotypeCaller in GVCF mode
gatk HaplotypeCaller \
    -R $GENOME \
    -I sample1.markdup.bam \
    -O sample1.g.vcf.gz \
    -ERC GVCF \
    --dbsnp $DBSNP \
    --native-pair-hmm-threads 4

echo "GVCF: sample1.g.vcf.gz"

工作流

步骤 1：碱基质量评分重校准（BQSR）

在变异检测前纠正碱基质量评分中的系统误差。

GENOME="GRCh38.fa"
KNOWN_SITES="dbsnp_146.hg38.vcf.gz Mills_and_1000G_gold_standard.indels.hg38.vcf.gz"
KNOWN_FLAGS=$(printf -- '--known-sites %s ' $KNOWN_SITES)

# Step 1a: Build recalibration table
gatk BaseRecalibrator \
    -R $GENOME \
    -I sample1.markdup.bam \
    $KNOWN_FLAGS \
    -O sample1.recal.table

# Step 1b: Apply recalibration
gatk ApplyBQSR \
    -R $GENOME \
    -I sample1.markdup.bam \
    --bqsr-recal-file sample1.recal.table \
    -O sample1.bqsr.bam

echo "BQSR BAM: sample1.bqsr.bam"
samtools flagstat sample1.bqsr.bam | head -3

步骤 2：使用 HaplotypeCaller 检测变异（GVCF 模式）

对每个样本分别运行变异检测，生成用于联合分型的中间 GVCF。

# HaplotypeCaller in GVCF mode (recommended for cohort analysis)
gatk HaplotypeCaller \
    -R GRCh38.fa \
    -I sample1.bqsr.bam \
    -O gvcfs/sample1.g.vcf.gz \
    -ERC GVCF \
    --dbsnp dbsnp_146.hg38.vcf.gz \
    --native-pair-hmm-threads 4

# For WES: specify target intervals
# gatk HaplotypeCaller ... -L exome_targets.interval_list --interval-padding 100

echo "GVCF: gvcfs/sample1.g.vcf.gz"
zcat gvcfs/sample1.g.vcf.gz | grep -v "^#" | wc -l

步骤 3：使用 GenomicsDBImport 汇总 GVCF

合并各样本的 GVCF，以便高效进行联合分型。

# Create sample map file: sample_name\tpath_to_gvcf
printf "ctrl_1\tgvcfs/ctrl_1.g.vcf.gz\n" > sample_map.txt
printf "ctrl_2\tgvcfs/ctrl_2.g.vcf.gz\n" >> sample_map.txt
printf "treat_1\tgvcfs/treat_1.g.vcf.gz\n" >> sample_map.txt
printf "treat_2\tgvcfs/treat_2.g.vcf.gz\n" >> sample_map.txt

# Import GVCFs into GenomicsDB for each chromosome
for CHR in chr1 chr2 chr3; do
    gatk GenomicsDBImport \
        --sample-name-map sample_map.txt \
        --genomicsdb-workspace-path genomicsdb/${CHR} \
        -L $CHR \
        --reader-threads 4
done

echo "GenomicsDB created for $(ls genomicsdb/ | wc -l) chromosomes"

步骤 4：使用 GenotypeGVCFs 进行联合分型

在 GenomicsDB 中对所有样本同时进行分型。

# Joint genotype all samples
mkdir -p vcfs
for CHR in chr1 chr2 chr3; do
    gatk GenotypeGVCFs \
        -R GRCh38.fa \
        -V gendb://genomicsdb/${CHR} \
        --dbsnp dbsnp_146.hg38.vcf.gz \
        -O vcfs/cohort_${CHR}.vcf.gz
done

# Merge per-chromosome VCFs
gatk MergeVcfs \
    $(ls vcfs/cohort_chr*.vcf.gz | sed 's/^/-I /') \
    -O vcfs/cohort_all.vcf.gz

echo "Joint genotyping complete: vcfs/cohort_all.vcf.gz"
gatk CountVariants -V vcfs/cohort_all.vcf.gz

步骤 5：变异过滤（硬过滤）

对于 VQSR 效力不足的小型队列，应用硬过滤。

# Separate SNPs and indels
gatk SelectVariants -V vcfs/cohort_all.vcf.gz --select-type-to-include SNP -O vcfs/snps.vcf.gz
gatk SelectVariants -V vcfs/cohort_all.vcf.gz --select-type-to-include INDEL -O vcfs/indels.vcf.gz

# Apply hard filters: SNPs
gatk VariantFiltration \
    -V vcfs/snps.vcf.gz \
    --filter-expression "QD < 2.0" --filter-name "QD2" \
    --filter-expression "FS > 60.0" --filter-name "FS60" \
    --filter-expression "MQ < 40.0" --filter-name "MQ40" \
    --filter-expression "MQRankSum < -12.5" --filter-name "MQRankSum-12.5" \
    -O vcfs/snps_filtered.vcf.gz

# Apply hard filters: Indels
gatk VariantFiltration \
    -V vcfs/indels.vcf.gz \
    --filter-expression "QD < 2.0" --filter-name "QD2" \
    --filter-expression "FS > 200.0" --filter-name "FS200" \
    --filter-expression "ReadPosRankSum < -20.0" --filter-name "ReadPosRankSum-20" \
    -O vcfs/indels_filtered.vcf.gz

# Merge filtered SNPs + indels
gatk MergeVcfs -I vcfs/snps_filtered.vcf.gz -I vcfs/indels_filtered.vcf.gz \
    -O vcfs/cohort_filtered.vcf.gz

echo "PASS variants: $(bcftools view -f PASS vcfs/cohort_filtered.vcf.gz | grep -v '^#' | wc -l)"

步骤 6：使用 Python 解析 VCF 结果

提取变异，添加基因信息注释，并整理成 DataFrame。

import subprocess
import pandas as pd
import io

# Use bcftools query to extract fields from filtered VCF
result = subprocess.run(
    ["bcftools", "query",
     "-f", "%CHROM\t%POS\t%ID\t%REF\t%ALT\t%QUAL\t%FILTER\t%INFO/QD\t%INFO/FS\n",
     "-i", "FILTER='PASS'",
     "vcfs/cohort_filtered.vcf.gz"],
    capture_output=True, text=True
)

cols = ["CHR", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "QD", "FS"]
df = pd.read_csv(io.StringIO(result.stdout), sep="\t", names=cols)
df["QUAL"] = pd.to_numeric(df["QUAL"], errors="coerce")
df["QD"] = pd.to_numeric(df["QD"], errors="coerce")

print(f"PASS variants: {len(df)}")
print(f"SNPs:  {(df['REF'].str.len() == 1) & (df['ALT'].str.len() == 1)).sum()}")
print(f"Indels: {((df['REF'].str.len() > 1) | (df['ALT'].str.len() > 1)).sum()}")
print(df.head())
df.to_csv("pass_variants.tsv", sep="\t", index=False)

关键参数

参数	默认值	范围/选项	作用
-ERC	NONE	GVCF、BP_RESOLUTION	启用参考置信度模式；队列工作流使用 GVCF
--native-pair-hmm-threads	4	1–32	HaplotypeCaller 中 pair-HMM 使用的线程数（CPU 消耗最高的步骤）
-L	全基因组	区间文件或 chr	将变异检测限制在指定区间（外显子组靶区、BED 区域）
--dbsnp	—	VCF 路径	用于在输出中标注 rsID 的 dbSNP VCF
--stand-call-conf	30	0–100	输出变异调用所需的最低基因型质量分数
-G	StandardAnnotation	注释组	要应用的注释模块（HaplotypeCaller 使用 StandardHCAnnotation）
--sample-name-map	—	TSV 文件	GenomicsDBImport 的样本到 GVCF 映射
--reader-threads	1	1–16	GenomicsDBImport 读取时使用的线程数
--filter-expression	—	JEXL 表达式	硬过滤表达式（例如 "QD < 2.0"）
--java-options	—	-Xmx4g	Java 堆大小；处理大型基因组时使用 -Xmx16g

常用配方

配方 1：单样本变异检测（无队列）

# For a single sample, skip GenomicsDBImport and call directly
gatk HaplotypeCaller \
    -R GRCh38.fa \
    -I sample.bqsr.bam \
    -O sample.vcf.gz \
    --dbsnp dbsnp_146.hg38.vcf.gz \
    --native-pair-hmm-threads 8

# Hard filter the single-sample VCF
gatk VariantFiltration \
    -V sample.vcf.gz \
    --filter-expression "QD < 2.0" --filter-name "QD2" \
    --filter-expression "FS > 60.0" --filter-name "FS60" \
    -O sample_filtered.vcf.gz

echo "PASS variants: $(bcftools view -f PASS sample_filtered.vcf.gz | grep -v '^#' | wc -l)"

配方 2：使用 Snakemake 运行完整流程

# Snakefile — GATK BQSR + HaplotypeCaller
configfile: "config.yaml"
SAMPLES = config["samples"]
GENOME  = config["genome"]
DBSNP   = config["dbsnp"]
KNOWN   = config["known_sites"]

rule all:
    input: expand("gvcfs/{sample}.g.vcf.gz", sample=SAMPLES)

rule bqsr:
    input: bam="markdup/{sample}.markdup.bam"
    output: bam="bqsr/{sample}.bqsr.bam"
    shell:
        """
        gatk BaseRecalibrator -R {GENOME} -I {input.bam} \
             --known-sites {KNOWN} -O bqsr/{wildcards.sample}.recal.table &&
        gatk ApplyBQSR -R {GENOME} -I {input.bam} \
             --bqsr-recal-file bqsr/{wildcards.sample}.recal.table \
             -O {output.bam}
        """

rule haplotypecaller:
    input: bam="bqsr/{sample}.bqsr.bam"
    output: gvcf="gvcfs/{sample}.g.vcf.gz"
    threads: 4
    shell:
        "gatk HaplotypeCaller -R {GENOME} -I {input.bam} -O {output.gvcf} "
        "-ERC GVCF --dbsnp {DBSNP} --native-pair-hmm-threads {threads}"

预期输出

输出	格式	说明
*.g.vcf.gz	GVCF	每个样本的 GVCF，包含参考置信度区块；作为 GenomicsDBImport 的输入
cohort_all.vcf.gz	VCF	联合基因分型的多样本 VCF；未过滤
cohort_filtered.vcf.gz	VCF	过滤后的 VCF；通过过滤的变异其 FILTER=PASS
*.recal.table	BQSR	BaseRecalibrator 生成的碱基质量重校准表
*.bqsr.bam	BAM	重校准后的 BAM；用作 HaplotypeCaller 的输入
genomicsdb/	目录	用于联合基因分型的按染色体划分的 GenomicsDB 工作区

故障排除

问题	原因	解决方案
SAM/BAM file has no @RG header	BWA-MEM2 未提供读组信息	使用 -R "@RG\tID:...\tSM:...\tPL:ILLUMINA" 重新比对
Java OutOfMemoryError	堆大小不足	添加 --java-options "-Xmx16g" 或更大的值
HaplotypeCaller 运行非常慢	HMM 为单线程运行	添加 --native-pair-hmm-threads 8；使用区间列表按染色体并行处理
GVCF 输出为空	区间错误或区域内没有 reads	检查 samtools view -c sample.bam chrN 获取 reads 数量
VQSR 失败（< 10k 个变异）	用于训练的变异数量过少	对小型队列或外显子组使用硬过滤，而不是 VQSR
GenomicsDB 导入时因路径已存在而失败	GenomicsDB 工作区已存在	重新运行前删除现有工作区：rm -rf genomicsdb/chr1
IndexOutOfBoundsException	BAM 与参考基因组的染色体名称不匹配	确保基因组 FASTA 与 BAM 使用相同的染色体命名方式（chr1 与 1）
缺少 BCFtools/tabix 索引	未创建 Tabix 索引（.tbi）	运行 gatk IndexFeatureFile -I file.vcf.gz 或 tabix -p vcf file.vcf.gz

参考资料

• GATK 文档 — 官方流程指南和工具文档
• GATK 生殖系短变异最佳实践 — 分步操作规范
• GATK GitHub：broadinstitute/gatk — 源代码、发行版和资源包
• Van der Auwera GA 和 O'Connor BD（2020）Genomics in the Cloud — O'Reilly Media；全面介绍 GATK4 的指南