AssemblyMapper

Core Principle: AssemblyMapper provides coordinate projection of sequence variants across genomic (g.), coding transcript (c.), non-coding transcript (n.), and protein (p.) reference spaces for a designated genome assembly.

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

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

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

Exercise

The hgvs package is not available in this sandbox, so you will implement relevant_transcripts and the transcript-selection policy yourself against a small table supplied on stdin. Both halves of what AssemblyMapper adds over VariantMapper: restricting to one assembly, and refusing to pick a transcript silently: become visible when you write them out.

Write map_variants(text), taking the whole of stdin and returning a list of output lines. Above a line -- stdin gives assembly <name> for the assembly this mapper targets, one tx <transcript> <assembly> <genomic_accession> <start> <end> line per transcript–genome alignment, and zero or more mane <transcript> lines. Below it come g. variants. For each variant emit <variant> wrong_assembly when its genomic accession appears in the table only under some other assembly, <variant> intergenic when no transcript on the target assembly contains the position, <variant> unsupported when the string is not a g. variant with a numeric position, and otherwise <variant> relevant=<n> chosen=<transcript>, where n counts the overlapping transcripts and the chosen one is the MANE Select transcript if it is among them and the alphabetically first otherwise. Finish with projected=<n> counting the variants that got a transcript; that line is the only output when no variants follow the --.