Variant Calling

intermediate

Variant Calling Strategies

Germline vs Somatic Calling

🧬Germline (HaplotypeCaller)Single sample GVCF → joint genotyping → VQSR | Expected VAF 50% (het) or 100% (hom) | Rare disease, hereditary cancer

🎯Somatic (Mutect2)Tumour + Normal paired | Low VAF variants (1-50%) | Panel of normals | Oncology, cancer panels

🤖DeepVariant (ML)CNN on pileup images | >99.9% F1 on SNPs | Works on WES/WGS/FFPE | Alternative to GATK

GATK HaplotypeCaller - Best Practices

HaplotypeCaller uses a local de-novo assembly of the active region to call SNVs and indels simultaneously. Run in GVCF mode per sample, then joint-genotype across cohorts.

code

# Step 1: Per-sample GVCF
gatk HaplotypeCaller \
  -R hg38.fa \
  -I recal.bam \
  -O sample.g.vcf.gz \
  -ERC GVCF \
  -L intervals.interval_list \
  --dbsnp dbsnp_138.hg38.vcf.gz

# Step 2: Combine GVCFs (for cohort calling)
gatk CombineGVCFs \
  -R hg38.fa \
  -V sample1.g.vcf.gz \
  -V sample2.g.vcf.gz \
  -O cohort.g.vcf.gz

# Step 3: Joint genotyping
gatk GenotypeGVCFs \
  -R hg38.fa \
  -V cohort.g.vcf.gz \
  -O raw_variants.vcf.gz

Variant Filtration

VQSR (Variant Quality Score Recalibration) is ideal for large cohorts. For small cohorts/single samples, use hard filters.

code

# Hard filters (single sample / small cohort)
# SNP filters
gatk VariantFiltration \
  -V raw_variants.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" \
  -O filtered_snps.vcf.gz

# Indel filters
gatk VariantFiltration \
  --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"

QD (QualByDepth): variant quality normalised by depth; <2 = low confidence
FS (FisherStrand): strand bias; >60 (SNP) or >200 (indel) = likely artifact
MQ (RMSMappingQuality): mapping quality of supporting reads; <40 = poor alignment
MQRankSum: comparison of MQ between alt and ref reads; large negative = artifact
ReadPosRankSum: whether alt reads cluster at end of reads (artifact pattern)

DeepVariant (ML-based)

Google's DeepVariant uses a convolutional neural network trained on pileup images. Outperforms GATK on SNPs (F1 >99.9%) and is competitive on indels.

code

# Run with Docker (easiest)
docker run \
  -v /data:/data \
  google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WES \
  --ref=/data/hg38.fa \
  --reads=/data/sample.bam \
  --regions=/data/capture.bed \
  --output_vcf=/data/output.vcf.gz \
  --num_shards=16

Post-Alignment Processing Understanding VCF Files