How to Annotate the SNPs in Your Own Genome

By the end of this you will have a Parquet table with one row per variant per transcript, plus a short prioritized list you can read. The table has columns for consequence, gene and protein change. It also has gnomAD allele frequency split by ancestry group, ClinVar assertion and review status, and a set of in-silico scores. You need the following.

• a GRCh38 VCF from whole-genome sequencing (gVCF or joint-called single sample, both work)
• the matching reference FASTA
• about 250 GB of free disk
• 16–32 GB of RAM
• Docker
• a few hours, most of which is downloads

A single human genome differs from the reference at roughly four to five million sites, so nothing here is optional cleanup work. Annotation is the step that turns that pile into something with structure.

A worked example of what we are annotating: rs1815739 in ACTN3 is a C>T substitution that turns an arginine codon at position 577 into a stop codon, and the T allele is common worldwide. rs4988235 sits in an intron of MCM6 upstream of LCT and is the classic European lactase-persistence variant. rs671 in ALDH2 changes Glu504 to Lys and is common in East Asian populations. Three SNPs, three completely different annotation categories: nonsense, regulatory-by-linkage and missense. A good pipeline tells them apart without you knowing anything in advance.

1. Check what you have before annotating it

Reference build mismatches are the single most common way personal-genome annotation goes silently wrong. GRCh37 and GRCh38 coordinates differ for essentially every variant, and no tool will stop you from annotating hg19 positions against an hg38 cache.

bcftools view -h sample.vcf.gz | grep -E '^##(reference|contig=<ID=chr1,)'
bcftools stats sample.vcf.gz | grep -E '^SN|^TSTV'

Three things to read off. First, contig naming: chr1 (UCSC) or 1 (Ensembl). Pick one and stay with it, because your reference FASTA, gnomAD VCF, and ClinVar VCF must agree. Second, the SN block: a 30x WGS single sample should land around 4.0–5.0 million SNPs and 500,000–900,000 indels. If you see 600,000 total variants, that is exome or array data, not WGS, and half of this guide does not apply. Third, TSTV: genome-wide transition/transversion ratio should be near 2.0–2.1. Values under 1.8 mean the call set carries a lot of false positives, and every downstream filter you write will inherit them.

If you need to convert builds, use a liftover as a last resort and expect to lose a fraction of a percent of sites plus some silent strand flips. Re-aligning FASTQs to GRCh38 with bwa-mem2 or dragen-os is slower and correct.

2. Normalize and decompose

Annotation tools match variants to databases by position and allele string. Multi-allelic rows and non-left-aligned indels will not match ClinVar or gnomAD records that were written differently, and you get a false negative with no error message.

bcftools norm -m -any -f GRCh38.primary_assembly.genome.fa \
  -c s --threads 4 -Oz -o sample.norm.vcf.gz sample.vcf.gz
bcftools index -t sample.norm.vcf.gz

-m -any splits multi-allelics into one row per ALT. -f with -c s left-aligns indels and repairs REF alleles that disagree with the FASTA (it will warn loudly if many disagree, which again means build mismatch). Then drop the rows you will never look at:

bcftools view -e 'GT="0/0" || GT="./." || ALT="*"' \
  -Oz -o sample.clean.vcf.gz sample.norm.vcf.gz

If your file is a gVCF with reference blocks, this step takes it from tens of millions of records to about five million.

3. Run VEP offline with a cache

Our recommendation is Ensembl VEP running offline against a downloaded cache. We do not recommend a web service. Reasons: you keep the data local, runs are reproducible because the cache version is pinned, there is no upload size limit, and the plugin system lets you attach dbNSFP, SpliceAI, and custom VCFs in one pass. SnpEff is faster and has a simpler install, and it is a fine choice if all you want is consequence and gene. It does not have VEP’s plugin ecosystem. ANNOVAR’s licensing is awkward for anything shared. SNPnexus and other web tools are useful for looking up a handful of variants. They are not suited to five million of yours.

docker run --rm -v "$PWD":/data ensemblorg/ensembl-vep:release_113.0 \
  vep --offline --cache --dir_cache /data/vep_cache \
      --assembly GRCh38 --species homo_sapiens \
      --fasta /data/GRCh38.primary_assembly.genome.fa \
      --input_file /data/sample.clean.vcf.gz \
      --output_file /data/sample.vep.vcf.gz --vcf --compress_output bgzip \
      --fork 8 --buffer_size 20000 \
      --everything --hgvsg --canonical --mane --biotype --symbol \
      --numbers --total_length --shift_hgvs 1 \
      --plugin dbNSFP,/data/dbNSFP5.0a_grch38.gz,\
REVEL_score,CADD_phred,AlphaMissense_score,AlphaMissense_class,\
SIFT_pred,Polyphen2_HVAR_pred,GERP++_RS,clinvar_clnsig \
      --plugin SpliceAI,snv=/data/spliceai_scores.raw.snv.hg38.vcf.gz,\
indel=/data/spliceai_scores.raw.indel.hg38.vcf.gz \
      --stats_file /data/vep_summary.html

Notes on the flags that matter. --everything turns on a bundle including SIFT, PolyPhen, and the Ensembl regulatory feature overlap. Keep every transcript (the default) rather than using --pick: you will want to filter to MANE Select yourself later, and discarding alternates now means you cannot check whether a missense call on the canonical transcript is intronic on another. --shift_hgvs 1 makes HGVS notation 3’-shifted per convention, which matters for matching literature. --fork 8 with --buffer_size 20000 on a modern laptop gets a whole genome through in roughly 30–60 minutes. On memory: VEP with dbNSFP will use several GB per fork thread, so drop to --fork 4 on a 16 GB machine.

dbNSFP only covers nonsynonymous and splice-site SNVs, about 84 million possible substitutions. It will be empty for the ~98% of your variants that are intergenic or intronic, which is correct and expected.

4. Attach population frequencies

Frequency is the most useful single annotation you have. A variant present in 4% of any gnomAD population is not the explanation for a rare phenotype, regardless of what a prediction score says.

Download the gnomAD v4.1 joint genome+exome sites VCFs per chromosome, strip them to the INFO fields you need (the full files are enormous), and annotate:

bcftools annotate -a gnomad.joint.v4.1.sites.chr_all.trimmed.vcf.gz \
  -c 'INFO/AF:=INFO/AF_joint,INFO/AF_nfe:=INFO/AF_joint_nfe,\
INFO/AF_afr:=INFO/AF_joint_afr,INFO/AF_eas:=INFO/AF_joint_eas,\
INFO/AF_sas:=INFO/AF_joint_sas,INFO/AF_amr:=INFO/AF_joint_amr,\
INFO/nhomalt:=INFO/nhomalt_joint,INFO/faf95:=INFO/fafmax_faf95_max' \
  --threads 4 -Oz -o sample.vep.gnomad.vcf.gz sample.vep.vcf.gz

Use faf95 (the filtering allele frequency, the 95% confidence lower bound of the highest population frequency) rather than raw global AF for rarity filters. Raw AF pooled across populations hides the case where a variant is at 2% in one group and absent elsewhere.

Two cautions about frequency data. The reference panels themselves were built from ascertained samples, and allele-frequency spectra estimated from SNPs discovered in small discovery panels are systematically distorted toward common variants, an effect quantified early in the genome-wide polymorphism era ¹. And population labels are coarse: if your ancestry is not well represented in gnomAD, your “rare” variants will include a large number that are simply unsampled. The 1000 Genomes pilot made the scale of normal variation concrete, reporting that a typical individual carries on the order of 250 to 300 loss-of-function variants and 50 to 100 variants previously reported in association with inherited disease ². Almost none of those are causing anything. Frequency is how you find that out.

5. Add ClinVar, and read the review status

bcftools annotate -a clinvar_20260901.vcf.gz \
  -c 'INFO/CLNSIG,INFO/CLNREVSTAT,INFO/CLNDN,INFO/CLNSIGCONF,INFO/CLNVI' \
  --threads 4 -Oz -o sample.annotated.vcf.gz sample.vep.gnomad.vcf.gz
bcftools index -t sample.annotated.vcf.gz

CLNREVSTAT is the field people skip and should not. A criteria_provided,_single_submitter pathogenic call from one lab in 2014 is a different object than reviewed_by_expert_panel. Treat the star rating as a first-class filter: expert panel and multiple-submitters-no-conflict are worth attention, single submitter is a lead, no_assertion_criteria_provided is noise. Also keep CLNSIGCONF, which records conflicting interpretations that the top-level CLNSIG collapses.

The deeper point: pathogenicity assertions were mostly derived from patients ascertained because they were sick. Penetrance estimated in affected families overstates the risk to a person found to carry the same allele incidentally, and interpreting a sequenced genome in an unselected individual is a genuinely different problem from interpreting one in a symptomatic patient ³. If you find a ClinVar pathogenic variant in a clinically actionable gene, do not act on it from this pipeline. Take it to a clinical genetics service for confirmatory testing in an accredited lab. Research-grade WGS is not a diagnostic result.

6. Flatten to a table you can query

VEP writes all transcript consequences into one pipe-delimited CSQ INFO field. Explode it:

bcftools +split-vep sample.annotated.vcf.gz \
  -d -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t[%GT]\t[%DP]\t[%GQ]\t\
%INFO/AF\t%INFO/faf95\t%INFO/nhomalt\t%INFO/CLNSIG\t%INFO/CLNREVSTAT\t%INFO/CLNDN\t\
%SYMBOL\t%Gene\t%Feature\t%MANE_SELECT\t%CANONICAL\t%BIOTYPE\t%Consequence\t\
%IMPACT\t%EXON\t%HGVSc\t%HGVSp\t%REVEL_score\t%CADD_phred\t\
%AlphaMissense_score\t%AlphaMissense_class\t%SpliceAI_pred_DS_AG\t%SpliceAI_pred_DS_DL\n' \
  -A tab > variants.tsv

-d gives one row per transcript consequence. Then load into DuckDB, which handles a 30-million-row TSV on a laptop without complaint:

CREATE TABLE v AS
SELECT * FROM read_csv('variants.tsv', delim='\t', header=false,
  columns={'chrom':'VARCHAR','pos':'BIGINT','ref':'VARCHAR','alt':'VARCHAR',
           'qual':'DOUBLE','gt':'VARCHAR','dp':'INTEGER','gq':'INTEGER',
           'af':'DOUBLE','faf95':'DOUBLE','nhomalt':'INTEGER',
           'clnsig':'VARCHAR','clnrev':'VARCHAR','clndn':'VARCHAR',
           'symbol':'VARCHAR','gene':'VARCHAR','feature':'VARCHAR',
           'mane':'VARCHAR','canonical':'VARCHAR','biotype':'VARCHAR',
           'csq':'VARCHAR','impact':'VARCHAR','exon':'VARCHAR',
           'hgvsc':'VARCHAR','hgvsp':'VARCHAR','revel':'DOUBLE','cadd':'DOUBLE',
           'am_score':'DOUBLE','am_class':'VARCHAR','sai_ag':'DOUBLE','sai_dl':'DOUBLE'},
  nullstr='.');
COPY v TO 'variants.parquet' (FORMAT PARQUET, COMPRESSION ZSTD);

Parquet plus DuckDB is the right storage format here. Queries over five million variants return in under a second, and the file is a few hundred megabytes rather than tens of gigabytes.

7. Prioritize in tiers

Do not try to write one query that finds “the interesting variants.” Write several, each with an explicit definition, and look at all of them.

Tier 1, high-confidence clinical assertions:

SELECT symbol, hgvsc, hgvsp, clnsig, clnrev, clndn, gt, dp, gq, faf95
FROM v
WHERE clnsig LIKE '%athogenic%'
  AND clnsig NOT LIKE '%Conflicting%'
  AND clnrev IN ('reviewed_by_expert_panel',
                 'criteria_provided,_multiple_submitters,_no_conflicts')
  AND (faf95 IS NULL OR faf95 < 0.01)
  AND mane <> '' AND gq >= 30 AND dp >= 10
ORDER BY symbol;

Expect single digits to low tens of rows, mostly recessive carrier states.

Tier 2, rare predicted loss of function not yet in ClinVar:

SELECT symbol, csq, hgvsp, cadd, sai_ag, sai_dl, faf95, nhomalt, gt
FROM v
WHERE impact = 'HIGH' AND mane <> ''
  AND (faf95 IS NULL OR faf95 < 0.001)
  AND nhomalt < 5 AND gq >= 40
ORDER BY cadd DESC;

Tier 3, rare missense with concordant scores. Require agreement rather than any single tool: AlphaMissense class likely_pathogenic and REVEL above 0.7, say. Each predictor alone has a false-positive rate that makes a list of hundreds.

Tier 4, everything with a well-replicated common-variant annotation (the ACTN3, MCM6 and ALDH2 class). These are effect sizes on traits, not diagnoses, and they are the variants where public annotation quality varies the most.

One thing this scheme deliberately does not do: score single SNPs in isolation for anything where the functional unit is a haplotype. Human variation is strongly structured by linkage, with blocks of correlated SNPs inherited together ⁴. Pharmacogene star alleles such as CYP2D6 and CYP2C19 are defined by combinations of variants in phase, so run PharmCAT on the normalized VCF rather than reading individual rs IDs out of your table, and understand that unphased short-read data leaves some diplotype calls ambiguous. Anything downstream of that belongs in a conversation with a prescribing clinician.

Common problems

Everything annotates as intergenic. Your contig names do not match the cache. chr1 against an Ensembl cache silently produces no gene overlaps. Rename with bcftools annotate --rename-chrs.

ClinVar and gnomAD annotations are mostly empty on indels. You skipped normalization, or you normalized your sample but the annotation source is left-aligned differently. Run bcftools norm -f on both files with the same FASTA.

VEP reports a stop-gained on a transcript you have never heard of. You are looking at a non-MANE, often non-coding, transcript. Filter on mane <> '' for a first pass, then relax it for genes where MANE Select is not the tissue-relevant isoform.

Hundreds of “pathogenic” hits. You are almost certainly including one-star and no-assertion submissions, or not filtering on frequency. A true pathogenic variant at 3% allele frequency in gnomAD is a misclassification far more often than it is a discovery.

Genotypes disagree with a previous test. Compare depth and quality at the site. The call alone is not enough. Genotype confidence is a property of the platform, the depth, and the local sequence context, and individual genomes have been sequenced on technologies with very different error profiles ⁵. Homopolymers, segmental duplications, and the HLA and CYP2D6 regions are where short-read calls fail most. Anything you plan to act on gets orthogonal confirmation.

Your variant has no rs ID. That is fine and common. rs IDs come from dbSNP submissions, which grew out of systematic SNP discovery projects that sampled a limited number of chromosomes ⁶ and were deposited into the public domain by a consortium effort ⁷. Absence of an rs ID means nobody submitted that position. It does not mean the call is wrong.

Scores disagree with each other. REVEL, CADD, AlphaMissense, and SpliceAI are trained on different labels with different circularity problems, and several were trained partly on ClinVar itself, so agreement with a ClinVar label is not independent evidence. Use them to rank, never to conclude.

Footnotes

1. Andrew G. Clark, Melissa J. Hubisz, Carlos D. Bustamante, et al. Ascertainment bias in studies of human genome-wide polymorphism. Genome Research, 2005. https://doi.org/10.1101/gr.4107905 ↩

2. The 1000 Genomes Project Consortium. A map of human genome variation from population-scale sequencing. Nature, 2010. https://doi.org/10.1038/nature09534 ↩

3. A.J. Marian. Sequencing Your Genome: What Does it Mean?. Methodist DeBakey Cardiovascular Journal, 2014. https://doi.org/10.14797/mdcj-10-1-3 ↩

4. Benjamin A. Salisbury, Manish Pungliya, Julie Y. Choi, et al. SNP and haplotype variation in the human genome. Mutation Research - Fundamental and Molecular Mechanisms of Mutagenesis, 2003. https://doi.org/10.1016/s0027-5107(03)00014-9 ↩

5. Dmitry Pushkarev, Norma F Neff, Stephen R Quake. Single-molecule sequencing of an individual human genome. Nature Biotechnology, 2009. https://doi.org/10.1038/nbt.1561 ↩

6. David Altshuler, Victor J. Pollara, Chris R. Cowles, et al. An SNP map of the human genome generated by reduced representation shotgun sequencing. Nature, 2000. https://doi.org/10.1038/35035083 ↩

7. Arthur L. Holden. The SNP Consortium: Summary of a Private Consortium Effort to Develop an Applied Map of the Human Genome. BioTechniques, 2002. https://doi.org/10.2144/jun0203 ↩