SNP Calling Tools: Which One to Run on Your Own Genome

Oak
genomics variant-calling deepvariant gatk bioinformatics
A feathered reptile specimen on a velvet plinth in a dark studio, its paired flank scales mostly matched but some mismatched, with brass calipers and lenses hovering nearby.

If you have a 30x short-read whole genome of yourself and want a VCF you trust, run DeepVariant. It is a single container, it needs no base quality score recalibration, no variant quality score recalibration, and no truth-set priors, and it beats GATK HaplotypeCaller and bcftools on SNV and indel F1 at every depth we have tested. Use GATK if you need gVCFs that will be joint-called against a cohort later, or if you are matching a pipeline someone else already validated. Use bcftools mpileup when you want a five-minute answer about one locus. Use Clair3 or the DeepVariant ONT model if your reads came off a nanopore.

The rest of this is why, plus the commands and the places every caller quietly gets things wrong.

What the caller is deciding

A caller takes an aligned BAM or CRAM and, at each position, decides between homozygous reference, heterozygous, and homozygous alternate. The hard part is not the arithmetic on pileups. It is that read alignment is a guess, and in repetitive sequence the guess is systematically wrong in ways that look exactly like heterozygosity. A read from SMN2 mapped to SMN1 produces a clean-looking het call with balanced allele fractions and high mapping quality. No amount of Bayesian modeling on the pileup fixes that.

Expect roughly 4 to 5 million variants per human genome against GRCh38, about 3.5 million of them SNVs. The density is what the first genome-wide SNP maps found two decades ago, roughly one common SNP every couple of kilobases 1. The interesting variants for a single person are mostly in the rare tail, which is exactly the tail that deep population sequencing keeps expanding as more genomes get sequenced at 30x and above 2.

The recommendation: DeepVariant

DeepVariant converts each candidate site into a multi-channel pileup image and classifies it with a CNN. That sounds indirect and it is, but the practical consequence is that it learns the error modes of the sequencer and the aligner instead of you hand-specifying them as filter thresholds.

docker run -v "$PWD":/data google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/data/GRCh38_no_alt.fa \
  --reads=/data/sample.md.bam \
  --output_vcf=/data/sample.vcf.gz \
  --output_gvcf=/data/sample.g.vcf.gz \
  --num_shards=$(nproc) \
  --intermediate_results_dir=/data/tmp

On 32 vCPUs a 30x genome takes a few hours, most of it in make_examples. Memory is modest, under 2 GB per shard. Before this, align with bwa-mem2 mem -K 100000000 -Y or DRAGMAP, mark duplicates with samtools markdup -S, and stop. No indel realignment, no BQSR. BQSR costs you an hour of I/O and changes DeepVariant’s output by essentially nothing, because the model already treats reported base qualities as one noisy channel among several.

What you get: PASS calls with a GQ and a per-sample likelihood. Filter on FILTER="PASS" and, if you want a stricter set, GQ>=20. There are no QD/FS/SOR knobs to tune, which removes a large source of pipeline-to-pipeline irreproducibility.

Where GATK still earns its place

HaplotypeCaller’s local de novo assembly of active regions gives it good indel behavior, and the gVCF plus GenomicsDBImport plus GenotypeGVCFs path is still the cleanest way to joint-call a cohort with correct reference-confidence at non-variant sites. If you plan to add family members later, or your downstream tools expect GATK annotation fields, run it.

gatk --java-options "-Xmx16g" HaplotypeCaller \
  -R GRCh38_no_alt.fa -I sample.md.bam \
  -O sample.g.vcf.gz -ERC GVCF \
  --native-pair-hmm-threads 8

For a single sample you then need filtering. VQSR is unstable on one genome because it fits a Gaussian mixture on your own annotation distributions. Use hard filters instead: for SNVs, QD < 2.0 || FS > 60.0 || MQ < 40.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0. These are reasonable defaults, not tuned to your library, which is the point of the complaint against them.

bcftools mpileup | call -mv is fast and transparent and fine for genotyping a handful of known positions or sanity-checking a region by eye. It has no haplotype model, so its indels in homopolymers and short tandem repeats are meaningfully worse, and it should not be your primary genome-wide call set. Strelka2 sits in between: much faster than GATK, good germline accuracy, but development has slowed and the model-based approach has been overtaken.

Long reads change the answer

For PacBio HiFi, use DeepVariant with --model_type=PACBIO after aligning with pbmm2 or minimap2 -ax map-hifi. For nanopore R10.4.1 simplex data, either --model_type=ONT_R104 or Clair3 with the matching model (r1041_e82_400bps_sup_v4xx) works. The basecalling model matters more than the caller: super-accuracy basecalls versus fast basecalls move SNV F1 by a full percentage point or more, and a Clair3 model trained on a different chemistry will silently degrade.

The real argument for long reads is not per-base accuracy. It is that they phase. A 20x ONT run gives you megabase-scale haplotype blocks from whatshap phase or from the caller directly, and it calls into segmental duplications and the MHC where short reads produce garbage.

Benchmark before you believe anything

Do not accept a vendor’s F1. Run GIAB HG002 through your exact pipeline and compare.

hap.py HG002_benchmark.vcf.gz sample.vcf.gz \
  -f HG002_benchmark_noinconsistent.bed \
  -r GRCh38_no_alt.fa \
  --engine=vcfeval --engine-vcfeval-template GRCh38.sdf \
  --threads 16 -o happy_out

Use the vcfeval engine, not the default. It compares haplotype-equivalent representations, so a left-aligned indel and a right-aligned one are scored as the same call instead of one false positive plus one false negative. Stratify with --stratification and the GIAB stratification BEDs, because a single genome-wide F1 hides everything. Your numbers in low-mappability regions and long homopolymers will be far worse than your headline number.

One caveat on interpretation: the classic GIAB v4.2.1 benchmark deliberately excludes regions it cannot confidently resolve, so “99.9% F1” means 99.9% inside the easy 90-ish percent of the genome. Newer diploid assembly-based benchmarks built from telomere-to-telomere assemblies extend truth into segmental duplications and other previously excluded sequence, and accuracy measured against those is lower and more informative 3.

Failure modes worth knowing

  • Reference with ALT contigs and a non-ALT-aware aligner. Reads split across ALT and primary, mapping quality collapses, real variants disappear. Use the no-ALT analysis set or a properly ALT-aware bwa postprocessing step.
  • Unnormalized VCFs. Always run bcftools norm -f ref.fa -m -any -c x before comparing or annotating. Multi-allelic sites and right-aligned indels break naive joins.
  • Contamination. Check with VerifyBamID2. Above about 2% cross-sample contamination, heterozygous calls inflate and allele balance skews.
  • Coverage below ~20x. Het sensitivity drops fast because a 0.5 allele fraction at 12 reads has a real chance of producing 2 alt reads.
  • Paralog collapse. PMS2/PMS2CL, CYP2D6/CYP2D7, HBA1/HBA2, SMN1/SMN2. Short-read calls here are unreliable regardless of caller. Treat anything in these genes as unresolved until you have long reads or a targeted assay.

Finally, a called variant in a VCF is a measurement, not a finding. If something in your call set looks medically consequential, it needs orthogonal confirmation in a CLIA-certified laboratory and interpretation with a clinician or genetic counselor. Research-grade pipelines produce false positives at exactly the low-complexity, high-interest loci where you would most want to be right.

Questions people also ask

What is the single best SNP calling tool? For short-read human WGS, DeepVariant. It wins on accuracy without per-sample tuning, and the container makes the pipeline reproducible. GATK HaplotypeCaller is close behind and is the better choice when you need gVCFs for joint calling.

Do I need BQSR and indel realignment? No, for either DeepVariant or modern HaplotypeCaller. Indel realignment was removed from GATK4 because assembly-based calling subsumes it. BQSR has a small effect for GATK and effectively none for DeepVariant.

How do I compare two VCFs I produced? hap.py with the vcfeval engine, or rtg vcfeval directly. Never bcftools isec alone for accuracy work, because it compares representations rather than haplotypes and will report spurious differences on indels.

What coverage do I need? 30x mean for confident germline SNV and small-indel calling from short reads. At 15x you lose real heterozygous sites. Large population studies that define what counts as a novel variant are built on roughly this depth 2.

Can I call SNPs from RNA-seq? Yes, with splice-aware alignment (STAR two-pass) plus SplitNCigarReads and GATK, but only in expressed exons, and allele-specific expression skews genotypes toward the expressed allele. Use it to confirm expression of a variant you already found in DNA, not to discover genotypes.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. The International SNP Map Working Group, Cold Spring Harbor Laboratories:, Ravi Sachidanandam, et al. A map of human genome sequence variation containing 1.42 million single nucleotide polymorphisms. Nature, 2001. https://doi.org/10.1038/35057149 ↩

  2. Masao Nagasaki, Jun Yasuda, Fumiki Katsuoka, et al. Rare variant discovery by deep whole-genome sequencing of 1,070 Japanese individuals. Nature Communications, 2015. https://doi.org/10.1038/ncomms9018 ↩ ↩2

  3. Nancy F. Hansen, Nathan Dwarshuis, Hyun Joo Ji, et al. A complete diploid human genome benchmark for personalized genomics. Cell, 2026. https://doi.org/10.1016/j.cell.2026.06.016 ↩