Bio Gatk Variant Calling

Version Compatibility

Reference examples tested with: GATK 4.5+, bcftools 1.19+

Before using code patterns, verify installed versions match. If versions differ:

• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Note: flag DEFAULTS (e.g. --max-alternate-alleles, --heterozygosity, --standard-min-confidence-threshold-for-calling) drift across GATK 4.x releases and DRAGEN-mode ports models into the caller. Confirm any default with gatk <Tool> --help rather than trusting a memorized value.

GATK Variant Calling

"Call germline variants from my BAM with GATK" -> Detect SNPs/indels by locally reassembling candidate haplotypes in active regions and genotyping reads against them.

• CLI: gatk HaplotypeCaller (germline diploid), gatk Mutect2 (somatic / mitochondria / mosaic)

Why HaplotypeCaller, not a pileup caller

A pileup/position-based caller (bcftools mpileup, the old UnifiedGenotyper) genotypes each column of the alignment independently, trusting the mapper's per-read placement. Near indels and clustered variation that placement is a per-read greedy optimum, not a locus-consistent one, so the same indel gets (mis)placed differently in different reads and is systematically misrepresented. HaplotypeCaller (HC) discards the local alignment and re-derives it: in any region showing signal it performs local de-novo reassembly of candidate haplotypes, then realigns every read against those haplotypes. The indel is then represented ONCE, in an assembled haplotype, not independently in each read. This is why HC (and DRAGEN, and DeepVariant) beat pileup callers on indels and complex loci, and why HC is the reference implementation everything else is benchmarked against.

How HaplotypeCaller works (decision-relevant mechanism)

1. Active-region determination. HC computes a per-locus activity score (a fast reference-vs-non-reference genotype-likelihood contrast on the pileup), smooths it with a Gaussian kernel, and thresholds it (--active-probability-threshold, default 0.002) to seed active regions, padded so the assembler sees flanking reference. Non-active bases still get a reference-confidence emission in GVCF mode. Pathologically high-depth or highly repetitive regions can fail to assemble -- exactly where DeepVariant/DRAGEN pull ahead.
2. Local reassembly. For each active region HC builds a de-Bruijn-like graph from k-mers of the reference plus overlapping reads (default k = 10 and 25, merged), prunes low-weight (error) edges, and enumerates best-supported haplotypes. Each haplotype is Smith-Waterman-realigned to the reference to translate assembled sequence into concrete SNP/indel events.
3. PairHMM read-vs-haplotype likelihoods. For every (read, haplotype) pair HC computes P(read | haplotype) with a Pair Hidden Markov Model that integrates over ALL alignments (forward algorithm), not just the best one -- the statistically correct way to weight support when the alignment itself is uncertain. This is the dominant compute cost; production runs use vectorized (AVX/AVX-512) kernels via --pair-hmm-implementation and --native-pair-hmm-threads. DRAGEN moves this same kernel onto an FPGA.
4. Genotype likelihoods. Per-haplotype likelihoods are marginalized to per-allele likelihoods, then genotype likelihoods (PLs) are computed under the assumed ploidy via the Bayesian DePristo/GATK model, emitting GT/AD/DP/GQ/PL plus site annotations (QD, FS, MQ, MQRankSum, ReadPosRankSum, SOR).
5. The diploid assumption. Genotyping defaults to diploid (--sample-ploidy 2), which hard-codes true alleles at fraction 0, 0.5, or 1.0. Anything violating that (pooled, polyploid, CNV, mosaic, hemizygous chrX/Y) needs an explicit --sample-ploidy or a somatic caller -- see Edge Cases.

Pipeline decision tree

What is the analysis context?
├── Single sample, want DRAGEN-like accuracy, open-source -> HaplotypeCaller --dragen-mode (hard-filter on QUAL)
├── Cohort < ~2000, human -> per-sample -ERC GVCF -> joint genotype -> VQSR/VETS (see joint-calling)
├── Cohort > ~2000, human -> ReblockGVCF + GnarlyGenotyper "Biggest Practices", or DeepVariant + GLnexus (see joint-calling)
├── Non-human / non-model organism -> hard filtering (no VQSR training resources)
├── Targeted panel / small exome -> hard filtering (too few variants for VQSR)
├── Non-diploid / pooled / sex chromosomes / mitochondria -> set --sample-ploidy or use Mutect2 (see Edge Cases)
└── Somatic / mosaic variants -> Mutect2 (not HaplotypeCaller)

Single-sample calling

Goal: Call germline SNPs and indels from one sample.

Approach: Run HaplotypeCaller directly to VCF; add annotations, intervals, or a calling-confidence floor as needed.

# Basic call
gatk HaplotypeCaller -R reference.fa -I sample.bam -O sample.vcf.gz

# Exome / panel: restrict to capture targets (much faster, fewer off-target artifacts)
gatk HaplotypeCaller -R reference.fa -I sample.bam -L targets.interval_list -O sample.vcf.gz

# Add standard annotations explicitly (usually emitted by default; force when a downstream filter needs them)
gatk HaplotypeCaller -R reference.fa -I sample.bam -O sample.vcf.gz \
    -A Coverage -A QualByDepth -A FisherStrand -A StrandOddsRatio \
    -A MappingQualityRankSumTest -A ReadPosRankSumTest

MarkDuplicates is required before calling for all modes. Whether BQSR precedes it is a real decision -- see below.

The GVCF reference-confidence model (-ERC GVCF)

Goal: Make per-sample calls that can later be combined into a cohort without re-visiting BAMs.

Approach: Emit a GVCF that records confidence at EVERY position, then joint-genotype separately.

-ERC GVCF records, at every position (variant and non-variant), how confidently the site is homozygous reference. Two properties make it the backbone of scalable cohort calling:

• The <NON_REF> symbolic allele. Every record carries a symbolic ALT <NON_REF> ("any allele not yet observed") with AD/PL computed against it. At joint-genotyping time a variant discovered in ANOTHER sample can therefore be evaluated in THIS sample even though this sample looked reference -- the evidence against <NON_REF> supplies it. This is what makes per-sample GVCFs forward-compatible with alleles found later in the cohort, and lets joint genotyping distinguish confident hom-ref from no-data (the "squaring-off" of a ragged genotype matrix).
• GQ banding. Contiguous non-variant sites with similar GQ collapse into homRef blocks so a GVCF is not one line per base. -ERC BP_RESOLUTION disables banding (one line per base, larger files).

# Per-sample GVCF (do this once per sample; reusable when the cohort grows -- the N+1 win)
gatk HaplotypeCaller -R reference.fa -I sample.bam -O sample.g.vcf.gz -ERC GVCF

# Single-sample genotyping straight from one GVCF
gatk GenotypeGVCFs -R reference.fa -V sample.g.vcf.gz -O sample.vcf.gz

The N+1 problem: naive joint calling re-processes all N samples whenever the cohort changes; the GVCF captures the expensive assembly/likelihood work once per sample, so adding sample N+1 only re-runs the cheap consolidation + GenotypeGVCFs. Cohort consolidation (GenomicsDBImport vs CombineGVCFs), scaling to tens of thousands (ReblockGVCF, GnarlyGenotyper), and joint-genotyping mechanics live in variant-calling/joint-calling -- not duplicated here.

BQSR, DRAGSTR, and DRAGEN-GATK mode

Does BQSR still earn its place? (honestly unsettled). BaseRecalibrator/ApplyBQSR builds an empirical error model over base-call covariates (reported quality, read group, cycle, sequence context) to correct systematic, instrument-specific miscalibration. On older continuous-quality 4-color instruments (HiSeq, MiSeq) it mattered. On modern 2-color chemistry (NovaSeq/NextSeq) qualities are emitted in only ~4 coarse bins, leaving little smooth structure to recalibrate; empirically the callset is largely unchanged with vs without BQSR, with the delta concentrated in borderline-GQ sites. Broad keeps BQSR in Best Practices for pipeline consistency; several large pipelines drop it on binned data. Treat it as caller/instrument-dependent, not mandatory -- there is no consensus universal recommendation.

gatk BaseRecalibrator -R reference.fa -I sample.bam \
    --known-sites dbsnp.vcf.gz --known-sites Mills_and_1000G_gold_standard.indels.vcf.gz \
    -O recal.table
gatk ApplyBQSR -R reference.fa -I sample.bam --bqsr-recal-file recal.table -O sample.recal.bam

Where the indel-accuracy lever actually moved: DRAGSTR. Indel errors scale with tandem-repeat context, so the bigger gain is STR-aware indel modeling, not base-quality recalibration. DRAGEN-GATK adds DRAGSTR: a per-sample auto-calibration (CalibrateDragstrModel -> --dragstr-params-path) that models a-priori indel error/variant probability as a function of STR period (repeat-unit length) and length (copy number), adjusting the PairHMM indel gap priors before genotyping.

--dragen-mode. DRAGEN is Illumina's FPGA-accelerated map-align-call engine (a genome in ~20-25 min; wins the difficult-to-map precisionFDA V2 regions using alt-aware mapping). Illumina and Broad co-developed DRAGEN-GATK, porting DRAGEN's error models into open-source GATK to a functionally equivalent pipeline (Regier et al. 2018: pipelines are functionally equivalent when their call differences are smaller than sequencing-replicate differences). HaplotypeCaller --dragen-mode enables DRAGSTR plus BQD (Base Quality Dropout, systematic local quality collapse) and FRD (Foreign Read Detection, contaminating/mismapped reads), and replaces classic BQSR -- error modeling moves inside the caller. QUAL is well-calibrated, so hard-filtering on QUAL is sufficient without VQSR.

# Single-sample DRAGEN mode (no separate BQSR); GVCF variant for cohorts
gatk HaplotypeCaller -R reference.fa -I sample.markdup.bam -O sample.g.vcf.gz -ERC GVCF --dragen-mode

# DRAGEN-mode hard filter: GATK-recommended QUAL floor for DRAGEN-GATK output
gatk VariantFiltration -R reference.fa -V sample.vcf.gz -O sample.filtered.vcf.gz \
    --filter-expression "QUAL < 10.4139" --filter-name "DRAGENHardQUAL"  # Broad's documented DRAGEN-mode QUAL cutoff

DRAGEN-ML's advertised FP/FN reductions and speed figures are vendor-reported (Illumina), trained on GIAB truth; treat GIAB benchmark dominance with the overfitting caveat that DL/ML callers may not transfer identically to non-GIAB ancestries.

Allele-specific annotations (AS_)

Standard annotations (QD, FS, MQ, ...) lump all reads at a site together, so at a multiallelic site a real allele co-located with an error-driven allele shares one site-level pass/fail. AS_ annotations (request with -G AS_StandardAnnotation during GVCF calling / genotyping) compute each metric per allele (AS_QD, AS_FS, AS_SOR, AS_MQ, AS_MQRankSum, AS_ReadPosRankSum), letting AS_VQSR (-AS) filter each allele independently. The benefit grows with cohort size, because multiallelic sites -- where a true allele and an artifact collide at one position -- proliferate as sample count rises. The GVCF workflow propagates the raw per-allele data through joint genotyping to enable this.

Filtering: choose the method, then see filtering-best-practices

Filtering separates real variants from artifacts AFTER calling; SNPs and indels always filter separately (different annotation distributions). Pick the method here; the full VariantRecalibrator/ApplyVQSR, VETS, and hard-filter recipes with threshold rationale live in variant-calling/filtering-best-practices.

VQSR needs many variants overlapping the truth resources to fit a stable multivariate density; it is unreliable on single exomes or panels -- those must hard-filter. GATK is deprecating VQSR toward VETS (isolation-forest backend); verify the current recommended path in the GATK release notes before committing a pipeline.