The SNP Analysis Tools We Recommend

If you have your own whole-genome sequencing data and want to work with the single-nucleotide variants in it, here is the stack we would use. Use bwa-mem2 for alignment to GRCh38 with the analysis set FASTA and DeepVariant for small-variant calling. Use bcftools for every normalization and filtering step after that. Use hap.py against a Genome in a Bottle truth set to measure how good your calls are and Ensembl VEP for annotation. Use plink2 if you are working at the scale of hundreds of thousands of genotypes rather than one genome. Everything in that list is free and runs locally. The tools that show up in search results for “SNP analysis tools online” are mostly visualization front ends or legacy academic software, and none of them is the thing you want between your FASTQs and a conclusion.

Below is the reasoning for each choice, the parameters that matter, and the failure modes that will bite you.

What a SNP call is, and why tool choice matters

A single-nucleotide polymorphism is a position in the genome where a single base differs between individuals, and the density of these markers is what makes them the workhorse of human genetics ¹. In your own data, a SNP call is a probabilistic statement: given a pile of short reads that a mapper assigned to this locus, what genotype best explains the observed bases and their quality scores? Two things can go wrong. The reads can be in the wrong place, or the genotype model can misread a correct pile-up.

The second problem does not produce most of the error in a personal genome. The first problem does. That is why alignment and reference choice deserve more of your attention than the choice between two mature callers.

Reference genome: use GRCh38, and keep an eye on T2T

Use GRCh38 for anything you intend to interpret. It fixes a large number of errors present in GRCh37, adds alternate haplotype scaffolds, and measurably changes alignment and variant calls in regions that matter ². Almost every clinical database, gnomAD frequency table, and annotation cache you will want is keyed to GRCh38 coordinates now.

Specifically, download the “analysis set” FASTA, the one where the pseudoautosomal regions on chrY are hard-masked and the decoy contigs are included. Aligning to the plain primary assembly without decoys sends reads from unplaced repetitive sequence into the wrong autosomal locations and creates false heterozygous calls.

The T2T-CHM13 assembly is the better reference on technical grounds. It is gapless, and mapping to it reduces reference bias and improves alignment accuracy at disease-associated loci relative to GRCh38 ³. We would not make it your primary coordinate system yet, because the annotation ecosystem lags. A reasonable compromise is to call on GRCh38 for interpretation and re-call on T2T-CHM13 when a specific region looks suspicious. Free web platforms now support all three coordinate systems side by side, which makes that kind of cross-check cheap ⁴.

Do not try to solve the reference question with liftOver. Coordinate conversion loses variants in exactly the regions where the two assemblies disagree, which are the regions you were worried about.

The calling pipeline, with the flags that matter

For 30x short-read WGS, this is the sequence we run.

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:1\tSM:sample\tPL:ILLUMINA' \
  GRCh38_full_analysis_set_plus_decoy_hla.fa r1.fq.gz r2.fq.gz \
  | samtools sort -@ 8 -o sample.bam -
samtools index sample.bam

-K 100000000 fixes the chunk size so results are deterministic regardless of thread count. -Y uses soft clipping for supplementary alignments, which downstream structural-variant tools expect.

Then call:

run_deepvariant --model_type=WGS \
  --ref=GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --reads=sample.bam --output_vcf=sample.vcf.gz \
  --output_gvcf=sample.g.vcf.gz --num_shards=32

We prefer DeepVariant over GATK HaplotypeCaller for a single human genome for two reasons. It needs no base quality score recalibration step, and its convolutional model on pile-up images gives lower SNV error rates on standard benchmarks without any hand-tuned hard filters. GATK is still the right answer if you are joint-calling a cohort and want the full VQSR or filtering machinery, or if you must match an existing pipeline exactly. bcftools mpileup | bcftools call remains the fastest option and is fine for genotyping known sites, but it loses ground on indels and in repetitive regions.

Whichever caller you use, normalize before you compare anything:

bcftools norm -m -any -f ref.fa -c w sample.vcf.gz \
  | bcftools norm -d exact -Oz -o sample.norm.vcf.gz

Left-alignment and multiallelic splitting are not optional. Two VCFs that describe the same variant differently will look like a disagreement to every downstream tool.

For joint calling across several samples, use GLnexus on the gVCFs rather than GATK’s GenomicsDBImport. It is faster and its default DeepVariant configuration is well tested.

Note that this whole pipeline is tuned for DNA. If you are calling variants from RNA sequencing, the aligner must be splice-aware and the tool comparison is a different exercise entirely ⁵. Copy-number variation is also a separate problem with separate tools, and tool agreement there is considerably worse than for SNVs ⁶.

Measure your own accuracy before you trust a single call

This is the step most people skip. Take a Genome in a Bottle sample, run your exact pipeline on it, and compare against the truth VCF with the stratification BEDs:

hap.py HG002_truth.vcf.gz sample.norm.vcf.gz \
  -f HG002_confident.bed -r ref.fa \
  --engine vcfeval --threads 16 -o happy_out

You will get precision and recall separately for SNVs and indels. On a good 30x short-read pipeline, SNV F1 in the high-confidence regions lands around 99.8 percent. That number is the ceiling on everything you do afterwards, and it does not apply uniformly. Recall in segmental duplications, homopolymers, and the HLA region is much worse. Run mosdepth to get per-region coverage and treat any locus below roughly 15x as unresolved rather than reference.

Also run VerifyBamID2 for contamination. A few percent of cross-sample contamination shows up as a wave of low-allele-fraction heterozygous calls, which is easy to mistake for somatic mosaicism.

Annotation: VEP, run offline

For turning a VCF into something interpretable, use Ensembl VEP with a local cache:

vep --cache --offline --assembly GRCh38 --fork 8 \
  --everything --vcf --compress_output bgzip \
  --plugin CADD,whole_genome_SNVs.tsv.gz \
  -i sample.norm.vcf.gz -o sample.vep.vcf.gz

VEP over SnpEff, because VEP’s transcript handling and plugin system are better maintained. The same is true of its canonical and MANE flags. Whichever you use, pick one transcript model and state it. Consequence calls differ between RefSeq and Ensembl transcripts often enough to change a conclusion.

Then filter on population frequency from gnomAD before you look at anything. Of the roughly four to five million variants in a personal genome, the overwhelming majority are common and uninformative for rare-disease-style reasoning. Population genetic analysis of amino acid changes shows that a substantial fraction of nonsynonymous variants are under negative selection, which is the statistical basis for the deleteriousness scores you will see in the VEP output ⁷. Those scores are priors over classes of variants. They are not statements about you.

One trap: at a small but real fraction of positions, the reference genome carries the minor allele, so the “variant” call is the common state in the population and your reference match is the rare one ⁸. Frequency filters keyed to the reference allele rather than to gnomAD allele frequencies will get these backwards.

Interpretation of any specific finding in a clinically relevant gene needs a clinical geneticist or a genetic counselor working from an accredited confirmatory test. Research-grade calls from your own pipeline are not a diagnosis and should not be treated as one.

Genotype-scale work and array data

If your data is a consumer array rather than sequencing, your tools change. Work in PLINK 2 binary format (.pgen/.pvar/.psam), which is far faster than the old .bed trio and handles dosages. Arrays type roughly 600,000 to 1,000,000 sites directly, so nearly everything you see in a polygenic score comes from imputation against a reference panel. Run imputation with Minimac4 or the TOPMed imputation server, and keep the per-variant R² values. Dropping anything below R² of 0.8 removes most imputation-driven noise.

Array genotypes are excellent for the purpose they were built for, which is association and population-level inference ⁹. Allele frequency spectra from dense SNP data support demographic inference and population size history ¹⁰. They are poor substrate for asking about a rare variant in a specific gene, because the site is probably not on the chip.

Questions people also ask

Should I use hg19 or hg38?

Use hg38 (GRCh38). It corrects real errors in the earlier assembly and changes variant calls in regions that matter for interpretation ². Keep hg19 only if you must interface with an old dataset, and convert coordinates for lookup only, never for re-analysis.

What methods can be used to detect SNPs?

Three families in practice: hybridization arrays, targeted sequencing of a locus or panel, and whole-exome or whole-genome sequencing. Arrays genotype a fixed set of known sites cheaply. Sequencing discovers novel variants and is the only option if you do not know in advance which position you care about.

How accurate is SNP testing?

For 30x short-read WGS with a modern pipeline, SNV precision and recall in high-confidence regions are both around 99.8 percent. That still leaves thousands of errors per genome. They are concentrated in repeats, homopolymers, and low-coverage regions. Accuracy outside the benchmark regions is substantially lower and is not well characterized.

What is the best free software for DNA sequence analysis?

For personal genome work, the tools are the ones below.

• samtools/bcftools
• bwa-mem2
• DeepVariant
• hap.py
• VEP
• plink2

All are free, actively maintained, and run on a single workstation. Use IGV for visual inspection of any call you plan to act on.

Where do free web tools fit?

They are useful for cross-referencing a shortlist of variants across coordinate systems and annotation sources without setting up infrastructure ⁴. We would not run primary analysis through a browser, both for reproducibility and because uploading your genome to a third party is a decision you should make deliberately.

Footnotes

1. I. C. Gray. Single nucleotide polymorphisms as tools in human genetics. Human Molecular Genetics, 2000. https://doi.org/10.1093/hmg/9.16.2403 ↩

2. Yan Guo, Yulin Dai, Hui Yu, et al. Improvements and impacts of GRCh38 human reference on high throughput sequencing data analysis. Genomics, 2017. https://doi.org/10.1016/j.ygeno.2017.01.005 ↩ ↩²

3. Ilaria Cherchi, Francesco Orlando, Orsetta Quaini, et al. T2T-CHM13 reference genome reduces mapping bias and enhances alignment accuracy at disease-associated variants. 2025. https://doi.org/10.64898/2025.12.17.694618 ↩

4. Noam Hadar, Vadim Dolgin, Katya Oustinov, et al. VARista: a free web platform for streamlined whole-genome variant analysis across T2T, hg38, and hg19. Human Genetics, 2024. https://doi.org/10.1007/s00439-024-02671-4 ↩ ↩²

5. S. Akila Parvathy Dharshini, Y.-H. Taguchi, M. Michael Gromiha. Identifying suitable tools for variant detection and differential gene expression using RNA-seq data. Genomics, 2020. https://doi.org/10.1016/j.ygeno.2019.12.011 ↩

6. Migle Gabrielaite, Mathias Husted Torp, Malthe Sebro Rasmussen, et al. A comparison of tools for copy-number variation detection in germline whole exome and whole genome sequencing data. 2021. https://doi.org/10.1101/2021.04.30.442110 ↩

7. Adam R. Boyko, Scott H. Williamson, Amit R. Indap, et al. Assessing the Evolutionary Impact of Amino Acid Mutations in the Human Genome. PLoS Genetics, 2008. https://doi.org/10.1371/journal.pgen.1000083 ↩

8. Mahmoud Koko, Mohammed O. E. Abdallah, Mutaz Amin, et al. Challenges imposed by minor reference alleles on the identification and reporting of clinical variants from exome data. BMC Genomics, 2018. https://doi.org/10.1186/s12864-018-4433-3 ↩

9. Petra Nowotny, Jennifer M Kwon, Alison M Goate. SNP analysis to dissect human traits. Current Opinion in Neurobiology, 2001. https://doi.org/10.1016/s0959-4388(00)00261-0 ↩

10. Xiaoming Liu, Yun-Xin Fu. Exploring population size changes using SNP frequency spectra. Nature Genetics, 2015. https://doi.org/10.1038/ng.3254 ↩