pythonforbio.
[WASM idle]
HGVS Nomenclature01/11

AssemblyMapper

Starting sandbox…
Beginnerlesson

AssemblyMapper

In one line: the high-level projector — moves a variant between genomic, transcript and protein coordinates for a chosen assembly, and normalizes as it goes.

Construction

import hgvs.assemblymapper, hgvs.dataproviders.uta

hdp = hgvs.dataproviders.uta.connect(pooling=True)

am = hgvs.assemblymapper.AssemblyMapper(
    hdp,
    assembly_name="GRCh38",        # or "GRCh37"
    alt_aln_method="splign",       # for RefSeq transcripts
    normalize=True,                 # normalize output (default True)
    replace_reference=True,         # fill in ref bases from the sequence
    prevalidation_level="INTRINSIC",  # None | INTRINSIC | EXTRINSIC
    in_par_assume="X",              # pseudoautosomal region → assume X
)

Build it once. Each construction opens database connections and caches transcript data. Constructing it per variant is the most common performance mistake in hgvs code.

The parameters worth understanding:

  • replace_reference=True fills in the reference bases from the actual sequence. This turns c.76_78del into c.76_78delCTT and — more importantly — surfaces mismatches between what the variant claims and what the sequence says.
  • normalize=True applies 3'-shifting to the output. Almost always what you want. See hgvs Normalizer.
  • in_par_assume resolves the pseudoautosomal regions, where a locus exists on both X and Y with identical sequence. Without it, PAR variants raise.

The projection methods

am.c_to_g(var_c)                 # transcript coding → genomic
am.g_to_c(var_g, tx_ac)          # genomic → a SPECIFIC transcript
am.c_to_p(var_c)                 # coding → protein (a PREDICTION)
am.c_to_n(var_c); am.n_to_c(var_n)
am.g_to_n(var_g, tx_ac); am.n_to_g(var_n)
am.g_to_t(var_g, tx_ac)          # dispatches to g_to_c or g_to_n
am.t_to_g(var_t)
am.t_to_p(var_t)

am.relevant_transcripts(var_g)   # which transcripts overlap this genomic variant?
am.assembly_name; am.data_version()

g_to_c requires a transcript accession — that is the point. A genomic variant has one position; its coding consequence depends entirely on which transcript you ask about, and a gene can have twenty. There is no "the" transcript.

The transcript selection problem

var_g = hp.parse_hgvs_variant("NC_000017.11:g.7676154G>C")
am.relevant_transcripts(var_g)
# ['NM_000546.6', 'NM_001126112.3', 'NM_001126113.3', ...]

Which one do you report? This is a policy decision, not a technical one, and it must be made explicitly:

Strategy Rationale
MANE Select current best practice; one agreed transcript per gene, RefSeq+Ensembl matched
RefSeq Select NCBI's representative transcript
Longest CDS crude but deterministic
Clinically established some genes have a legacy transcript everyone uses (e.g. specific BRCA1 and CFTR transcripts)
All of them most complete; report every consequence

MANE Select is the right default in 2026. Where a clinical convention exists for a gene, it overrides MANE — because a variant reported against a different transcript than the literature uses will not match anything.

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

def preferred_tx(gene, available):
    mane = MANE_BY_GENE.get(gene)
    if mane in available:
        return mane
    return sorted(available)[0] if available else None

Record which transcript you used in every output row. A c. variant without its accession is not interpretable.

GRCh37 vs GRCh38

am37 = AssemblyMapper(hdp, assembly_name="GRCh37", alt_aln_method="splign")
am38 = AssemblyMapper(hdp, assembly_name="GRCh38", alt_aln_method="splign")

am37.c_to_g(var_c)      # NC_000017.10:g.7579472G>C
am38.c_to_g(var_c)      # NC_000017.11:g.7676154G>C

Same variant, ~97,000 bp apart. A genomic coordinate without an assembly is meaningless, and mixing assemblies within one analysis is a silent, serious error.

Note the accession itself encodes the assembly: NC_000017.10 is GRCh37 chr17, NC_000017.11 is GRCh38 chr17. That versioning is a genuine safety feature — unlike chr17, which tells you nothing.

Projecting through a transcript is not liftover. Going g.(GRCh37) → c. → g.(GRCh38) works only where the transcript aligns to both assemblies consistently. For genuine coordinate conversion, use CrossMap or UCSC liftOver with the official chain files.

Protein projection is a prediction

var_p = am.c_to_p(var_c)
str(var_p)          # 'NP_000537.3:p.(Pro72Arg)'

The parentheses are HGVS's marker for a predicted consequence. c_to_p applies the reference codon table to the reference sequence; it does not know about:

  • nonsense-mediated decay
  • alternative splicing that skips the exon
  • the actual translated product
  • readthrough or reinitiation

For a frameshift it will predict a new stop position, which is a mechanical calculation, not an observation. Treat these as annotations, not facts, and preserve the parentheses in your reports.

Error handling

from hgvs.exceptions import (
    HGVSError, HGVSDataNotAvailableError, HGVSInvalidVariantError,
    HGVSInvalidIntervalError, HGVSUsageError, HGVSUnsupportedOperationError,
)

def project(am, var_c):
    try:
        return am.c_to_g(var_c), None
    except HGVSDataNotAvailableError as e:
        return None, f"no transcript data: {e}"       # UTA lacks this accession
    except HGVSInvalidIntervalError as e:
        return None, f"position out of range: {e}"     # position beyond the transcript
    except HGVSInvalidVariantError as e:
        return None, f"invalid variant: {e}"           # ref base mismatch
    except HGVSError as e:
        return None, f"hgvs error: {e}"

Catch the specific exception types — each means a different data problem, and lumping them together loses the diagnostic information you need.

HGVSUnsupportedOperationError from c_to_p typically means a variant type the protein predictor cannot handle (some complex indels, variants spanning the start codon).

Batch pattern

def project_all(df, am, hp):
    rows = []
    for r in df.itertuples():
        row = {"input": r.hgvs, "g": None, "p": None, "error": None}
        try:
            v = hp.parse_hgvs_variant(r.hgvs)
            row["g"] = str(am.c_to_g(v))
            row["p"] = str(am.c_to_p(v))
        except HGVSError as e:
            row["error"] = f"{type(e).__name__}: {e}"
        rows.append(row)
    out = pd.DataFrame(rows)
    print(out["error"].notna().value_counts())
    return out

Never let one bad variant kill a batch of 50,000. Collect the errors, report the counts, and inspect the patterns — they usually cluster around one systematic cause.

AssemblyMapper vs VariantMapper

import hgvs.variantmapper
vm = hgvs.variantmapper.VariantMapper(hdp)
vm.c_to_g(var_c, "NC_000017.11")        # you supply the genomic accession

VariantMapper is the lower-level engine; AssemblyMapper wraps it with assembly awareness, transcript filtering and normalization. Use AssemblyMapper unless you need the extra control.

Common mistakes

  • Constructing AssemblyMapper per variant.
  • Not recording which transcript was used.
  • Mixing GRCh37 and GRCh38.
  • Treating p.(...) as observed.
  • Using transcript projection as liftover.
  • Catching bare Exception, losing the diagnostic type.
  • Assuming one transcript per gene.
  • Forgetting replace_reference=True, so ref-base mismatches go undetected.
  • Ignoring PAR regions on X/Y.

See also

Projection c to g to p · UTA and Data Providers · hgvs Normalizer · hgvs Validator · hgvs Pitfalls · Applied - Annotating a Variant Table

lesson example

No output yet — run the code to populate this drawer.

start of pathHGVS Nomenclature