Applied - Annotating a Variant Table

In one line: the full pipeline — VCF in, HGVS projection, population frequency, gene annotation, prioritised table out.

Uses: Applied - Reading VCF with pandas · AssemblyMapper · Merging and Joining · hgvs Normalizer · groupby · Applied - GFF and Genomic Intervals

This note ties together most of the vault. It is the end-to-end shape of a clinical or research variant analysis.

The pipeline

VCF  →  pandas frame  →  HGVS g.  →  normalize  →  pick transcript
     →  project to c./p.  →  join gnomAD  →  join gene annotation
     →  prioritise  →  report

Step 1: VCF to frame

import pandas as pd, numpy as np
vcf = read_vcf("cohort.vcf.gz")            # see [[Applied - Reading VCF with pandas]]

assert not vcf["ALT"].str.contains(",").any(), \
    "multi-allelic sites present — run: bcftools norm -m -any -f ref.fa"

Split and left-align multi-allelics upstream with bcftools. Doing it in pandas is possible and error-prone; bcftools norm -m -any -f reference.fa is one command and correct.

Step 2: build genomic HGVS

from bioutils.assemblies import make_ac_name_map
AC = {v: k for k, v in make_ac_name_map("GRCh38").items()}

def to_hgvs_g(chrom, pos, ref, alt):
    ac = AC[str(chrom).removeprefix("chr")]
    if len(ref) == 1 and len(alt) == 1:
        return f"{ac}:g.{pos}{ref}>{alt}"
    # indels: VCF anchors on the preceding base; strip it
    if len(ref) > len(alt) and ref.startswith(alt):          # deletion
        d = ref[len(alt):]
        s = pos + len(alt)
        return f"{ac}:g.{s}_{s + len(d) - 1}del" if len(d) > 1 else f"{ac}:g.{s}del"
    if len(alt) > len(ref) and alt.startswith(ref):          # insertion
        i = alt[len(ref):]
        return f"{ac}:g.{pos + len(ref) - 1}_{pos + len(ref)}ins{i}"
    return f"{ac}:g.{pos}_{pos + len(ref) - 1}delins{alt}"   # complex

vcf["hgvs_g"] = [to_hgvs_g(*r) for r in
                 zip(vcf["CHROM"], vcf["POS"], vcf["REF"], vcf["ALT"])]

The NC_ accession, not chr17 — the accession version encodes the assembly, which is exactly what chr17 fails to do.

Normalize immediately after building, because VCF is left-aligned and HGVS is 3'-most:

vcf["hgvs_g"] = vcf["hgvs_g"].map(
    lambda s: str(hn.normalize(hp.parse_hgvs_variant(s)))
)

Skipping this step is the number one cause of failed ClinVar lookups. See hgvs Pitfalls.

Step 3: transcript selection

MANE = pd.read_csv("MANE.GRCh38.v1.4.summary.txt.gz", sep="\t")
MANE_BY_GENE = dict(zip(MANE["symbol"], MANE["RefSeq_nuc"]))
MANE_SET = set(MANE["RefSeq_nuc"])

CLINICAL_OVERRIDE = {
    "BRCA1": "NM_007294.4",
    "CFTR":  "NM_000492.4",
}

def pick_transcript(gene, available):
    if gene in CLINICAL_OVERRIDE and CLINICAL_OVERRIDE[gene] in available:
        return CLINICAL_OVERRIDE[gene]
    mane = MANE_BY_GENE.get(gene)
    if mane in available:
        return mane
    for tx in sorted(available):
        if tx in MANE_SET:
            return tx
    return sorted(available)[0] if available else None

A genomic variant has one c. description per transcript, and a gene can have twenty. Choosing is policy, not computation. MANE Select is the right default; gene-specific clinical conventions override it, because a variant reported on a different transcript than the literature uses will not match anything.

Step 4: project

from hgvs.exceptions import HGVSError

def annotate_one(hgvs_g, gene=None):
    out = {"hgvs_c": None, "hgvs_p": None, "transcript": None,
           "consequence_region": None, "error": None}
    try:
        var_g = hp.parse_hgvs_variant(hgvs_g)
        txs = am.relevant_transcripts(var_g)
        if not txs:
            out["consequence_region"] = "intergenic"; return out

        tx = pick_transcript(gene, txs)
        out["transcript"] = tx
        var_c = am.g_to_c(var_g, tx)
        out["hgvs_c"] = str(var_c)

        off = getattr(var_c.posedit.pos.start, "offset", 0)
        base = var_c.posedit.pos.start.base
        out["consequence_region"] = (
            "intronic" if off != 0 else
            "5_utr" if base < 0 else
            "3_utr" if str(var_c.posedit.pos.start).startswith("*") else
            "coding"
        )

        try:
            out["hgvs_p"] = str(am.c_to_p(var_c))
        except HGVSError:
            out["hgvs_p"] = "p.?"
    except HGVSError as e:
        out["error"] = f"{type(e).__name__}: {e}"
    return out

ann = pd.DataFrame([annotate_one(g, gene) for g, gene
                    in zip(vcf["hgvs_g"], vcf.get("gene", [None]*len(vcf)))])
vcf = vcf.join(ann)
print(vcf["error"].notna().value_counts())

Note the region classification reads offset and datum, not just base — c.87+3 is intronic and c.87 is not. See SequenceVariant.

Step 5: join population frequencies

gnomad = pd.read_parquet("gnomad_v4_af.parquet")

before = len(vcf)
vcf = vcf.merge(
    gnomad[["chrom", "pos", "ref", "alt", "af", "af_popmax", "nhomalt"]],
    left_on=["CHROM", "POS", "REF", "ALT"],
    right_on=["chrom", "pos", "ref", "alt"],
    how="left", validate="many_to_one",
)
assert len(vcf) == before, "merge changed the row count"
print(f"{vcf['af'].isna().mean():.1%} not in gnomAD")

validate="many_to_one" and the row-count assertion. A non-unique join key silently multiplies rows, and every downstream count is then wrong. See Merging and Joining.

Both sides must be normalized identically, or the join misses. gnomAD is left-aligned VCF-style; if you normalized to HGVS 3'-most, join on the original VCF coordinates (as above) rather than the HGVS string.

Step 6: join gene annotation

genes = pd.read_parquet("gene_annotation.parquet")
genes["gene_id"] = genes["gene_id"].str.split(".").str[0]      # strip versions

vcf = vcf.merge(genes[["gene_symbol", "gene_id", "biotype",
                       "pli", "oe_lof_upper", "disease", "inheritance"]],
                left_on="gene", right_on="gene_symbol",
                how="left", validate="many_to_one")

pLI and LOEUF (oe_lof_upper) are gnomAD constraint metrics — how depleted a gene is of loss-of-function variants in the general population. LOEUF < 0.35 marks a highly constrained gene, which raises the prior that a LoF variant there is consequential.

Step 7: prioritise

IMPACT_ORDER = ["MODIFIER", "LOW", "MODERATE", "HIGH"]
vcf["impact"] = pd.Categorical(vcf["impact"], categories=IMPACT_ORDER, ordered=True)

vcf["is_rare"] = vcf["af_popmax"].fillna(0) < 0.001        # ← fillna!
vcf["is_pass"] = vcf["FILTER"] == "PASS"
vcf["is_coding"] = vcf["consequence_region"] == "coding"
vcf["is_constrained"] = vcf["oe_lof_upper"].fillna(1.0) < 0.35

vcf["priority"] = (
    vcf["is_pass"].astype(int) * 2
    + vcf["is_rare"].astype(int) * 3
    + (vcf["impact"] >= "MODERATE").astype(int) * 3
    + vcf["is_constrained"].astype(int) * 2
    + vcf["gene"].isin(PANEL).astype(int) * 4
)

shortlist = (vcf.query("is_pass and is_rare and impact >= 'MODERATE'")
                .sort_values("priority", ascending=False))

fillna(0) on the allele frequency is critical. A variant absent from gnomAD has NA frequency; NA fails < 0.001, so without the fill you silently exclude the rarest variants — the exact opposite of the intent. See Missing Data in pandas.

The ordered categorical makes impact >= "MODERATE" work. Without ordered=True it raises; without the categorical it compares alphabetically and "MODIFIER" > "MODERATE" is True. See Categorical dtype.

This score is a triage heuristic, not a classification. ACMG/AMP criteria are what determine pathogenicity, and they require evidence this table does not contain.

Step 8: filter accounting

def report(df, mask, label):
    print(f"{label:32s} {mask.sum():>8,} / {len(df):>8,} ({mask.mean():6.1%})")
    return df.loc[mask]

print(f"{'input':32s} {len(vcf):>8,}")
d = report(vcf, vcf["is_pass"], "PASS")
d = report(d, d["is_rare"], "rare (popmax < 0.1%)")
d = report(d, d["impact"] >= "MODERATE", "moderate/high impact")
d = report(d, d["gene"].isin(PANEL), "on panel")

Filter accounting belongs in every pipeline and every methods section. It turns "we found nothing" into "step 3 removed 99.8% of variants, which is unexpected".

Step 9: report

COLS = ["gene", "hgvs_c", "hgvs_p", "hgvs_g", "transcript",
        "consequence_region", "impact", "af_popmax", "nhomalt",
        "oe_lof_upper", "QUAL", "DP", "priority", "disease", "inheritance"]

shortlist[COLS].to_csv("prioritised_variants.tsv", sep="\t", index=False)

with pd.ExcelWriter("report.xlsx") as w:
    shortlist[COLS].to_excel(w, sheet_name="shortlist", index=False)
    vcf[COLS].to_excel(w, sheet_name="all_variants", index=False)
    pd.DataFrame([PROVENANCE]).to_excel(w, sheet_name="provenance", index=False)

The provenance sheet is not optional. hgvs version, UTA snapshot, SeqRepo snapshot, assembly, alignment method, MANE version, gnomAD version. Without them the table cannot be regenerated or audited. See Applied - Reproducible Environment.

Visual summary

fig, axes = plt.subplots(2, 2, figsize=(11, 8), layout="constrained")

sns.countplot(data=vcf, x="consequence_region", ax=axes[0,0],
              order=vcf["consequence_region"].value_counts().index)
sns.histplot(data=vcf, x="af_popmax", bins=50,
             log_scale=(True, False), ax=axes[0,1])
sns.countplot(data=vcf, x="impact", ax=axes[1,0], order=IMPACT_ORDER)
sns.scatterplot(data=vcf, x="oe_lof_upper", y="priority",
                hue="impact", alpha=0.5, ax=axes[1,1])

Caveats worth stating out loud

One transcript's consequence is not the whole story. A variant benign on MANE Select may be damaging on another expressed transcript.
p.(...) is predicted, not observed.
Absent from gnomAD ≠ rare. It may be in a poorly-covered region. Check gnomAD's coverage track.
Filter status reflects the caller's model, not truth.
Constraint metrics are gene-level, not variant-level.
This is triage, not diagnosis.

Common mistakes

Not normalizing before joining to gnomAD or querying ClinVar.
Multi-allelic sites left unsplit.
NA frequencies failing the rarity filter.
Merges without validate=.
Not recording the transcript per row.
Version suffixes breaking gene-ID joins.
Alphabetical impact ordering.
Mixing GRCh37 and GRCh38.
No provenance in the output.
Presenting a priority score as a classification.

Applied - Annotating a Variant Table

Applied - Annotating a Variant Table

The pipeline

Step 1: VCF to frame

Step 2: build genomic HGVS

Step 3: transcript selection

Step 4: project

Step 5: join population frequencies

Step 6: join gene annotation

Step 7: prioritise

Step 8: filter accounting

Step 9: report

Visual summary

Caveats worth stating out loud

Common mistakes

See also