domains / science / gatk-variant-calling / gatk-variant-calling

gatk-variant-calling

Skill card

combined

The two skills cover the same core GATK HaplotypeCaller workflow but diverge significantly in scope and depth. The first copy (FreedomIntelligence) is a procedural 'how-to' guide that includes specific commands for VQSR and Hard Filtering, which are often considered part of the 'calling' pipeline in older best practices. The second copy (GPTomics) is a conceptual and decision-oriented guide that explains the *mechanism* of HaplotypeCaller (local reassembly, PairHMM), provides a decision tree for choosing between tools (HC vs Mutect2 vs DRAGEN), and explicitly excludes filtering depth, pointing to other skills. The first copy lacks the mechanistic explanation and the decision tree for edge cases (ploidy, somatic), while the second copy lacks the specific VQSR/Hard Filter command examples that the first provides. Combining them yields a complete skill that covers both the 'why/when' (decision tree, mechanism) and the 'how' (commands for calling and filtering).

Repositories

Other copies

Combined skill

---
name: bio-gatk-variant-calling
description: Comprehensive GATK HaplotypeCaller workflow. Covers the local-reassembly mechanism, decision tree for tool selection (HC vs Mutect2 vs DRAGEN), single-sample and cohort GVCF workflows, and post-calling filtering (VQSR and Hard Filtering). Use when calling germline variants or deciding on the appropriate GATK pipeline.
tool_type: cli
primary_tool: gatk
---

# GATK Variant Calling

GATK HaplotypeCaller is the gold standard for germline variant calling. This skill covers the GATK Best Practices workflow, including mechanism, decision-making, and filtering.

## Prerequisites

BAM files should be preprocessed:
1. Mark duplicates
2. Base quality score recalibration (BQSR) - optional but recommended

## Why HaplotypeCaller?

HaplotypeCaller (HC) outperforms pileup callers on indels and complex loci because it performs **local de-novo reassembly** of candidate haplotypes in active regions, rather than genotyping each column independently. It uses a PairHMM to compute read-vs-haplotype likelihoods, integrating over all alignments. This makes it robust to mis-mapped reads near indels.

## Pipeline Decision Tree

Use this tree to select the correct tool and mode:

```
What is the analysis context?
├── Single sample, want DRAGEN-like accuracy, open-source -> HaplotypeCaller --dragen-mode
├── Cohort < ~2000, human -> per-sample -ERC GVCF -> joint genotype -> VQSR/VETS
├── Cohort > ~2000, human -> ReblockGVCF + GnarlyGenotyper or DeepVariant + GLnexus
├── 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
└── Somatic / mosaic variants -> Mutect2 (not HaplotypeCaller)
```

## Single-Sample Calling

### Basic HaplotypeCaller

```bash
gatk HaplotypeCaller \
    -R reference.fa \
    -I sample.bam \
    -O sample.vcf.gz
```

### With Standard Annotations

```bash
gatk HaplotypeCaller \
    -R reference.fa \
    -I sample.bam \
    -O sample.vcf.gz \
    -A Coverage \
    -A QualByDepth \
    -A FisherStrand \
    -A StrandOddsRatio \
    -A MappingQualityRankSumTest \
    -A ReadPosRankSumTest
```

### Target Intervals (Exome/Panel)

```bash
gatk HaplotypeCaller \
    -R reference.fa \
    -I sample.bam \
    -L targets.interval_list \
    -O sample.vcf.gz
```

## GVCF Workflow (Recommended for Cohorts)

The GVCF workflow enables joint genotyping across samples for better variant calls.

### Step 1: Generate GVCFs per Sample

```bash
gatk HaplotypeCaller \
    -R reference.fa \
    -I sample.bam \
    -O sample.g.vcf.gz \
    -ERC GVCF
```

### Step 2: Combine GVCFs

**Option A: GenomicsDBImport (Recommended for large cohorts)**

```bash
# Create sample map file: sample_map.txt
gatk GenomicsDBImport \
    --genomicsdb-workspace-path genomicsdb \
    --sample-name-map sample_map.txt \
    -L intervals.interval_list
```

**Option B: CombineGVCFs (Smaller cohorts)**

```bash
gatk CombineGVCFs \
    -R reference.fa \
    -V sample1.g.vcf.gz \
    -V sample2.g.vcf.gz \
    -O cohort.g.vcf.gz
```

### Step 3: Joint Genotyping

```bash
# From GenomicsDB
gatk GenotypeGVCFs \
    -R reference.fa \
    -V gendb://genomicsdb \
    -O cohort.vcf.gz

# From combined GVCF
gatk GenotypeGVCFs \
    -R reference.fa \
    -V cohort.g.vcf.gz \
    -O cohort.vcf.gz
```

## Post-Calling Filtering

### Variant Quality Score Recalibration (VQSR)

Machine learning-based filtering using known variant sites. Requires many variants (WGS preferred).

**SNP Recalibration:**

```bash
gatk VariantRecalibrator \
    -R reference.fa \
    -V cohort.vcf.gz \
    --resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap.vcf.gz \
    --resource:omni,known=false,training=true,truth=false,prior=12.0 omni.vcf.gz \
    --resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G.vcf.gz \
    --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp.vcf.gz \
    -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR \
    -mode SNP \
    -O snp.recal \
    --tranches-file snp.tranches

gatk ApplyVQSR \
    -R reference.fa \
    -V cohort.vcf.gz \
    -O cohort.snp_recal.vcf.gz \
    --recal-file snp.recal \
    --tranches-file snp.tranches \
    --truth-sensitivity-filter-level 99.5 \
    -mode SNP
```

**Indel Recalibration:**

```bash
gatk VariantRecalibrator \
    -R reference.fa \
    -V cohort.snp_recal.vcf.gz \
    --resource:mills,known=false,training=true,truth=true,prior=12.0 Mills.vcf.gz \
    --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp.vcf.gz \
    -an QD -an MQRankSum -an ReadPosRankSum -an FS -an SOR \
    -mode INDEL \
    --max-gaussians 4 \
    -O indel.recal \
    --tranches-file indel.tranches

gatk ApplyVQSR \
    -R reference.fa \
    -V cohort.snp_recal.vcf.gz \
    -O cohort.vqsr.vcf.gz \
    --recal-file indel.recal \
    --tranches-file indel.tranches \
    --truth-sensitivity-filter-level 99.0 \
    -mode INDEL
```

### Hard Filtering (When VQSR Not Suitable)

For small datasets, exomes, or single samples where VQSR fails.

**Extract SNPs and Indels:**

```bash
gatk SelectVariants -R reference.fa -V cohort.vcf.gz --select-type-to-include SNP -O snps.vcf.gz
gatk SelectVariants -R reference.fa -V cohort.vcf.gz --select-type-to-include INDEL -O indels.vcf.gz
```

**Apply Hard Filters:**

```bash
# Filter SNPs
gatk VariantFiltration \
    -R reference.fa \
    -V snps.vcf.gz \
    -O snps.filtered.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" \
    --filter-expression "ReadPosRankSum < -8.0" --filter-name "ReadPosRankSum-8" \
    --filter-expression "SOR > 3.0" --filter-name "SOR3"

# Filter Indels
gatk VariantFiltration \
    -R reference.fa \
    -V indels.vcf.gz \
    -O indels.filtered.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" \
    --filter-expression "SOR > 10.0" --filter-name "SOR10"
```

## Edge Cases

- **Non-diploid samples:** Use `--sample-ploidy` to adjust for polyploidy or hemizygosity.
- **Somatic/Mosaic:** Use `gatk Mutect2` instead of HaplotypeCaller.
- **Mitochondria:** Use `gatk Mutect2` with `--mitochondria-mode` or specific ploidy settings.
- **Sex Chromosomes:** Handle X/Y/PAR regions carefully, often requiring separate ploidy settings or specialized tools.