Bio Variant Annotation

Version Compatibility

Reference examples tested with: bcftools 1.19+, VEP 110+, SnpEff 5.2+, ANNOVAR 2020Jun07+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags
• Note: SnpEff and SnpSift use single-dash -version, not --version

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Note: gnomAD v2 is GRCh37; v3/v4 are GRCh38. MANE transcripts exist only on GRCh38. Confirm the build of every annotation source matches the VCF before annotating.

Variant Annotation

"Annotate my variants with functional and clinical information" -> Map each variant onto a transcript model, classify its coding/splice consequence, and attach population frequency and pathogenicity evidence.

• CLI: vep (Ensembl), snpEff/SnpSift, table_annovar.pl (ANNOVAR), bcftools csq/annotate
• Python: cyvcf2 to parse VEP CSQ / SnpEff ANN strings; bcftools +split-vep to flatten them

The governing principle: annotation is not deterministic

A variant's consequence is not a property of the variant. It is a property of the tuple (variant, transcript model, engine, engine version, parameter set). Change any element and the reported consequence, the HGVS string, and downstream the ACMG PVS1 eligibility can all change. The single most damaging naive belief in clinical genomics is that a VCF line has one true annotation. It does not. On any discordant result the first question is never "what does the variant do" but "which transcript and which tool+version produced that call."

The decision is therefore NOT to find the "right" tool. It is to PIN every axis and record it on the report: genome build, transcript set (prefer MANE Select on GRCh38), engine + version, predictor + version, gnomAD version. Reproducibility comes from pinning, not from picking. Never compare a variant annotated on RefSeq to one annotated on Ensembl.

Three independent axes of non-determinism: the transcript SET (RefSeq vs Ensembl/GENCODE vs MANE), the transcript-SELECTION heuristic (canonical, worst-consequence, --pick, MANE Select), and the ENGINE itself (VEP/SnpEff/ANNOVAR encode different splice-region widths, up/downstream windows, HGVS-shifting rules, and consequence-severity orderings). Discordance concentrates in indels, splice-region, and multi-transcript genes; loss-of-function calls that drive PVS1 are among the least concordant across engines (McLaren 2016 Genome Biol 17:122; Cingolani 2012 Fly 6(2):80-92; Wang 2010 Nucleic Acids Res 38(16):e164).

Normalize before annotation, and never hand-derive HGVS from POS

Normalization is mandatory: the same variant represented differently produces different annotations.

# -m-any splits multiallelic records to biallelic so each ALT gets its own annotation
bcftools norm -f reference.fa -m-any input.vcf.gz -Oz -o normalized.vcf.gz

The HGVS 3'-rule vs VCF left-align clash is a genuine, still-live trap. VCF normalization requires indels left-aligned (most 5' on the forward genomic strand; Tan 2015 Bioinformatics 31(13):2202-2204). HGVS mandates the opposite: the 3'-rule places an indel in a repeat at the most 3' position with respect to the transcript. For a plus-strand gene, transcript-3' is the rightmost genomic position, the opposite end from VCF left-alignment; for a minus-strand gene the two can coincide by accident of strand. Net effect: a correctly left-aligned VCF POS and a correct HGVS c. string for the same indel can point to different repeat units.

Rules:

• Left-align + normalize the VCF BEFORE annotation (idempotent representation for matching/merging).
• Rely on the engine's HGVS generator to apply the transcript 3'-shift (VEP does; verify SnpEff/ANNOVAR per version). Never hand-derive an HGVS string from POS.
• When matching a patient variant to a ClinVar/literature HGVS, normalize BOTH to the same representation first. dup vs ins describe the same event but do not string-match. See variant-calling/variant-normalization.
• Adjacent SNVs in one codon (an MNV split by the caller) can be individually benign but jointly change the amino acid; per-SNV annotation silently misannotates. Use phase-aware annotation.

Choosing the transcript set

Decision: report on MANE Select (plus MANE Plus Clinical where assigned), not worst-consequence. Worst-consequence reports a canonical splice change in a minor non-expressed isoform as "splice" even when MANE Select is intronic, manufacturing false PVS1 candidates; restricting to MANE Select alone can miss a variant that only hits a MANE Plus Clinical isoform -- which is exactly why that tier exists. The Ensembl "canonical" transcript is frequently NOT the MANE Select one, so migrating a pipeline to MANE changes some reported c. coordinates (expected, not erroneous). MANE is GRCh38-only; GRCh37 pipelines must lift over or maintain their own per-gene transcript pins.

Why --pick is dangerous for clinical use. VEP by default reports all consequences for all overlapping transcripts. --pick collapses to one block per variant using an ordered heuristic whose defaults are canonical status, biotype, consequence rank, then transcript length, then finally accession order. The late tiebreakers are not clinically motivated: when transcripts tie, --pick can let transcript length or alphanumeric ENST order decide which consequence a real patient gets, can pick per-variant (so variant A and variant B in one gene land on different transcripts, destroying coordinate consistency), and can silently hide a PVS1-eligible consequence behind a benign one. Defensible configurations: pin MANE Select (with Plus Clinical), or pin the lab's validated per-gene list. If --pick-style collapse is used at all, constrain it so MANE leads and length/accession never decide:

# --pick_order forces MANE first; length/accession can no longer choose the reported transcript
vep -i norm.vcf --vcf --cache --offline --assembly GRCh38 \
    --mane_select --pick --pick_order mane_select,canonical,biotype,rank -o out.vcf

Consequence, impact, and NMD (the PVS1 hinge)

Anchor the vocabulary to Sequence Ontology (SO) terms (missense_variant, stop_gained, frameshift_variant, splice_donor_variant, splice_acceptor_variant, start_lost, stop_lost, inframe_deletion, synonymous_variant, ...), which VEP emits natively. SnpEff/ANNOVAR map to mostly-equivalent terms but differ on splice-region width and finer intronic terms (VEP adds splice_donor_5th_base_variant, splice_polypyrimidine_tract_variant), so term-level matching across engines is impossible.

Impact buckets are NOT evidence. SnpEff HIGH/MODERATE/LOW/MODIFIER and ANNOVAR exonic;splicing groupings are triage conveniences. HIGH lumps stop_gained, frameshift, and canonical splice together, but whether any of these earns PVS1 depends on NMD, exon location, and the gene's LOF mechanism -- none of which the bucket knows. Treating "HIGH impact" as "PVS1 met" is a classic error.

The NMD 50-nt / last-exon rule governs PVS1 strength (Abou Tayoun 2018 Hum Mutat 39(11):1517-1524). A premature termination codon (PTC) more than ~50-55 nt upstream of the last exon-exon junction triggers nonsense-mediated decay -> no protein -> strong LOF. A PTC in the last exon, or within ~50 nt of the final junction (3' end of the penultimate exon), escapes NMD -- the truncated protein IS made, so full-strength PVS1 on the "NMD -> no protein" logic is unjustified (downgrade to PVS1_Strong/Moderate/Supporting depending on how much protein / which domains are lost). Single-exon genes have no junctions, so a nonsense variant escapes NMD by construction; PTCs near the start can reinitiate downstream. PVS1 also requires that LOF is the established disease mechanism -- a null allele in a gain-of-function/dominant-negative gene must NOT fire PVS1. The ACMG combining lives in variant-calling/clinical-interpretation; this skill supplies the consequence + NMD inputs it needs.

Annotation engines: run commands

Fix the transcript set and tool+version in the pipeline and record them; that, not tool choice, is what makes annotation reproducible.

Ensembl VEP

Goal: Annotate consequence, HGVS, impact, frequencies, and plugin predictions against Ensembl/MANE.

Approach: Run offline against the cache with explicit assembly and transcript selection; add predictor plugins rather than relying on the built-in SIFT/PolyPhen.

# --everything enables --hgvs --symbol --canonical --af --af_gnomade --af_gnomadg --sift b
# --polyphen b --pubmed etc. Prefer --mane_select over the enabled --canonical for reporting.
vep -i norm.vcf.gz -o out.vcf --vcf --cache --offline \
    --species homo_sapiens --assembly GRCh38 --everything --mane_select --fork 4

# Calibrated predictors as plugins (one calibrated tool per evidence type -- see below)
vep -i norm.vcf.gz -o out.vcf --vcf --cache --offline \
    --plugin dbNSFP,dbNSFP4.3a.gz,REVEL_score,CADD_phred \
    --plugin AlphaMissense,file=AlphaMissense_hg38.tsv.gz \
    --plugin SpliceAI,snv=spliceai_snv.vcf.gz,indel=spliceai_indel.vcf.gz

SnpEff / SnpSift

Goal: Fast batch effect annotation plus database cross-referencing.

Approach: snpEff ann against a prebuilt genome database, then chain SnpSift annotate/filter.

# Human GRCh38 DB expands to 3-4 GB in memory; give the JVM >= 8 GB or it OOMs/thrashes
snpEff -Xmx8g ann GRCh38.105 norm.vcf > out.vcf
snpEff -Xmx8g ann GRCh38.105 norm.vcf | SnpSift annotate clinvar.vcf.gz > annotated.vcf

ANNOVAR

Goal: Table-driven gene/frequency/pathogenicity annotation.

Approach: table_annovar.pl with paired -protocol/-operation lists (g=gene, f=filter, r=region).

table_annovar.pl norm.vcf humandb/ -buildver hg38 -out annotated -remove \
    -protocol refGene,gnomad30_genome,clinvar_20230416,dbnsfp42a \
    -operation g,f,f,f -nastring . -vcfinput

bcftools csq / annotate

Goal: Lightweight consequence prediction (csq) and database field transfer (annotate) without a full engine.

Approach: csq maps variants to a GFF3 and emits a BCSQ field; annotate -c copies ID/INFO columns from a position-matched source.

bcftools csq -f reference.fa -g genes.gff3.gz norm.vcf.gz -Oz -o csq.vcf.gz   # adds BCSQ
bcftools annotate -a dbsnp.vcf.gz -c ID norm.vcf.gz -Oz -o rsid.vcf.gz        # copy rsIDs

See usage-guide.md for BED/TAB annotation, field removal, --set-id, chromosome renaming, and database download recipes.

Pathogenicity predictors: one calibrated tool, not a stack

Two structural problems pervade this literature. (1) Circularity: most predictors train on ClinVar/HGMD labels, so benchmarking or ACMG-calibrating on those same databases is partly self-referential; a headline "AUC 0.9x" is optimistic on truly novel variants. (2) Ensembles ingest each other: REVEL is a random forest over 13 component scores including SIFT and PolyPhen, so "REVEL agrees with PolyPhen" is not independent corroboration -- PolyPhen is inside REVEL. Independence between evidence lines is the load-bearing assumption of the ACMG points system; stacking correlated predictors silently over-calls pathogenic.

Therefore: use exactly ONE calibrated predictor per evidence type (missense; splicing), applied at its calibrated strength.

SIFT/PolyPhen alone are near-worthless now. In the ClinGen SVI calibration (Pejaver 2022 Am J Hum Genet 109(12):2163-2177) neither reached even Supporting strength for PP3; both call a large fraction of all missense "damaging" (low positive predictive value on rare variants); and both are components of REVEL, so quoting them alongside it double-counts. Legacy pipelines surfacing "SIFT: deleterious, PolyPhen: probably damaging" prominently are decorative, not evidentiary.

CADD >= 20 is not "pathogenic." In Pejaver 2022 raw CADD did not reach Supporting for PP3, and the developer-recommended CADD >= 20 mapped to Moderate evidence for benign -- an inversion of how CADD 20 is casually used. Reserve CADD for its intended non-coding/genome-wide ranking.

Calibrated REVEL thresholds (Pejaver 2022). PP3_Supporting >= 0.644 and BP4_Supporting <= 0.290 are the well-reproduced values. The Moderate/Strong REVEL cutoffs (commonly quoted as PP3_Moderate >= 0.773, PP3_Strong >= 0.932; BP4_Moderate <= 0.183, BP4_Strong <= 0.016) come from the supplementary tables and are not uniformly reproduced -- verify against the Pejaver 2022 supplement / current ClinGen SVI recommendation table before hard-coding, rather than treating them as fixed. PP3 and BP4 are mutually exclusive by construction; only tools reaching >= Strong in the calibration qualify.

SpliceAI delta scores. Per variant, SpliceAI emits four deltas (acceptor gain/loss, donor gain/loss), each 0-1; the max is the headline. Developer-recommended interpretation: 0.2 high recall, 0.5 recommended, 0.8 high precision. Caveats: (i) know whether the pipeline uses the masked or raw model; (ii) the default scoring window is +/-50 bp -- deep-intronic/pseudoexon effects need a widened window (up to +/-10 kb) or are missed; (iii) a delta is a prediction, and converting it to PS3/PP3 strength needs the ClinGen splicing calibration, not the raw cutoffs; (iv) it does not report the RESULT (exon skip vs intron retention), which is what determines PVS1. Pangolin (Zeng 2022 Genome Biol 23:103) is an emerging tissue-aware alternative -- check current ClinGen splicing guidance. AlphaMissense and the splicing calibrations are still-evolving; verify the current ClinGen SVI approved-tool list before standardizing on one.