jaechang-hits/gatk-variant-calling

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

Step 6: Parse VCF Results with Python

Extract variants, annotate with gene info, and prepare a 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)

Key Parameters

| Parameter | Default | Range/Options | Effect | |-----------|---------|---------------|--------| | -ERC | NONE | GVCF, BP_RESOLUTION | Emit reference confidence mode; GVCF for cohort workflows | | --native-pair-hmm-threads | 4 | 1–32 | Threads for pair-HMM in HaplotypeCaller (most CPU-intensive step) | | -L | whole genome | interval file or chr | Restrict calling to intervals (exome targets, BED regions) | | --dbsnp | — | VCF path | dbSNP VCF for rsID annotation in output | | --stand-call-conf | 30 | 0–100 | Min genotype quality score to emit a variant call | | -G | StandardAnnotation | annotation group | Annotation modules to apply (StandardHCAnnotation for HaplotypeCaller) | | --sample-name-map | — | TSV file | Sample-to-GVCF mapping for GenomicsDBImport | | --reader-threads | 1 | 1–16 | Threads for GenomicsDBImport reading | | --filter-expression | — | JEXL expression | Hard filter expression (e.g., "QD < 2.0") | | --java-options | — | -Xmx4g | Java heap size; use -Xmx16g for large genomes |

Common Recipes

Recipe 1: Single-Sample Variant Calling (No Cohort)

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

Recipe 2: Run Full Pipeline with 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}"

Expected Outputs

| Output | Format | Description | |--------|--------|-------------| | *.g.vcf.gz | GVCF | Per-sample GVCF with reference confidence blocks; input to GenomicsDBImport | | cohort_all.vcf.gz | VCF | Joint-genotyped multi-sample VCF; unfiltered | | cohort_filtered.vcf.gz | VCF | Filtered VCF; FILTER=PASS for passing variants | | *.recal.table | BQSR | Base quality recalibration table from BaseRecalibrator | | *.bqsr.bam | BAM | Recalibrated BAM; use as HaplotypeCaller input | | genomicsdb/ | Directory | GenomicsDB workspace per chromosome for joint genotyping |

Troubleshooting

| Problem | Cause | Solution | |---------|-------|----------| | SAM/BAM file has no @RG header | Missing read group from BWA-MEM2 | Re-align with -R "@RG\tID:...\tSM:...\tPL:ILLUMINA" | | Java OutOfMemoryError | Insufficient heap size | Add --java-options "-Xmx16g" or more | | HaplotypeCaller very slow | Single-threaded HMM | Add --native-pair-hmm-threads 8; use interval lists to parallelize by chr | | Empty GVCF output | Wrong interval or no reads in region | Check samtools view -c sample.bam chrN for read counts | | VQSR fails (< 10k variants) | Too few variants for training | Use hard filters instead of VQSR for small cohorts or exomes | | GenomicsDB import fails on existing path | GenomicsDB workspace already exists | Delete existing workspace: rm -rf genomicsdb/chr1 before re-running | | IndexOutOfBoundsException | Chromosome name mismatch between BAM and reference | Ensure genome FASTA and BAM use same chr naming (chr1 vs 1) | | BCFtools/tabix index missing | Tabix index (.tbi) not created | Run gatk IndexFeatureFile -I file.vcf.gz or tabix -p vcf file.vcf.gz |

References

• GATK documentation — official pipeline guides and tool documentation
• GATK Best Practices for germline short variants — step-by-step protocol
• GATK GitHub: broadinstitute/gatk — source code, releases, and resource bundle
• Van der Auwera GA & O'Connor BD (2020) Genomics in the Cloud — O'Reilly Media; comprehensive GATK4 guide