bio-gatk-variant-calling

What this skill does

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

# GATK Variant Calling

GATK HaplotypeCaller performs local de novo assembly of haplotypes in active regions, making it more accurate than pileup-based callers (bcftools) for indels and complex variants.

## Pipeline Decision Tree

```
What is the analysis context?
├── Single sample, highest accuracy -> DRAGEN-GATK mode (hard filter on QUAL)
├── Cohort <1000, human -> Standard Best Practices (GVCF -> joint genotype -> VQSR)
├── Cohort >1000, human -> GATK "Biggest Practices" or DeepVariant + GLnexus
├── Non-human organism -> Hard filtering (no VQSR training resources)
├── Targeted panel / small exome -> Hard filtering (too few variants for VQSR)
└── Somatic variants -> Mutect2 (not HaplotypeCaller)
```

## DRAGEN-GATK Mode (Current Recommendation)

DRAGEN-GATK mode is the recommended approach for single-sample germline calling. It incorporates three innovations that improve accuracy:
- **BQD (Base Quality Dropout)**: Corrects systematic base quality errors at read ends
- **FRD (Foreign Read Detection)**: Identifies contaminating reads at the genotyping level
- **DragSTR**: Improved short tandem repeat handling with a Markov model

Key differences from standard mode: BQSR is no longer needed (error correction happens during calling), and QUAL scores are more calibrated, making simple hard filtering on QUAL sufficient without VQSR.

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

For DRAGEN-GVCF mode (cohort workflows):
```bash
gatk HaplotypeCaller \
    -R reference.fa \
    -I sample.bam \
    -O sample.g.vcf.gz \
    -ERC GVCF \
    --dragen-mode
```

Note: DRAGEN mode is available only for single-sample calling. Joint genotyping with GenotypeGVCFs does not use DRAGEN mode but benefits from the improved per-sample GVCFs.

## Prerequisites

BAM files should be preprocessed:
1. Mark duplicates (always required)
2. Base quality score recalibration (BQSR) -- recommended for standard mode, NOT needed for DRAGEN-GATK mode

## Single-Sample Calling

**Goal:** Call germline SNPs and indels from a single sample using HaplotypeCaller.

**Approach:** Run local de novo assembly of haplotypes in active regions to detect variants with optional annotation enrichment.

**"Call variants from my BAM file using GATK"** -> Perform local haplotype assembly and genotyping on aligned reads using HaplotypeCaller.

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

### Adjust Calling Confidence

```bash
gatk HaplotypeCaller \
    -R reference.fa \
    -I sample.bam \
    -O sample.vcf.gz \
    --standard-min-confidence-threshold-for-calling 20
```

## GVCF Workflow (Recommended for Cohorts)

**Goal:** Enable joint genotyping across a cohort by generating per-sample genomic VCFs.

**Approach:** Call each sample in GVCF mode (-ERC GVCF), combine into a GenomicsDB or merged GVCF, then jointly genotype.

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 (GenomicsDBImport)

```bash
# Create sample map file
# sample_map.txt:
# sample1    /path/to/sample1.g.vcf.gz
# sample2    /path/to/sample2.g.vcf.gz

gatk GenomicsDBImport \
    --genomicsdb-workspace-path genomicsdb \
    --sample-name-map sample_map.txt \
    -L intervals.interval_list
```

### Alternative: CombineGVCFs (smaller cohorts)

```bash
gatk CombineGVCFs \
    -R reference.fa \
    -V sample1.g.vcf.gz \
    -V sample2.g.vcf.gz \
    -V sample3.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
```

## Variant Quality Score Recalibration (VQSR)

**Goal:** Apply machine learning-based variant filtering using known truth/training sets.

**Approach:** Build a Gaussian mixture model from annotation values at known sites, then apply a sensitivity threshold to classify variants.

### When to Use VQSR vs Hard Filtering

| Context | Filtering Method | Rationale |
|---------|-----------------|-----------|
| Human WGS cohort, >30 samples | VQSR | Enough variants for model training; truth sets available |
| Human WGS, single sample, DRAGEN mode | Hard filter on QUAL | DRAGEN QUAL scores are well-calibrated |
| Human exome, >30 samples | VQSR (with `--max-gaussians 4` for indels) | Fewer variants but usually sufficient |
| <30 exomes or gene panels | Hard filtering | Too few variants for reliable GMM training |
| Non-model organisms | Hard filtering | No HapMap/1000G training resources |
| Somatic calling (Mutect2) | FilterMutectCalls | Dedicated somatic filtering, not VQSR |

Sensitivity tranche selection: 99.5% for SNPs (0.5% true positives may be filtered) and 99.0% for indels (indels have higher error rates). Adjust based on the analysis goal -- use 99.7% for discovery, 99.0% for clinical stringency.

### SNP Recalibration

```bash
# Build SNP model
# Annotations: QD, MQ are most informative; FS/SOR detect strand bias; RankSum tests detect systematic ref/alt differences
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

# Apply SNP filter (99.5% sensitivity retains 99.5% of true positives)
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
```

### Allele-Specific VQSR (Recommended for Large Cohorts)

Standard VQSR evaluates all alleles at a site together. Allele-specific (AS) filtering evaluates each allele independently -- critical at multiallelic sites where a true variant and an artifact may co-occur. Increasingly important as cohort size grows (more multiallelic sites).

```bash
# Requires AS annotations from GenotypeGVCFs (-G AS_StandardAnnotation)
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 -AS \
    -O snp.as