Running a Variant Effect Predictor on Your Own Whole Genome

By the end of this you will have a single tab-delimited file (and a Parquet copy) where every alternate allele in your genome carries a predicted molecular consequence against a specific transcript set, an allele frequency from gnomAD, and scores from a handful of predictors that disagree with each other in informative ways. You need: a germline VCF from a whole-genome sequence (30x or better, ideally from a pipeline built on GATK or DeepVariant), the exact reference FASTA that VCF was called against, roughly 400 GB of free disk if you want the full plugin data set, 16 GB of RAM, and Docker or a working Perl toolchain. Budget an evening for downloads and about 30-90 minutes of wall clock for the annotation run itself on 4-8 cores. Everything here is measurement and interpretation. Annotation is not diagnosis, and a variant that looks alarming in a TSV needs a clinical laboratory and a genetics clinician before it means anything about your health.

1. Normalize the VCF before you annotate anything

VEP’s consequence calls depend on the exact representation of an indel. AGCT -> A and GCT -> are the same deletion written two ways, and multi-allelic records hide alleles from some downstream tools. Split and left-align first:

bcftools norm -m -any -f GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --check-ref ws \
  -Oz -o me.norm.vcf.gz me.vcf.gz
bcftools index -t me.norm.vcf.gz
bcftools stats me.norm.vcf.gz | head -30

--check-ref ws warns on reference mismatches and swaps REF/ALT where the ALT matches the reference. If you get thousands of warnings, your FASTA is the wrong build or the wrong analysis set (the difference between the “no-alt” analysis set and the full GRCh38 primary assembly is enough to break things). Confirm the build from the VCF header contig lengths: chr1 is 248,956,422 bp on GRCh38 and 249,250,621 on GRCh37. Do not guess.

Filter to PASS if your caller emits filtered records, but keep a copy of everything. Hard filters in a personal genome throw away real variants, and you will want to go back.

2. Install VEP and pull a cache

Two viable paths. Docker is the one we use, because VEP has a long Perl dependency chain (Bio::DB::HTS, DBD::mysql, Bio::EnsEMBL::XS) that breaks in interesting ways on macOS.

docker pull ensemblorg/ensembl-vep:release_113.0
mkdir -p $HOME/vep_data && chmod a+rwx $HOME/vep_data

docker run -t -i -v $HOME/vep_data:/data ensemblorg/ensembl-vep \
  perl INSTALL.pl -a cf -s homo_sapiens -y GRCh38 -c /data --CACHE_VERSION 113

The native route, if you prefer it:

git clone https://github.com/Ensembl/ensembl-vep.git
cd ensembl-vep && perl INSTALL.pl -a acfp -s homo_sapiens -y GRCh38 -g all

-a acfp installs the API, cache, FASTA, and all available plugins. The -g all grabs every plugin’s code (not its data files, which you fetch yourself).

Three cache flavors exist for human. homo_sapiens is Ensembl/GENCODE transcripts. homo_sapiens_refseq is RefSeq. homo_sapiens_merged is both, with a SOURCE field telling you which database each transcript came from. Pick merged if you ever need to compare against a clinical report, which will usually quote a RefSeq accession like NM_000059.4. Pick the plain Ensembl cache if you want a smaller download and consistent transcript models. The GRCh38 merged cache is around 25-30 GB uncompressed.

On the question people search for: GENCODE and Ensembl are not separate annotation sets. The GENCODE consortium produces the gene annotation, and Ensembl ships it as the Ensembl gene set for human and mouse. GENCODE release 44 and Ensembl release 110 carry the same transcripts. RefSeq is a genuinely different annotation, maintained at NCBI, and it disagrees with GENCODE on exon boundaries and transcript existence often enough to change a consequence call.

3. Run a baseline annotation

Start without plugins so you know the pipeline works:

docker run -v $HOME/vep_data:/data -v $PWD:/work ensemblorg/ensembl-vep vep \
  --input_file /work/me.norm.vcf.gz \
  --output_file /work/me.vep.tsv \
  --tab --compress_output bgzip \
  --cache --dir_cache /data --merged --offline \
  --assembly GRCh38 --fasta /data/homo_sapiens/113_GRCh38/Homo_sapiens.GRCh38.dna.toplevel.fa.gz \
  --everything \
  --mane --canonical --biotype --symbol --hgvs --hgvsg \
  --fork 8 --buffer_size 25000 \
  --stats_file /work/me.vep.html

What matters here:

• --offline forbids any database lookup. Without it VEP will silently fall back to the Ensembl MySQL server and your run will take days. Always set it.
• --fasta is required for HGVS notation and for indel shifting. Point it at the same sequence you normalized against.
• --everything is shorthand for a bundle including --sift b --polyphen b --af --af_gnomade --af_gnomadg --max_af --pubmed --var_synonyms --variant_class and more. It roughly triples output width. Worth it on a personal genome, where you only run this a few times.
• --fork 8 is close to the sweet spot. Scaling above 12 forks gives diminishing returns because the cache reads become the bottleneck.
• --buffer_size trades memory for speed. 25000 uses several GB with forking. Drop to 5000 if you get OOM kills.

A 30x WGS VCF has roughly 4-5 million small variants. Expect 30-60 minutes at 8 forks with the merged cache, longer once plugins are attached.

4. Add the plugins that carry real information

The base VEP consequence is a categorical statement about transcript structure: missense, stop_gained, splice_donor, intron_variant. It says nothing about whether a missense change matters. That comes from plugins.

Download each data file and tabix-index it. The stack we use:

AlphaMissense. Precomputed pathogenicity scores for ~71 million possible human missense substitutions. Grab AlphaMissense_hg38.tsv.gz from Zenodo, then:

tabix -s 1 -b 2 -e 2 -f -S 1 AlphaMissense_hg38.tsv.gz

CADD v1.7. whole_genome_SNVs.tsv.gz is about 80 GB plus its index, and the InDels file adds several more. CADD is the only one here that scores non-coding variants genome-wide, which is its main value.

SpliceAI. Precomputed scores from Illumina, distributed under a license that requires registration. Two VCFs, SNV and indel, hg38. SpliceAI’s delta scores are the single most useful non-coding annotation for a personal genome, because splice-disrupting variants outside the canonical ±2 sites are common and VEP alone calls them intron_variant.

LOFTEE. Flags loss-of-function calls as high or low confidence. Roughly a third to a half of naive LoF calls in a healthy genome fail LOFTEE’s filters (last-exon truncations, rescue splice sites, ancestral alleles). On GRCh38 you need the grch38 branch of the plugin repo plus human_ancestor.fa.gz and the GERP bigwig.

dbNSFP 4.x, if you want 40+ missense predictors in one column set. Large, redundant, useful mainly when you want to see whether predictors agree.

Full command:

vep --input_file me.norm.vcf.gz --output_file me.vep.tsv --tab --compress_output bgzip \
  --cache --dir_cache /data --merged --offline --assembly GRCh38 \
  --fasta /data/Homo_sapiens.GRCh38.dna.toplevel.fa.gz \
  --everything --mane --hgvs --hgvsg --fork 8 --buffer_size 25000 \
  --plugin AlphaMissense,file=/data/AlphaMissense_hg38.tsv.gz \
  --plugin CADD,snv=/data/whole_genome_SNVs.tsv.gz,indels=/data/gnomad.genomes.r4.0.indel.tsv.gz \
  --plugin SpliceAI,snv=/data/spliceai_scores.raw.snv.hg38.vcf.gz,indel=/data/spliceai_scores.raw.indel.hg38.vcf.gz \
  --plugin LoF,loftee_path:/plugins/loftee,human_ancestor_fa:/data/human_ancestor.fa.gz,gerp_bigwig:/data/gerp_conservation_scores.homo_sapiens.GRCh38.bw \
  --plugin dbNSFP,/data/dbNSFP4.7a_grch38.gz,REVEL_score,MetaRNN_score,ClinPred_score,gnomAD_exomes_AF \
  --plugin UTRAnnotator,/plugins/uORF_5UTR_GRCh38_PUBLIC.txt \
  --dir_plugins /plugins

--dir_plugins must contain the .pm files. With Docker, mount your plugin directory to /plugins.

One thing to hold onto about these scores: most of them were trained or benchmarked on clinical variant databases, which are enriched for genes people already study. Systematic experimental measurement of missense effects (deep mutational scanning across many proteins) gives a cleaner yardstick, and predictors calibrated against that kind of functional evidence behave differently from ones tuned on ClinVar labels ¹. When AlphaMissense and REVEL disagree by a wide margin on a variant in a gene with little literature, treat the disagreement as the signal.

5. Choose one consequence per variant, deliberately

A single missense SNV in a gene with 15 transcripts produces 15 output rows by default. You need a policy.

• --pick gives one row per variant, chosen by VEP’s ordered criteria (MANE Select, then canonical, then biotype, then consequence severity, then transcript support level, then length).
• --pick_allele_gene gives one row per allele per gene. This is what we use. It keeps overlapping-gene cases visible without exploding the table.
• --flag_pick_allele keeps every row but adds a PICK=1 column. Largest output, most flexibility. Use this if you have the disk and intend to query with DuckDB anyway.

Add --mane_select filtering only if you want to force the clinical default transcript. MANE Select is the transcript RefSeq and Ensembl agree on, and roughly 19,000 protein-coding genes have one. Genes without a MANE Select are exactly the genes where transcript choice changes the answer.

6. Load it into something queryable

A 5-million-variant annotated TSV with --everything is 3-8 GB gzipped and around 100 columns. Do not try to open it in pandas without types.

import duckdb

duckdb.sql("""
  CREATE TABLE vep AS
  SELECT * FROM read_csv('me.vep.tsv.gz',
      delim='\t', header=true, skip=0,
      nullstr='-', sample_size=200000, all_varchar=true)
""")

duckdb.sql("COPY vep TO 'me.vep.parquet' (FORMAT PARQUET, COMPRESSION zstd)")

VEP’s tab output writes ## header lines before the column header, so strip them first with zcat me.vep.tsv.gz | grep -v '^##' | gzip > clean.tsv.gz. Read everything as varchar, then cast the numeric columns explicitly, because - is used as the null marker and will poison type inference.

Then the queries you want:

SELECT SYMBOL, Consequence, HGVSp, am_pathogenicity, CADD_PHRED,
       SpliceAI_pred_DS_AG, SpliceAI_pred_DS_DL, MAX_AF, LoF, LoF_filter
FROM vep
WHERE TRY_CAST(MAX_AF AS DOUBLE) < 0.001
  AND (TRY_CAST(am_pathogenicity AS DOUBLE) > 0.9
       OR LoF = 'HC'
       OR GREATEST(TRY_CAST(SpliceAI_pred_DS_AG AS DOUBLE),
                   TRY_CAST(SpliceAI_pred_DS_DL AS DOUBLE)) > 0.5)
ORDER BY CADD_PHRED DESC;

In a typical genome that returns a few hundred rows. Most are benign. This is where you stop doing computational biology and start needing a clinical lab, because turning a row into a statement about your health requires ACMG/AMP classification, family history, phenotype, and orthogonal confirmation by an accredited method.

VEP also ships filter_vep, which is convenient for quick passes without SQL:

filter_vep -i me.vep.tsv.gz -o rare_lof.tsv \
  --filter "MAX_AF < 0.001 and LoF is HC and BIOTYPE is protein_coding"

7. Know what the annotation does not cover

Roughly 1-2% of your genome is protein-coding exon. Everything else gets a consequence term like intergenic_variant or regulatory_region_variant that carries almost no predictive weight. The ENCODE pilot established early that functional elements extend far beyond coding sequence, including pervasive transcription and many regulatory regions with no protein product ². Long non-coding RNA annotation in particular remains incomplete, with catalogs differing substantially in which transcripts they include, so a variant’s consequence against a lncRNA depends heavily on which annotation release you ran ³.

Two partial mitigations. First, sequence-based constraint models trained across genomes now give genome-wide functional-constraint estimates that cover non-coding positions, which is a different kind of evidence from CADD’s supervised approach ⁴. Second, SpliceAI covers the one non-coding mechanism with a clean readout. Beyond that, be skeptical of any confident non-coding claim.

A related limitation on the coding side: a high-impact LoF call is only as meaningful as the evidence that the gene makes a protein at all. Community efforts to establish protein-level evidence for every predicted human protein-coding gene are still finishing the long tail ⁵. A stop-gained in a gene with thin protein evidence deserves less weight than one in a well-characterized gene.

8. Cross-check against your expression data

If you have RNA-seq from the same person, join on gene. A predicted high-impact LoF in a gene with zero TPM in every tissue you sequenced is a different object from one in a gene expressed at 200 TPM. Nonsense-mediated decay is directly observable: for a heterozygous stop-gained, check allele-specific expression at that position in the RNA BAM. If the ALT allele is present in DNA at ~50% and in RNA at ~10%, the transcript is being degraded, which supports the LoF call. This is the single highest-value cross-check available to someone with both layers, and it requires no new software beyond bcftools mpileup on the RNA alignment.

Common problems

Run takes hours and shows network activity. You forgot --offline. VEP is hitting the public Ensembl database.

ERR: Forked process failed. Almost always memory. Reduce --buffer_size to 5000 and --fork to 4. The per-fork memory scales with buffer size times number of plugins.

Cache version mismatch. --cache_version must match the installed cache directory name. A VEP 113 binary with a 110 cache will either fail or silently use stale transcripts. Check with ls $HOME/vep_data/homo_sapiens/.

Everything is intergenic_variant. Chromosome naming mismatch. Ensembl caches use 1, 2, X. Many VCFs use chr1. VEP handles this in most cases but not all. Force it: bcftools annotate --rename-chrs chr_map.txt.

Plugin returns empty columns. The data file is not tabix-indexed, or the index is older than the data file, or the build does not match (hg19 AlphaMissense against a GRCh38 VCF returns nothing rather than erroring). Test on a 1000-line VCF slice first, always.

RefSeq and Ensembl give different consequences. Expected. With the merged cache, check the SOURCE column and compare. If the two disagree on whether a variant is missense or synonymous, that is usually a difference in exon boundary or reading frame between the two annotations, and it is a real ambiguity, not a bug.

gnomAD frequency is missing for an indel. Representation mismatch. VEP matches on normalized position and allele, and gnomAD’s representation may differ if your normalization used a different reference. Re-run bcftools norm with the exact gnomAD reference.

HGVS strings missing. You did not pass --fasta, or the FASTA is not bgzipped with a .fai and .gzi index.

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. Tayfun Özçelık, Barış Kayaalp, Kerem Çil, et al. Prediction of human missense variant effects from functional evidence. Research Square, 2025. https://doi.org/10.21203/rs.3.rs-7536763/v1 ↩

2. Ewan Birney, Paul Flicek, Damian Keefe, et al. Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project. Nature, 2007. https://doi.org/10.1038/nature05874 ↩

3. Barbara Uszczyńska-Ratajczak, Julien Lagarde, Adam Frankish, et al. Towards a complete map of the human long non-coding RNA transcriptome. Nature Reviews Genetics, 2018. https://doi.org/10.1038/s41576-018-0017-y ↩

4. Chengzhong Ye, Gonzalo Benegas, Carlos Albors, et al. Predicting genome-wide functional constraints with GPN-Star. Nature, 2026. https://doi.org/10.1038/s41586-026-11005-5 ↩

5. Gilbert S. Omenn, Sandra Orchard, Lydie Lane, et al. The 2024 Report on the Human Proteome from the HUPO Human Proteome Project. Journal of Proteome Research, 2024. https://doi.org/10.1021/acs.jproteome.4c00776 ↩