Bio Applied Variant Annotation

Applied Variant Annotation: VEP, snpEff, gnomAD, ClinVar, dbNSFP

When to Use

• You have a called, normalized VCF and need consequence terms (missense, frameshift, splice), HGVS notation, and IMPACT (HIGH/MODERATE/LOW/MODIFIER) per variant.
• You need to attach population allele frequency (gnomAD) to distinguish rare candidates from common polymorphisms.
• You need ClinVar clinical significance and dbNSFP in-silico predictor scores (CADD, REVEL, SIFT, PolyPhen) joined onto each variant record.
• You are comparing VEP vs. snpEff output, or need to parse the CSQ/ANN INFO string programmatically downstream.
• Explicitly NOT for variant calling (see bio-applied-variant-calling-and-snp-analysis) or ACMG clinical classification (see bio-applied-clinical-genomics) — this skill only produces the annotated VCF those steps consume.

Version Compatibility

• Ensembl VEP 114-116 (2026 releases; cache/FASTA must match the VEP release and genome build, e.g. GRCh38 release 116).
• snpEff/SnpSift 5.2 (build GRCh38.p14 or matching database).
• dbNSFP 4.9a academic (hg38); gnomAD v4.1 (joint exomes+genomes, GRCh38).
• Python: cyvcf2 >=0.31 for parsing annotated VCFs (pip install cyvcf2).

Prerequisites

• Install: conda install -c bioconda ensembl-vep=116 snpeff=5.2 cyvcf2 (or mamba); VEP also needs --cache --dir_cache <path> downloaded via vep_install, and a genome FASTA for --hgvs.
• Input VCF must already be left-normalized and split (multiallelics decomposed) — see bio-variant-calling-variant-normalization — annotators evaluate one ALT per record correctly only after normalization.
• Concepts: VCF INFO fields and headers (bio-variant-calling-vcf-basics), variant calling upstream (bio-applied-variant-calling-and-snp-analysis).

Core Workflows

Goal: Annotate a VCF with consequence/HGVS/impact plus gnomAD AF, ClinVar significance, and dbNSFP predictor scores in one VEP run.

Approach: Build the vep CLI command as an argv list (never a shell string, to avoid quoting bugs with plugin args) combining cache-based consequence prediction, the --custom ClinVar VCF, and the dbNSFP plugin; run it with subprocess.run.

import subprocess


def run_vep(
    input_vcf: str, output_vcf: str, cache_dir: str, fasta: str,
    clinvar_vcf: str, dbnsfp_gz: str, fork: int = 4,
) -> subprocess.CompletedProcess:
    """Annotate a VCF with Ensembl VEP: consequence, HGVS, impact,
    gnomAD exome/genome AF, ClinVar significance, and dbNSFP scores.

    Args:
        input_vcf: normalized, decomposed input VCF (one ALT per record).
        output_vcf: path for the VEP-annotated VCF (--vcf output format).
        cache_dir: local VEP cache dir (from `vep_install --CACHEDIR`).
        fasta: reference FASTA matching the cache's genome build (for --hgvs).
        clinvar_vcf: bgzipped+tabixed ClinVar VCF (same build as cache).
        dbnsfp_gz: bgzipped, tabixed dbNSFP file for the dbNSFP plugin.
    Returns:
        The completed subprocess.CompletedProcess (raises on non-zero exit).
    """
    cmd = [
        "vep",
        "--input_file", input_vcf, "--output_file", output_vcf,
        "--vcf", "--force_overwrite",
        "--cache", "--offline", "--dir_cache", cache_dir, "--fasta", fasta,
        "--species", "homo_sapiens", "--assembly", "GRCh38",
        "--everything", "--hgvs", "--symbol", "--canonical", "--numbers", "--domains",
        "--af_gnomade", "--af_gnomadg",
        "--custom", f"file={clinvar_vcf},short_name=ClinVar,format=vcf,type=exact,"
                    "fields=CLNSIG%CLNREVSTAT%CLNDN",
        "--plugin", f"dbNSFP,{dbnsfp_gz},SIFT_score,Polyphen2_HDIV_score,CADD_phred,REVEL_score",
        "--fork", str(fork), "--no_stats",
    ]
    return subprocess.run(cmd, check=True, capture_output=True, text=True)

Goal: Programmatically parse the VEP CSQ (or snpEff ANN) INFO string per variant so downstream Python code can filter/report.

Approach: Read the ##INFO=<ID=CSQ,...Format: a|b|c...> header line to get field order (VEP's format is not fixed across runs — it depends on which flags were used), then zip that order against each |-delimited transcript block; a variant can carry multiple transcript annotations comma-separated.

from cyvcf2 import VCF


def parse_vep_csq(vcf_path: str, max_gnomad_af: float = 0.01) -> list[dict]:
    """Parse VEP CSQ annotations from an annotated VCF.

    Reads the CSQ Format string from the VCF header (field order depends on
    which VEP flags were used, so it must never be hardcoded), then expands
    each variant's (possibly multi-transcript) CSQ string into records.

    Args:
        vcf_path: path to a VEP-annotated VCF (bgzipped or plain).
        max_gnomad_af: keep only records at/below this gnomAD exome AF
                       (rare-variant triage; None disables filtering).
    Returns:
        List of dicts, one per transcript annotation, with keys from the
        CSQ header plus 'chrom', 'pos', 'ref', 'alt'.
    """
    vcf = VCF(vcf_path)
    csq_fields = None
    for h in vcf.header_iter():
        info = h.info(extra=True)
        if info.get("ID") == "CSQ":
            csq_fields = info["Description"].split("Format: ")[1].strip('"').split("|")
            break
    if csq_fields is None:
        raise ValueError("No CSQ INFO field found -- was this VCF annotated by VEP?")

    records = []
    for variant in vcf:
        csq_raw = variant.INFO.get("CSQ")
        if not csq_raw:
            continue
        for transcript_csq in csq_raw.split(","):
            values = transcript_csq.split("|")
            rec = dict(zip(csq_fields, values))
            af = rec.get("gnomADe_AF") or ""
            if max_gnomad_af is not None and af not in ("", "."):
                if float(af) > max_gnomad_af:
                    continue
            rec.update(chrom=variant.CHROM, pos=variant.POS, ref=variant.REF, alt=variant.ALT[0])
            records.append(rec)
    return records

Goal: Cross-check or substitute VEP with snpEff when a lighter, Java-only annotator is preferred (e.g. non-human genomes with a snpEff-built database but no VEP cache).

Approach: Run snpEff ann to produce ANN fields, then use SnpSift to join ClinVar and dbNSFP annotations onto the same VCF as separate filter/annotate passes, piping VCF text through each stage.

import subprocess


def annotate_with_snpeff(
    input_vcf: str, genome_db: str, clinvar_vcf: str, dbnsfp_txt_gz: str,
) -> str:
    """Run snpEff consequence annotation, then SnpSift annotate (ClinVar)
    and SnpSift dbnsfp (predictor scores), streaming VCF text through
    each step instead of writing intermediate files.

    Args:
        input_vcf: normalized, decomposed input VCF.
        genome_db: snpEff database name (e.g. 'GRCh38.p14').
        clinvar_vcf: bgzipped+tabixed ClinVar VCF with an ID column snpEff/SnpSift can match.
        dbnsfp_txt_gz: bgzipped, tabixed dbNSFP text file.
    Returns:
        Annotated VCF text (ANN field from snpEff plus joined INFO fields).
    """
    snpeff = subprocess.run(
        ["snpEff", "-v", genome_db, input_vcf], check=True, capture_output=True, text=True,
    )
    annotated = subprocess.run(
        ["SnpSift", "annotate", "-id", clinvar_vcf, "-"],
        input=snpeff.stdout, check=True, capture_output=True, text=True,
    )
    scored = subprocess.run(
        ["SnpSift", "dbnsfp", "-db", dbnsfp_txt_gz,
         "-f", "SIFT_score,Polyphen2_HDIV_score,CADD_phred,REVEL_score", "-"],
        input=annotated.stdout, check=True, capture_output=True, text=True,
    )
    return scored.stdout

Pitfalls

• Multiple transcript blocks per variant: CSQ/ANN pack one entry per overlapping transcript; picking "the" consequence means choosing a policy (canonical transcript, MANE Select, or worst-consequence-across-transcripts) — don't silently take the first field.
• Un-normalized input breaks left-alignment matching: annotate only after normalizing multiallelics (bcftools norm -m-) and indels, or VEP/snpEff and your gnomAD/ClinVar joins will disagree on which ALT they're describing.
• gnomAD exomes vs. genomes are different denominators: gnomADe_AF (exomes) and gnomADg_AF (genomes) can differ substantially at low-coverage or non-exonic sites — check both, and prefer gnomADg_AF for non-coding variants.
• CSQ field order is not fixed: it depends on exactly which --everything/plugin flags were used for that run; always parse it from the VCF header's Format: string, never hardcode column positions.
• VEP cache/FASTA build mismatch silently gives wrong HGVS/consequences: cache, FASTA, and --assembly must all be the same genome build (e.g. all GRCh38) or annotations will be plausible-looking but wrong.
• IMPACT is not pathogenicity: HIGH impact (e.g. stop_gained) is a consequence-severity heuristic, not a clinical call — pathogenicity classification is a separate downstream step (ACMG/AMP).

See Also

• bio-applied-variant-calling-and-snp-analysis — upstream variant calling that produces the VCF this skill annotates.
• bio-applied-clinical-genomics — downstream ACMG/AMP pathogenicity classification using ClinVar/gnomAD evidence gathered here.
• bio-applied-functional-annotation — genome/gene-level functional annotation (as opposed to per-variant consequence).
• bio-core-biological-databases — general patterns for querying Ensembl/NCBI databases used as annotation sources.