Annotating Your Own VCF: A Working Pipeline

Oak
vcf annotation vep genomics bcftools
A feathered reptile-bird hybrid on a dark studio plinth, its spine feathers each clipped with small colored enamel tags.

By the end of this you will have two artifacts: an annotated, bgzipped VCF with per-transcript consequence records, and a flat Parquet table with one row per variant-transcript pair carrying gene symbol, consequence, MANE transcript ID, gnomAD allele frequency, ClinVar classification, and in-silico predictor scores. The table is a few hundred megabytes and queryable with DuckDB in under a second. You need: a GRCh38 VCF from short-read whole-genome sequencing (roughly 4–5 million variants for one person, mostly SNVs with 400k–600k indels), the matching reference FASTA, Docker or a conda environment, about 100 GB of free disk for the VEP cache and annotation databases, and 16 GB of RAM. Everything here is measurement and interpretation. Nothing in an annotated VCF is a diagnosis, and any variant you plan to act on needs confirmation in a clinical lab and a conversation with a genetic counselor or physician.

1. Read the VCF header before you touch anything

A VCF is a tab-delimited text file with two parts. Lines starting with ## are meta-information: the file format version, contig names and lengths, and the definitions of every INFO, FORMAT, and FILTER key used below. One line starting with #CHROM names the columns. Everything after is one variant record per line, with fixed columns CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO, then FORMAT and one column per sample. Annotation means writing new keys into the INFO column (and defining them in the header), which is why header hygiene matters: a tool that writes CSQ without a matching ##INFO=<ID=CSQ,...> line will break every downstream parser.

Start here:

bcftools view -h sample.vcf.gz | grep -E '^##(reference|contig=<ID=(chr)?1,|source|INFO=<ID=(AF|DP))'
bcftools stats sample.vcf.gz > stats.txt
grep -E '^SN' stats.txt

You are checking four things: which reference build, whether contigs are chr1 or 1, which caller produced the file, and the ratio of SNPs to indels. A transition/transversion ratio (bcftools stats reports it as TSTV) well below 1.9 for a whole genome usually means the file contains unfiltered or low-quality calls, and annotating those wastes hours and pollutes every downstream count.

2. Pin the reference build and contig naming

Almost every painful annotation bug is an assembly mismatch. GRCh37 and GRCh38 coordinates differ across most of the genome, and annotating GRCh37 positions against a GRCh38 cache produces plausible-looking nonsense: real gene names, real consequences, wrong variants. Check the build with a known position rather than trusting the header. rs334 (the HBB sickle variant) is at chr11:5,227,002 on GRCh38 and chr11:5,248,232 on GRCh37.

If contigs are 1, 2, X, rename them to match your reference:

bcftools annotate --rename-chrs chr_map.txt -Oz -o renamed.vcf.gz sample.vcf.gz
bcftools index -t renamed.vcf.gz

where chr_map.txt is two whitespace-separated columns (1 chr1, MT chrM, and so on). Do not lift over between builds unless you have no alternative. If you must, use CrossMap or Picard LiftoverVcf with the UCSC chain file, keep the reject file, and expect to lose 1–2% of variants, disproportionately in segmental duplications. We would rather re-call from the BAM/CRAM against GRCh38 than annotate a lifted file.

3. Normalize: split multiallelics, then left-align

This is the step people skip and then spend a week debugging. A record like chr1 100 . A AT,ATT is one line with two alleles. Databases like ClinVar and gnomAD store one allele per record, so a naive join misses both. Separately, an indel can be written at several equivalent positions in a repeat, and two files that disagree about which position to use will never match.

bcftools norm -m -any renamed.vcf.gz \
  | bcftools norm -f GRCh38.primary_assembly.genome.fa -c w -Oz -o norm.vcf.gz -
bcftools index -t norm.vcf.gz

-m -any splits multiallelic records into one ALT per line. -f with the reference does left-alignment and trimming. -c w warns (rather than silently continuing) when the REF allele in the file does not match the reference FASTA, which is your canary for a build mismatch. Run it and read the stderr summary: it reports how many records were realigned and how many REF mismatches it saw. Any REF mismatch count above a handful means stop and fix step 2.

Two caveats. Splitting multiallelics rewrites the GT field (1/2 becomes 0/1 on two lines), so downstream zygosity logic must read the split file, not the original. And if your caller emitted MNPs or complex alleles, bcftools norm -a will atomize them, which improves database matching but destroys the phase relationship between the component SNVs. For coding variants in the same codon that matters, because two adjacent substitutions can produce a different amino acid than either alone. Keep both versions of the file if you care about that.

4. Annotate with VEP against a local cache

Ensembl VEP is the tool we use. It resolves each variant against a transcript set and reports the consequence per transcript, which is the core of functional annotation: the same variant can be a missense change in one isoform, intronic in another, and upstream of a third 1. Download the merged cache (Ensembl plus RefSeq) and run fully offline. Online mode is 100x slower and will not finish a whole genome.

vep \
  --input_file norm.vcf.gz --format vcf \
  --output_file annotated.vcf.gz --vcf --compress_output bgzip \
  --cache --offline --dir_cache /data/vep --merged \
  --assembly GRCh38 --fasta GRCh38.primary_assembly.genome.fa \
  --species homo_sapiens \
  --fork 8 --buffer_size 20000 \
  --hgvsg --hgvsc --hgvsp --symbol --canonical --mane --biotype --numbers \
  --pick_allele_gene --pick_order mane_select,mane_plus_clinical,canonical,rank \
  --sift b --polyphen b \
  --plugin dbNSFP,/data/dbNSFP4.7a_grch38.gz,REVEL_score,CADD_phred,AlphaMissense_score,AlphaMissense_class,gnomAD_exomes_AF \
  --plugin SpliceAI,snv=/data/spliceai_scores.raw.snv.hg38.vcf.gz,indel=/data/spliceai_scores.raw.indel.hg38.vcf.gz \
  --stats_file vep_summary.html

Notes on the flags that matter:

  • --fork 8 with --buffer_size 20000 gets a whole genome down to roughly 30–60 minutes on a modern 8-core machine. Larger buffers increase memory roughly linearly.
  • --pick_allele_gene keeps one consequence per allele per gene instead of all transcripts. Without it, a genome produces 40M+ CSQ entries and your Parquet file balloons. With --pick_order starting at mane_select, you get the transcript clinical labs use. We recommend this over --everything, which turns on flags you did not ask for and slows the run substantially.
  • --mane adds the MANE Select and MANE Plus Clinical flags. Use these as your canonical reporting transcript. If you want the full per-transcript picture, run a second pass without --pick_allele_gene and write it to a separate file.
  • --hgvsg requires --fasta. Without it, HGVS genomic notation is silently omitted.

VEP writes everything into a single pipe-delimited CSQ INFO field whose column order is defined in the header. Do not parse it with cut. Use bcftools +split-vep, which reads the header to find field positions.

5. Add population frequency and clinical databases

VEP plugins cover predictor scores. For allele frequency and clinical assertions, annotating from indexed VCFs with bcftools annotate is faster and easier to keep current, since ClinVar ships a new GRCh38 VCF weekly.

# gnomAD v4 genomes, joint AF and popmax
bcftools annotate -a gnomad.v4.1.sites.chr_all.vcf.gz \
  -c INFO/gnomAD_AF:=INFO/AF,INFO/gnomAD_AF_grpmax:=INFO/AF_grpmax,INFO/gnomAD_nhomalt:=INFO/nhomalt \
  -Oz -o with_af.vcf.gz annotated.vcf.gz
bcftools index -t with_af.vcf.gz

# ClinVar
bcftools annotate -a clinvar_20260901.vcf.gz \
  -c INFO/CLNSIG,INFO/CLNREVSTAT,INFO/CLNDN,INFO/CLNHGVS \
  -Oz -o final.vcf.gz with_af.vcf.gz
bcftools index -t final.vcf.gz

bcftools annotate matches on CHROM, POS, REF, and ALT, which is exactly why step 3 was non-negotiable. The := syntax renames the source field so gnomAD’s AF does not collide with your caller’s AF. Silent collisions here are a real failure mode: bcftools will overwrite the existing key without complaint if you do not rename.

Sanity check the join rate. Roughly 99% of a healthy individual’s SNVs should have a gnomAD frequency. If you see 60%, your normalization did not match gnomAD’s (gnomAD v4 is normalized and split, so a mismatch means your file is not).

bcftools query -f '%INFO/gnomAD_AF\n' final.vcf.gz | awk '$1=="."{n++}END{print n, NR, n/NR}'

6. Flatten to Parquet and query it

The VCF is the archival artifact. For analysis you want columns.

bcftools +split-vep final.vcf.gz \
  -f '%CHROM\t%POS\t%REF\t%ALT\t%FILTER\t%INFO/gnomAD_AF\t%INFO/gnomAD_nhomalt\t%INFO/CLNSIG\t%INFO/CLNREVSTAT\t%INFO/CLNDN\t%SYMBOL\t%Gene\t%Feature\t%MANE_SELECT\t%Consequence\t%IMPACT\t%HGVSc\t%HGVSp\t%REVEL_score\t%CADD_phred\t%AlphaMissense_score\t%SpliceAI_pred_DS_AG\t%SpliceAI_pred_DS_AL\t%SpliceAI_pred_DS_DG\t%SpliceAI_pred_DS_DL[\t%GT\t%DP\t%GQ]\n' \
  -d -A tab \
  | gzip > variants.tsv.gz

-d expands multiple consequences onto separate lines. -A tab sets the separator for multi-value fields. Then:

-- duckdb
CREATE TABLE v AS
SELECT * FROM read_csv('variants.tsv.gz', delim='\t', header=false, nullstr='.',
  columns={'chrom':'VARCHAR','pos':'BIGINT','ref':'VARCHAR','alt':'VARCHAR','filter':'VARCHAR',
           'gnomad_af':'DOUBLE','nhomalt':'BIGINT','clnsig':'VARCHAR','clnrevstat':'VARCHAR','clndn':'VARCHAR',
           'symbol':'VARCHAR','gene':'VARCHAR','feature':'VARCHAR','mane':'VARCHAR','csq':'VARCHAR','impact':'VARCHAR',
           'hgvsc':'VARCHAR','hgvsp':'VARCHAR','revel':'DOUBLE','cadd':'DOUBLE','am':'DOUBLE',
           'ds_ag':'DOUBLE','ds_al':'DOUBLE','ds_dg':'DOUBLE','ds_dl':'DOUBLE',
           'gt':'VARCHAR','dp':'INTEGER','gq':'INTEGER'});
COPY v TO 'variants.parquet' (FORMAT PARQUET, COMPRESSION ZSTD);

Now the whole genome is one file you can slice interactively. This is the point of the exercise: integrated annotation plus a query layer is what turns a variant call set into something you can ask questions of, rather than a file you grep 2.

7. Filter down to a readable set

Four or five million rows is not a report. A standard first pass:

SELECT chrom, pos, ref, alt, symbol, csq, hgvsp, gnomad_af, clnsig, revel, cadd, gt, dp, gq
FROM v
WHERE filter = 'PASS'
  AND mane IS NOT NULL
  AND (gnomad_af IS NULL OR gnomad_af < 0.001)
  AND (impact = 'HIGH'
       OR (impact = 'MODERATE' AND revel > 0.7)
       OR greatest(ds_ag, ds_al, ds_dg, ds_dl) > 0.5)
ORDER BY gnomad_af NULLS FIRST;

That typically lands in the low hundreds of rows for one genome. A second, separate query on ClinVar alone is worth running, because a common variant can still carry a pathogenic assertion:

SELECT * FROM v
WHERE clnsig LIKE '%Pathogenic%' AND clnsig NOT LIKE '%Conflicting%'
  AND clnrevstat IN ('criteria_provided,_multiple_submitters,_no_conflicts',
                     'reviewed_by_expert_panel','practice_guideline')
  AND gt != '0/0';

Filter on CLNREVSTAT. ClinVar contains single-submitter assertions with no assertion criteria, and treating those as equivalent to expert-panel review is the single most common way to scare yourself over nothing.

Then check read evidence for anything that survives. Annotation pipelines that carry genotype quality and depth alongside the functional call let you discard artifacts before they reach a human reader 3. A “pathogenic” call at DP 6 with GQ 20 in a homopolymer is a sequencing error until proven otherwise. Load the CRAM in IGV and look at the actual reads.

8. Pharmacogenomics is a separate pipeline

Star-allele calling for genes like CYP2D6, CYP2C19, and DPYD depends on haplotype combinations and structural variation, not on single-variant consequence annotation. VEP will tell you a CYP2C19 variant is synonymous and be useless to you. Run PharmCAT, which takes a preprocessed VCF and produces named diplotypes with annotations drawn from curated pharmacogenomic guidelines 4. PharmCAT ships its own preprocessor script that normalizes and restricts the VCF to its position list. Run it on the file from step 3, not on the VEP output.

Pharmacogenomic results describe metabolizer phenotype. Any medication decision based on them belongs to a prescribing clinician.

9. Structural variants need their own path

Short-read SV callers (Manta, Delly, GRIDSS) emit VCFs where ALT is a symbolic allele like <DEL> and the interval is defined by POS and INFO/END. None of the above applies: bcftools norm has nothing to left-align, and exact-match database joins fail because two callers rarely agree on breakpoints to the base. Use AnnotSV or SvAnna for gene overlap and dosage sensitivity, and match against gnomAD-SV with reciprocal-overlap criteria (50% is the common threshold) rather than exact position. Visual inspection is worth more here than for SNVs, and there is a solid set of purpose-built viewers for structural variants from whole-genome data 5. Genome-wide overviews such as Circos-style plots built directly from VCF are useful for spotting large-scale patterns before you start filtering 6. If you have long reads, the SV picture is substantially better resolved, and the tooling and reference resources differ enough that you should follow a long-read-specific workflow rather than adapting this one 7.

Common problems

REF mismatch warnings from bcftools norm. Your FASTA is not the reference the caller used. Common causes: analysis set versus primary assembly (differs at the PAR and decoy contigs), or a hg19/GRCh37 mix. Compare md5 of a chromosome sequence, or just re-download the exact FASTA named in the VCF header.

gnomAD annotation rate far below 99%. Multiallelic records were not split, or your file is not left-aligned. Re-run step 3 and check that bcftools norm reports realigned records.

VEP output missing HGVS notation. --hgvsc/--hgvsp need --fasta to be supplied when running offline. It fails quietly.

CSQ field parsed into the wrong columns. Something between VEP and your parser changed field order, or a record has an empty trailing field. Always use bcftools +split-vep, which reads the field order from the header rather than assuming it.

Two INFO keys named AF. Your caller’s alternate allele fraction and gnomAD’s population frequency will collide silently. Rename on annotation with := and verify with bcftools view -h | grep 'ID=AF'.

Chromosome ordering errors on merge. If you split by chromosome and concatenate, bcftools concat requires consistent contig order across inputs and matching headers. Use bcftools concat -a with indexed files when in doubt.

Runs that are slow for no obvious reason. VEP with --everything and no --fork on a whole genome can take most of a day. Set --fork to your physical core count and pick your fields explicitly. If you are annotating many genomes, the whole flow parallelizes cleanly by chromosome shard and is a natural fit for a cloud batch pipeline rather than one long-running machine 8.

A HIGH-impact variant in a gene you recognize. Check the read pileup, check CLNREVSTAT, check whether the transcript is MANE Select, and check whether the gene is a known false-positive-prone region. Even in well-run diagnostic programs, the path from a candidate variant to a confirmed finding runs through clinical review and segregation evidence, not through a filter query 9. Take the specific variant, the HGVS notation, and the transcript ID to a clinical geneticist.

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. Peter Robinson, Manuel Holtgrewe. Variant Annotation. Computational Exome and Genome Analysis, 2017. https://doi.org/10.1201/9781315154770-15 ↩

  2. F. Anthony San Lucas, Gao Wang, Paul Scheet, et al. Integrated annotation and analysis of genetic variants from next-generation sequencing studies with variant tools. Bioinformatics, 2011. https://doi.org/10.1093/bioinformatics/btr667 ↩

  3. Zarko Manojlovic. Germline VCF Annotator: a lightweight pipeline for processing germline VCFs with robust variant extraction and read evidence quality control. 2026. https://doi.org/10.64898/2026.04.06.716730 ↩

  4. Teri E. Klein, Marylyn D. Ritchie. PharmCAT: A Pharmacogenomics Clinical Annotation Tool. Clinical Pharmacology & Therapeutics, 2017. https://doi.org/10.1002/cpt.928 ↩

  5. Toshiyuki T. Yokoyama, Masahiro Kasahara. Visualization tools for human structural variations identified by whole-genome sequencing. Journal of Human Genetics, 2019. https://doi.org/10.1038/s10038-019-0687-0 ↩

  6. E Drori, D Levy, P Smirin-Yosef, et al. CircosVCF: circos visualization of whole-genome sequence variations stored in VCF files. Bioinformatics, 2017. https://doi.org/10.1093/bioinformatics/btw834 ↩

  7. Unknown. Bioinformatics for human long-read whole genome sequencing. Nature Reviews Methods Primers, 2026. https://doi.org/10.1038/s43586-026-00528-w ↩

  8. Jeffrey G Reid, Andrew Carroll, Narayanan Veeraraghavan, et al. Launching genomics into the cloud: deployment of Mercury, a next generation sequence analysis pipeline. BMC Bioinformatics, 2014. https://doi.org/10.1186/1471-2105-15-30 ↩

  9. Dong Sun, Robert H Henderson, Emma Clement, et al. Analysis of familial exudative vitreoretinopathy (FEVR) cases in the UK 100 000 genomes project increases diagnostic rate and implicates heterozygous CTNND1 mutations in FEVR. Journal of Medical Genetics, 2025. https://doi.org/10.1136/jmg-2025-111083 ↩