Applied - Annotating a Variant Table

Core Principle: End-to-end variant annotation pipelines transform raw VCF records into prioritized clinical tables by integrating HGVS nomenclature, gene models, and population allele frequencies.

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.

See also

Applied - Reading VCF with pandas · Applied - Variant Normalization Pipeline · AssemblyMapper · Projection c to g to p · Merging and Joining · hgvs Pitfalls · Applied - Reproducible Environment · Applied - GFF and Genomic Intervals

Exercise

Build the whole annotation pipeline in one function. hgvs is not installed here, so the projection step is replaced by a join on the plain chrom:pos:ref:alt key; every other stage of the lesson -- quality filter, population-frequency join, panel join, prioritisation, filter accounting -- is present.

Write annotate(text), returning a list of output lines. The input arrives in three sections separated by lines containing exactly --: the variant calls (chrom, pos, ref, alt, gene, qual, tab separated), a population frequency table (chrom:pos:ref:alt and its allele frequency), and a gene panel (gene and its role). A section may be empty. Ignore blank lines and any variant row that does not have exactly six fields. Drop calls with qual below 30, then join the frequencies and drop anything whose frequency exceeds 0.01; a key absent from the frequency table is novel, not common, and must survive. Join the panel on gene name; a gene missing from the panel has role none.

Assign each surviving call a tier: 1 for a panel gene with no frequency entry, 2 for a panel gene with frequency at or below 0.001, 3 for any other panel gene, and 4 for a gene that is not on the panel. Print input=<n>, then low_qual=<n>, then common=<n>, then kept=<n>. Follow those four lines with one line per surviving call, ordered by tier ascending, then quality descending, then key ascending, holding the tier, the key, the gene, the role and the frequency separated by tabs; format the frequency with exactly four decimal places, or print novel when there was no entry. If nothing survives, print no variants pass after the accounting lines instead.