Applied - Variant Normalization Pipeline

Core Principle: Clinical variant normalization pipelines standardize heterogeneous variant descriptions through syntax validation, reference allele verification, canonical 3' indel shifting, and cross-assembly projection.

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>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

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.

See also

hgvs Normalizer · hgvs Validator · hgvs Parser · hgvs Pitfalls · Applied - Annotating a Variant Table · UTA and Data Providers · Applied - Reproducible Environment

Exercise

hgvs is not installed in this sandbox, so implement the two stages the lesson delegates to it: the safe formatting cleanup, and the normalization itself, which is left-shifting and trimming against a reference sequence. Doing it by hand is what shows why a regex must never touch semantics.

Write normalize(text), returning a list of output lines. The first input line is the reference sequence, 1-based. Every later non-blank line is one raw variant record in the form pos:ref>alt, possibly carrying formatting noise: surrounding whitespace, internal whitespace, an HTML-escaped > in place of the >, an en-dash – used as the separator, an accession prefix such as NC_000017.11:g. that must be stripped, and lower-case bases. Clean all of that, then uppercase. A record that still does not parse, that leaves either allele empty, whose alleles contain anything but A, C, G and T, or whose two alleles are identical and so describe no change at all, is a parse_error. A record whose reference allele runs past the end of the reference sequence, or whose position is below one, is out_of_range. A record whose reference allele does not match the reference sequence at that position is a ref_mismatch.

Normalise the rest. Repeatedly, while both alleles are non-empty and end in the same base, drop that base from both; if this empties either allele, step the position back by one and prepend the reference base at the new position to both. Should a trim empty an allele while the position is already one, the variant has no preceding base to anchor to: count it as out_of_range and drop it. Then, while both alleles are longer than one base and start with the same base, drop that base from both and advance the position by one. Print input=<n>, then ok=<n> parse_error=<n> ref_mismatch=<n> out_of_range=<n>, then unique=<n> duplicates=<n> where duplicates are the successfully normalised records beyond the first for each distinct result. Then one line per distinct normalised variant, sorted by position, then reference allele, then alternate allele, holding the position, the two alleles and the number of input records that produced it, tab separated. If nothing normalises, print no variants normalized after the counts instead.