/gatk-variant-calling

name: "gatk-variant-calling"
description: "GATK Best Practices for germline SNP/indel calling from WGS/WES BAMs. Per-sample HaplotypeCaller GVCFs, GenomicsDBImport, GenotypeGVCFs joint calling, VQSR or hard filters. Requires BWA-MEM2-aligned, markdup, BQSR BAMs. Use DeepVariant for a faster DL alternative; GATK is the NIH/ENCODE standard."
license: "BSD-3-Clause"

GATK — Germline Variant Calling Pipeline

Overview

GATK (Genome Analysis Toolkit) implements the GATK Best Practices workflow for calling SNPs and indels from Illumina WGS and WES data. The pipeline runs HaplotypeCaller per sample (producing GVCF files), consolidates GVCFs with GenomicsDBImport, performs joint genotyping with GenotypeGVCFs, and filters variants with VQSR (Variant Quality Score Recalibration) or hard filters. GATK requires BWA-MEM2-aligned, duplicate-marked, and base quality score recalibrated (BQSR) BAM files as input. It integrates with Picard tools, samtools, and bcftools for pre- and post-processing. The GATK4 workflow is the NIH/ENCODE standard for germline variant calling in research and clinical genomics.

When to Use

• Calling germline SNPs and indels from WGS or WES samples for population genetics or clinical variant analysis
• Running joint genotyping across multiple samples for cohort-scale studies (families, case-control)
• Applying base quality score recalibration (BQSR) to improve variant calling accuracy before HaplotypeCaller
• Generating GVCF files for scalable cohort expansion: add new samples without reprocessing existing ones
• Producing variant call sets for downstream annotation with Ensembl VEP, ANNOVAR, or SnpEff
• Use **DeepVariant** (Google) instead for a faster deep-learning approach with comparable accuracy
• Use **bcftools call** instead for rapid variant calling without assembly-based local realignment

Prerequisites

• **Software**: GATK4, Java 17+, samtools, BWA-MEM2
• **Reference files**: genome FASTA + known variants VCF (dbSNP, 1000G, Mills indels)
• **Input**: duplicate-marked, sorted BAM with `@RG` read group headers (from BWA-MEM2)

> **Check before installing**: The tool may already be available in the current environment (e.g., inside a `pixi` / `conda` env). Run `command -v gatk` first and skip the install commands below if it returns a path. When running inside a pixi project, invoke the tool via `pixi run gatk` rather than bare `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)

Pre-flight Interview

Settle these with the user before writing any analysis code.

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 hangs on D3 because variant recalibration is a cohort method - it has nothing to train on for a single sample, so the filtering choice is only meaningful once the cohort question is answered.

Quick Start

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"

Workflow

Step 1: Base Quality Score Recalibration (BQSR)

Correct systematic errors in base quality scores before variant calling.

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

Step 2: Call Variants with HaplotypeCaller (GVCF Mode)

Run per-sample variant calling, producing an intermediate GVCF for joint genotyping.

# 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.vc