/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)

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

Step 3: Consolidate GVCFs with GenomicsDBImport

Merge per-sample GVCFs for efficient joint genotyping.

# 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"

Step 4: Joint Genotyping with GenotypeGVCFs

Genotype all samples simultaneously across the 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

Step 5: Variant Filtration (Hard Filters)

Apply hard filters for small cohorts where VQSR is underpowered.

# Separate SNPs and indels
gatk SelectVariants -V vcfs/cohort_all.vcf.gz --select-type-to