DNA SNP Analysis: What a SNP Is and How to Analyze Your Own

A SNP (single nucleotide polymorphism) is one position in the genome where a single base differs between people: you might carry a C where the reference sequence has a T. “DNA SNP analysis” means determining which base you carry at such positions (genotyping), then asking what is known about the alleles you carry. In practice that is a file problem: sequencing reads go in, a VCF of genotypes comes out, and almost all the interesting work happens after that, in filtering and annotation. The first large public catalog of these sites mapped 1.42 million SNPs across the genome and estimated roughly one SNP per 1.9 kb of sequence.¹ Current dbSNP builds catalog on the order of a billion RefSNP records, most of them rare.

What counts as a SNP, precisely

Strictly, a SNP is a single-base substitution segregating in a population, conventionally at ≥1% minor allele frequency. A single-base change seen once in one person is a SNV (single nucleotide variant). In tooling nobody is careful about this: GATK, DeepVariant and bcftools all emit SNVs, and the VCF does not know or care about allele frequency until you annotate it.

A SNP has an identifier (rsID, e.g. rs1801133), a position on a specific reference build (chr1:11796321 on GRCh38 is not the same coordinate as on GRCh37), a reference allele, and one or more alternate alleles. Your genotype at that site is 0/0 (homozygous reference), 0/1 (heterozygous), or 1/1 (homozygous alternate). Three fields, and most confusion in consumer genetics comes from getting one of them wrong.

Your genome carries roughly 4–5 million variant sites relative to GRCh38, of which about 3.5–4 million are SNVs and the rest indels and larger events.²

Array genotypes versus whole-genome calls

A consumer array (23andMe, AncestryDNA) reports 600,000 to 1 million preselected positions. It measures only sites the chip designer chose, it cannot report a site that was not printed on the chip, and it produces no-calls at a few percent of probes. The output is a tab-separated text file: rsid chromosome position genotype, with genotypes like AG.

Whole-genome sequencing at 30x calls essentially every accessible base, including the ~99.9% of SNPs the array skipped and the rare variants that are far more likely to be functionally consequential. The cost is that you now own the calling problem, including its errors.³

If you have both, reconcile them. Array files report alleles on the probe’s design strand, which for Illumina chips is TOP/BOT orientation. That is not necessarily the forward genomic strand. An A/G on the chip can be a T/C in your VCF at the same locus. Compare on the plus strand after normalization or you will generate phantom discordance.

The pipeline we would run

Start from FASTQ or an unaligned BAM. Do not start from someone else’s VCF.

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

Use a reference with decoys and alt contigs (the 1000 Genomes analysis set). Decoys absorb reads from sequences absent in the primary assembly and cut false heterozygous calls in paralogous regions.

For calling, we use DeepVariant over GATK HaplotypeCaller for single-sample germline work:

run_deepvariant --model_type=WGS --ref=GRCh38.fa \
  --reads=s1.bam --output_vcf=s1.vcf.gz --num_shards=32

DeepVariant reaches roughly 99.8% F1 on SNVs against GIAB truth sets at 30x, needs no VQSR training, and does not require a cohort. GATK is the better choice if you are joint-calling a family, where gVCFs and GenotypeGVCFs let you distinguish homozygous reference from no-coverage.

Always normalize before you compare anything:

bcftools norm -m -any -f GRCh38.fa -Oz -o s1.norm.vcf.gz s1.vcf.gz
bcftools index s1.norm.vcf.gz

-m -any splits multiallelic records into one line per ALT, and -f left-aligns indels. Two VCFs of the same genome that disagree are usually just un-normalized. Most of the time they are not discordant.

Filtering

If you use GATK without VQSR, the standard hard filters for SNVs are a reasonable starting point. They are QD < 2.0, FS > 60.0, MQ < 40.0, MQRankSum < -12.5 and ReadPosRankSum < -8.0. On top of that, we drop genotypes with GQ < 20 or DP < 10 before drawing any conclusion about a single site.

Coverage is not uniform, and mean depth hides this. GC-rich promoters and first exons routinely run at a fraction of genome-wide mean depth. Library preparation choices, particularly the fragmentation method, materially change how evenly clinically relevant genes get covered.⁴ Before you trust a call, look at the actual depth at that position:

samtools depth -r chr1:11796321-11796321 -a s1.bam

And look at the reads. samtools tview or IGV takes thirty seconds and catches the majority of bad calls: all alternate reads on one strand, all at the same read position, or the site sitting inside a homopolymer run.

Annotation

Calls become interpretable only after annotation. We use Ensembl VEP:

vep -i s1.norm.vcf.gz -o s1.vep.vcf --vcf --cache --assembly GRCh38 \
    --everything --plugin CADD,whole_genome_SNVs.tsv.gz \
    --custom gnomad.genomes.v4.sites.vcf.gz,gnomAD,vcf,exact,0,AF,AF_popmax \
    --custom clinvar.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT

The fields that carry most of the weight:

• gnomAD allele frequency, now spanning roughly 800,000 individuals. A variant at 5% frequency in your ancestry group is not a rare deleterious allele, whatever a website says about it.
• ClinVar CLNSIG together with CLNREVSTAT. A “pathogenic” assertion with one submitter and no assertion criteria is weak evidence. Two or more submitters with no conflicts, or expert panel review, is much stronger.
• Consequence terms such as missense_variant, stop_gained, splice_donor_variant and synonymous_variant. Most of your SNPs are intergenic or intronic and will never be interpretable individually.
• CADD PHRED and, for anything near a splice junction, SpliceAI delta scores.

Where SNP analysis goes wrong

Reference build mismatch. An rsID looked up on GRCh37 coordinates and applied to a GRCh38 VCF gives a plausible-looking wrong answer. Check which build your file uses before anything else.

Low-complexity and segmentally duplicated regions. Short reads map ambiguously in roughly 5% of the genome, and error rates there are far above genome-wide averages.³ Regions like PMS2, SMN1/SMN2, CYP2D6, and the HLA locus need specialized callers or long reads. A short-read SNP call in CYP2D6 is not trustworthy on its own.

Strand and allele orientation, as above, mostly on array data.

Over-reading single variants. Most variants found in a healthy person’s genome are of uncertain significance, and penetrance estimates derived from families ascertained through affected individuals overstate risk when the same variant is found in someone with no family history.²

What you can conclude, and what needs a clinician

From a well-filtered, annotated VCF you can reliably determine your genotype at any well-covered common SNP and your carrier status for many recessive alleles. You can also determine your ancestry principal components and inputs to polygenic scores. Those scores transfer poorly across ancestry groups and should be read as population-level statements. They are not statements about an individual.

You cannot use a research-grade VCF as a clinical result. Any variant that appears to matter medically needs orthogonal confirmation in a CLIA/CAP laboratory and interpretation with a genetics professional. That professional will weigh your family history and phenotype as well as the file. The file alone is not enough. We treat our own sequencing as a hypothesis generator and nothing more.

Questions people also ask

What does “DNA SNP” mean? It is shorthand for a single nucleotide polymorphism: one position in the DNA sequence where different people carry different bases. Mapping these positions genome-wide was the goal of the first human variation maps, which found roughly one SNP every 1.9 kb.¹

How many SNPs do I have? Compared with the GRCh38 reference, about 3.5–4 million single-nucleotide differences, plus several hundred thousand indels.² Tens of thousands fall in protein-coding sequence, and a few hundred are rare protein-truncating variants.

Is a SNP the same as a mutation? They describe the same physical change with different framing. “SNP” implies the variant is common in the population. “Mutation” implies it arose recently or causes disease. The neutral term used in files and papers is “variant”.

Can I do SNP analysis on 23andMe raw data? Partly. You can check genotypes at the ~700,000 sites on the chip, but the file has no depth or quality information, carries a small no-call rate, and omits nearly every rare variant. For anything beyond common-variant lookups, sequence.

What is the best single file format to keep? Keep the CRAM (aligned reads, about 15–20 GB at 30x with a reference) and a normalized, bgzipped VCF with its index. From the CRAM you can re-call with a better caller in three years. From a VCF alone you cannot.

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. Tuuli Lappalainen, Alexandra J. Scott, Margot Brandt, et al. Genomic Analysis in the Age of Human Genome Sequencing. Cell, 2019. https://doi.org/10.1016/j.cell.2019.02.032 ↩ ↩² ↩³

3. Daniel C. Koboldt, Li Ding, Elaine R. Mardis, et al. Challenges of sequencing human genomes. Briefings in Bioinformatics, 2010. https://doi.org/10.1093/bib/bbq016 ↩ ↩²

4. Vanessa Process, Madana M.R. Ambavaram, Sameer Vasantgadkar, et al. Optimization of DNA Fragmentation Techniques to Maximize Coverage Uniformity of Clinically Relevant Genes Using Whole Genome Sequencing. Diagnostics, 2025. https://doi.org/10.3390/diagnostics15182294 ↩