Applied - Variant Normalization Pipeline

In one line: taking a column of messy, clinician-typed variant strings and turning it into validated, normalized, comparable variants — with a full accounting of what failed and why.

Uses: hgvs Parser · hgvs Validator · hgvs Normalizer · UTA and Data Providers · SeqRepo · DataFrame · hgvs Pitfalls

The problem

A spreadsheet from a clinical lab:

NM_000546.6:c.215C>G
NM_000546.6 c.215C>G          ← space instead of colon
c.215C>G                       ← no accession
NM_000546:c.215C>G             ← no version
NM_000546.6:c.215C&gt;G          ← HTML-escaped
NM_000546.6:c.215C–G           ← en-dash from a word processor
p.Pro72Arg                     ← protein only
TP53 P72R                      ← legacy shorthand
IVS3+2T>C                      ← legacy intron notation

Expect a 5–20% raw parse failure rate. The job is to fix what is safely fixable, report the rest, and never guess at semantics.

Setup

import os, re, hgvs
import hgvs.parser, hgvs.dataproviders.uta
import hgvs.normalizer, hgvs.validator, hgvs.assemblymapper
from hgvs.exceptions import (HGVSError, HGVSParseError, HGVSInvalidVariantError,
                             HGVSInvalidIntervalError, HGVSDataNotAvailableError)
import pandas as pd

hp  = hgvs.parser.Parser()
hdp = hgvs.dataproviders.uta.connect(pooling=True)
hn  = hgvs.normalizer.Normalizer(hdp)
hv  = hgvs.validator.Validator(hdp)
am  = hgvs.assemblymapper.AssemblyMapper(
    hdp, assembly_name="GRCh38", alt_aln_method="splign",
    replace_reference=True, normalize=True)

PROVENANCE = {
    "hgvs_version": hgvs.__version__,
    "uta_version": hdp.data_version(),
    "seqrepo_dir": os.environ.get("HGVS_SEQREPO_DIR", "remote"),
    "assembly": am.assembly_name,
    "alt_aln_method": "splign",
}
print(PROVENANCE)

Build every object once. Constructing AssemblyMapper per variant is a 100× slowdown. Record the provenance — the UTA and SeqRepo versions determine your output. See hgvs Pitfalls.

Step 1: safe cleanup

def tidy(s, default_ac=None):
    """Fix formatting only. NEVER change semantics."""
    if not isinstance(s, str):
        return None
    s = s.strip()
    s = s.replace("&gt;", ">").replace("&lt;", "<")
    s = s.replace("–", "-").replace("—", "-")   # en/em dash
    s = s.replace("→", ">")                           # → arrow
    s = re.sub(r"\s+", "", s)
    s = re.sub(r"^([A-Z_]+\d+\.\d+)(c\.|g\.|n\.|p\.)", r"\1:\2", s)  # missing colon
    if ":" not in s and default_ac:
        s = f"{default_ac}:{s}"
    return s

Formatting only. Fixing an en-dash is safe. Rewriting c.76_77insG to c.76dup is a normalization decision requiring sequence data — that belongs to hgvs Normalizer, not a regex.

Never strip an accession version to make a lookup succeed. If UTA lacks that version, that is information.

Step 2: parse, validate, normalize

def process_one(s):
    out = {"cleaned": s, "parsed": None, "normalized": None,
           "status": None, "error": None}

    if not s:
        out["status"] = "empty"; return out

    try:
        var = hp.parse_hgvs_variant(s)
        out["parsed"] = str(var)
    except HGVSParseError as e:
        out["status"] = "parse_error"; out["error"] = str(e); return out

    try:
        hv.validate(var)
    except HGVSInvalidVariantError as e:
        out["status"] = "ref_mismatch"; out["error"] = str(e); return out
    except HGVSInvalidIntervalError as e:
        out["status"] = "out_of_range"; out["error"] = str(e); return out
    except HGVSDataNotAvailableError as e:
        out["status"] = "no_data"; out["error"] = str(e); return out
    except HGVSError as e:
        out["status"] = "invalid"; out["error"] = str(e); return out

    try:
        norm = hn.normalize(var)
        out["normalized"] = str(norm)
        out["status"] = "ok"
    except HGVSError as e:
        out["status"] = "normalize_error"; out["error"] = str(e)

    return out

Classify by exception type. "1,200 variants failed" tells you nothing; "1,150 no_data, 40 ref_mismatch, 10 out_of_range" tells you there are three separate problems with three separate fixes.

Step 3: run the batch

def normalize_table(df, col="variant", default_ac=None):
    cleaned = df[col].map(lambda s: tidy(s, default_ac))
    results = pd.DataFrame([process_one(s) for s in cleaned], index=df.index)
    out = df.join(results)

    print(out["status"].value_counts().to_string())
    ok = out["status"] == "ok"
    print(f"\nusable: {ok.sum():,}/{len(out):,} ({ok.mean():.1%})")

    changed = ok & (out["parsed"] != out["normalized"])
    print(f"changed on normalization: {changed.sum():,} ({changed.mean():.1%})")

    for st in out.loc[~ok, "status"].dropna().unique():
        print(f"\n--- {st} (examples) ---")
        print(out.loc[out["status"] == st, [col, "error"]].head(3).to_string())

    return out

result = normalize_table(variants, col="hgvs_c")

The "changed on normalization" count is diagnostic. A high proportion means whoever produced the file was not normalizing — worth telling them, and worth checking whether their downstream lookups were failing too.

Step 4: deduplicate on the normalized form

ok = result.query("status == 'ok'")

print(f"unique raw:        {ok['hgvs_c'].nunique():,}")
print(f"unique normalized: {ok['normalized'].nunique():,}")

collapsed = (ok.groupby("normalized")
               .agg(n_records=("hgvs_c", "size"),
                    raw_forms=("hgvs_c", lambda s: sorted(set(s))),
                    samples=("sample_id", lambda s: sorted(set(s))))
               .reset_index()
               .sort_values("n_records", ascending=False))

multi = collapsed[collapsed["raw_forms"].str.len() > 1]
print(f"\n{len(multi)} variants had multiple raw representations:")
print(multi.head()[["normalized", "raw_forms"]].to_string())

This is the payoff. Two labs describing the same deletion as c.100del and c.103del collapse into one row. Without normalization your cohort counts are wrong — you split one allele across several rows and understate its frequency.

Step 5: project and enrich

def project(norm_str):
    out = {"g": None, "p": None, "transcript": None, "proj_error": None}
    try:
        var = hp.parse_hgvs_variant(norm_str)
        out["transcript"] = var.ac
        out["g"] = str(am.c_to_g(var))
        try:
            out["p"] = str(am.c_to_p(var))
        except HGVSError as e:
            out["p"] = "p.?"; out["proj_error"] = f"protein: {e}"
    except HGVSError as e:
        out["proj_error"] = str(e)
    return out

proj = pd.DataFrame([project(s) for s in collapsed["normalized"]],
                    index=collapsed.index)
final = collapsed.join(proj)

Protein projection failing (intronic variants, complex indels) should not kill the genomic projection — hence the nested try.

Step 6: write with provenance

from datetime import datetime, timezone

final.attrs = PROVENANCE
final.to_parquet("normalized_variants.parquet")

meta = {
    **PROVENANCE,
    "run_at": datetime.now(timezone.utc).isoformat(),
    "n_input": len(variants),
    "n_ok": int((result["status"] == "ok").sum()),
    "n_unique": int(final["normalized"].nunique()),
    "status_counts": result["status"].value_counts().to_dict(),
}
Path("normalized_variants.meta.json").write_text(json.dumps(meta, indent=2))

An analysis that does not record its UTA and SeqRepo versions is not reproducible, regardless of how carefully you pinned your Python packages. See Applied - Reproducible Environment.

Handling the hard cases

Protein-only variants (p.Pro72Arg). You cannot recover the DNA change — several codon changes give the same amino acid substitution. Flag them and ask for the c. form.

Legacy IVS notation (IVS3+2T>C). Requires knowing the exon structure to convert. hgvs will not parse it. Either map it via the transcript's exon table, or escalate.

Gene symbol + protein shorthand (TP53 P72R). Requires a transcript choice and a reverse translation. Ambiguous by construction. Escalate rather than guess.

Transcript not in UTA. Check which versions exist before deciding:

avail = {t[3] for t in hdp.get_tx_for_gene(gene)}
print(sorted(a for a in avail if a.startswith(base_accession)))

If you substitute a different version, record that you did — do not silently swap.

The regression test

KNOWN = [
    ("NM_000546.6:c.215C>G",       "NC_000017.11:g.7676154G>C"),
    ("NM_000546.6:c.215C>G ",      "NC_000017.11:g.7676154G>C"),   # whitespace
    ("NM_000546.6:c.215C–G",  None),                          # en-dash → error
]

def test_pipeline():
    for raw, expected_g in KNOWN:
        r = process_one(tidy(raw))
        if expected_g:
            assert r["status"] == "ok", f"{raw}: {r['error']}"
            assert str(am.c_to_g(hp.parse_hgvs_variant(r["normalized"]))) == expected_g

Run it whenever you change hgvs version, UTA snapshot, SeqRepo snapshot, or assembly. Configuration errors that would otherwise reach a clinical report show up in seconds.

Performance

Step	Speed
Parse	~1,000–10,000/s (pure computation)
Validate	fast with local SeqRepo; slow remotely
Normalize	fast locally
Project	fast with local UTA; very slow remotely

Run UTA and SeqRepo locally for any batch work. The public instances are a shared courtesy service. See UTA and Data Providers and SeqRepo.

Common mistakes

Skipping normalization, so equivalent variants stay distinct.
Rewriting semantics during cleanup.
Stripping accession versions.
Catching bare Exception, losing the classification.
Not reporting failure counts and examples.
Deduplicating on the raw string.
Not recording UTA/SeqRepo versions.
The public UTA instance in a pipeline.
Guessing at ambiguous input rather than escalating it.

Applied - Variant Normalization Pipeline

Applied - Variant Normalization Pipeline

The problem

Setup

Step 1: safe cleanup

Step 2: parse, validate, normalize

Step 3: run the batch

Step 4: deduplicate on the normalized form

Step 5: project and enrich

Step 6: write with provenance

Handling the hard cases

The regression test

Performance

Common mistakes

See also