SNP Analysis Tools: What to Use on Your Own Genome

Oak
snp wgs variant-calling annotation bioinformatics
A laboratory machine pulls a glowing strand through a reading gate, flagging occasional colored beads into lit glass trays below.

If you have your own data, the short answer is: call variants with DeepVariant (or GATK HaplotypeCaller if you need a gVCF for joint calling), normalize the VCF with bcftools norm, annotate with Ensembl VEP running offline against a local cache, join against gnomAD and ClinVar with vcfanno, and do everything downstream in plink2, bcftools query, or a dataframe library. The web tools that rank for “SNP analysis tools” (SNPnexus, SNPscan, QualitySNPng) are annotation front-ends or legacy array viewers. They are fine for looking up a handful of positions. They are the wrong shape for a 4.5-million-variant personal VCF that you want to query repeatedly.

Below is the pipeline we would run, and the places where it quietly goes wrong.

What a SNP is, in terms of the file you hold

A SNP is a single base position where the population carries more than one allele. Your genome differs from the GRCh38 reference at roughly 4 to 5 million positions, the large majority of which are single-nucleotide substitutions and the rest short indels. The first genome-wide catalog of this variation described 1.42 million SNPs, about one every 1.9 kb of sequence, which was the state of the map before short-read sequencing made per-person genomes routine.1

In a VCF, a SNP is one line: CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SAMPLE. The genotype lives in the sample column as 0/1, 1/1, or 0/0, plus DP (depth), GQ (genotype quality), and AD (allele depth per allele). Almost every analysis you will do is a filter or a join on those fields. Learn bcftools query -f and you have covered 80% of the work:

bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%GT\t%DP\t%GQ]\n' sample.vcf.gz

Detection: how SNPs are called from raw data

There are two ways your genotypes were produced, and they have very different error profiles.

Genotyping arrays (Illumina Global Screening Array, and the chips behind consumer kits) hybridize your DNA to 500k to 1M fixed probes. They only report positions someone chose in advance. Call rates above 98% are normal, and per-genotype accuracy at common sites is high, but the arrays are poor at rare variants: a site with 1-in-10,000 population frequency will have almost no heterozygous training clusters, and a small cluster-calling error becomes a confident-looking false positive. This is why a rare “pathogenic” hit in a raw array file should be treated as a hypothesis and confirmed by sequencing in an accredited clinical laboratory, with a genetic counselor or physician interpreting the result.

Short-read whole-genome sequencing reads your DNA directly, typically 2×150 bp on Illumina reversible-terminator chemistry, which is the chemistry that made per-individual genomes practical at 30x depth.2 At 30x, a heterozygous site is covered by roughly 15 supporting reads, and the caller’s job is distinguishing that from a sequencing error stack or a mismapped paralog.

For SNV and small indel calling from an aligned CRAM, we use DeepVariant:

run_deepvariant \
  --model_type=WGS \
  --ref=GRCh38_no_alt.fa \
  --reads=sample.cram \
  --output_vcf=sample.vcf.gz \
  --output_gvcf=sample.g.vcf.gz \
  --num_shards=$(nproc)

It outperforms hand-tuned hard filters on indels in particular, and it has one knob instead of twenty. Use GATK HaplotypeCaller -ERC GVCF instead when you plan to joint-call across family members, which is the case where a gVCF earns its disk space. With GATK on a single sample you have no cohort for VQSR, so fall back to hard filters: QD < 2.0, FS > 60.0, MQ < 40.0, MQRankSum < -12.5, SOR > 3.0 for SNVs.

Sanity-check the output before you believe anything in it:

bcftools stats -F GRCh38_no_alt.fa sample.vcf.gz | grep -E '^SN|^TSTV'

Genome-wide Ti/Tv should land near 2.0 to 2.1. Exome-restricted, closer to 3.0. A Ti/Tv of 1.5 means your call set is contaminated with error and every downstream count is inflated.

If you want a real accuracy number rather than a vibe, run hap.py against a Genome in a Bottle truth set restricted to the high-confidence BED:

hap.py HG002_truth.vcf.gz sample.vcf.gz \
  -f HG002_confident.bed -r GRCh38_no_alt.fa \
  --engine=vcfeval -o bench

You cannot do this on your own genome (there is no truth set for you), but running it once on HG002 with your exact pipeline tells you what your F1 is on SNVs versus indels, and where it collapses. It collapses in low-complexity regions, homopolymers longer than about 10 bp, segmental duplications, and the HLA and killer-cell immunoglobulin-like receptor loci. Do not interpret calls there from short reads.

Normalization: the step everyone skips and then regrets

Two VCFs can represent the same variant differently. bcftools fixes this, and every join you do afterward depends on it:

bcftools norm -m -any -f GRCh38_no_alt.fa -Ou sample.vcf.gz \
  | bcftools norm -d exact -Oz -o sample.norm.vcf.gz
bcftools index -t sample.norm.vcf.gz

-m -any splits multi-allelic records into one ALT per line. -f left-aligns and trims indels against the reference. -d exact drops duplicates. After this, your variant key is chrom:pos:ref:alt and it is stable.

Do not key on rsID. rsIDs are merged, retired, and remapped between dbSNP builds, and a rsID can refer to a multi-allelic site where the allele you care about is not the one your file carries. Key on coordinates plus alleles, and carry the rsID along as a label only.

Also settle your reference build before anything else. GRCh37 and GRCh38 coordinates are not interchangeable, and CrossMap or Picard LiftoverVcf will silently drop or flip a small fraction of sites, disproportionately in the regions that were wrong in GRCh37 and therefore the regions you most want. If you have the raw reads, realign to GRCh38 rather than lifting.

Annotation: what the variant does

Run VEP offline. The web front-ends are rate-limited and you will want to re-run this many times.

vep -i sample.norm.vcf.gz -o annotated.vcf.gz \
  --cache --offline --assembly GRCh38 --vcf --compress_output bgzip \
  --everything --pick_allele_gene \
  --plugin CADD,whole_genome_SNVs.tsv.gz \
  --plugin SpliceAI,snv=spliceai_scores.snv.hg38.vcf.gz

--pick_allele_gene collapses the transcript explosion to one consequence per allele per gene, which is what you want unless you are specifically studying isoforms. SnpEff and ANNOVAR are reasonable substitutes. VEP wins on plugin ecosystem and on matching Ensembl’s canonical transcript definitions.

Then attach population frequency and clinical assertions with vcfanno, which is fast and configured by a small TOML file:

[[annotation]]
file = "gnomad.genomes.v4.1.sites.chr*.vcf.bgz"
fields = ["AF", "AF_popmax", "nhomalt"]
ops = ["self", "self", "self"]
names = ["gnomad_af", "gnomad_af_popmax", "gnomad_nhomalt"]

Frequency is your strongest filter. A variant claimed to be severe and highly penetrant that appears at 2% allele frequency in gnomAD with hundreds of homozygotes is almost certainly misannotated. When you read ClinVar, read CLNREVSTAT too: a single-submitter assertion with no assertion criteria is not the same evidence class as a three-star expert panel review. Interpretation of a clinically relevant variant belongs with a clinician and a certified laboratory, not with your pipeline.

The interpretive bottleneck is real and it is not a tooling problem. Sequencing costs fell far faster than our ability to say what a given variant does, and the gap between “called accurately” and “understood” is where most personal-genome analysis stalls.3

Querying and population-scale work

plink2 is the right tool once you are doing anything statistical. Convert once:

plink2 --vcf sample.norm.vcf.gz --make-pgen --out sample

For polygenic scores, take a scoring file from the PGS Catalog and use --score:

plink2 --pfile sample --score PGS000123.txt 1 4 6 header cols=+scsums

Two failure modes here. First, strand ambiguity: A/T and C/G SNPs are palindromic, so if the scoring file and your VCF disagree on strand you will add a weight with the wrong sign and never notice. Filter them out or resolve strand by allele-frequency comparison. Second, portability: score weights derived in one ancestry group frequently lose much of their predictive value in another, and a raw score in isolation means nothing without the reference distribution it is meant to be compared against.

For repeated ad hoc queries over a full annotated VCF, we convert to Parquet and use DuckDB. Sub-second scans over 5 million rows beat re-parsing a bgzipped VCF every time.

Questions people also ask

What methods can be used to detect SNPs? Hybridization arrays, targeted PCR assays (TaqMan, KASP), Sanger sequencing for single sites, and short- or long-read sequencing with a variant caller. Arrays and targeted assays only see positions chosen in advance. Sequencing sees everything the reads cover, at the cost of more compute and a harder error model.

How can SNPs be detected accurately? Sufficient depth (30x for germline WGS), a caller validated against Genome in a Bottle truth sets, normalization before any comparison, and restriction to high-confidence regions. Report accuracy as SNV and indel F1 separately, because indel accuracy is always worse.

Can a SNP be inherited? Yes. Germline SNPs are inherited from your parents and are present in essentially every cell. De novo variants arise in the germline of a parent and appear in the child but neither parent, at a rate of roughly 50 to 70 per genome per generation. Somatic variants arise in a tissue during life and are not inherited.

Are there any SNPs linked to autism? Common SNPs associated with autism have been identified in large genome-wide association studies, and each contributes a very small amount of risk. Rare de novo variants and copy-number changes carry larger individual effects. No SNP or combination of SNPs in a consumer or research file diagnoses autism, and any question in this area belongs with a clinician.

How are SNPs analyzed once called? Filter by quality and depth, annotate consequence and population frequency, then either look at specific loci or aggregate across many with a scoring or burden method. Almost all of it reduces to joins on a normalized chrom:pos:ref:alt key.

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. David Bentley, Shankar Balasubramanian, Harold Swerdlow, et al. Accurate whole human genome sequencing using reversible terminator chemistry. Nature, 2008. https://doi.org/10.1038/nature07517 ↩

  3. 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 ↩