Projection c to g to p

In one line: moving a variant between coordinate systems — the operation that makes a lab's c. variant comparable to a pipeline's g. variant.

The coordinate systems

        ┌─────────────────────────────────────────────┐
        │  g.   genomic — NC_000017.11:g.7676154G>C   │
        │       one position, assembly-specific        │
        └──────────────────┬──────────────────────────┘
                           │  needs the transcript's exon structure,
                           │  strand, and CDS boundaries
        ┌──────────────────▼──────────────────────────┐
        │  c.   coding — NM_000546.6:c.215C>G          │
        │       ONE PER TRANSCRIPT — a gene has many   │
        └──────────────────┬──────────────────────────┘
                           │  needs the codon table
        ┌──────────────────▼──────────────────────────┐
        │  p.   protein — NP_000537.3:p.(Pro72Arg)     │
        │       a PREDICTION                            │
        └─────────────────────────────────────────────┘

Also: n. for non-coding transcripts, m. for mitochondrial, r. for RNA.

The asymmetry that matters

c. → g. is one-to-one. A coding position on a specific transcript maps to exactly one genomic position.

g. → c. is one-to-many. A genomic position falls in however many transcripts overlap it. For a gene with 20 annotated transcripts, one genomic variant has 20 different c. descriptions, several of which may have different predicted consequences — missense in one, intronic in another, absent from a third.

am.c_to_g(var_c)                       # unambiguous
am.relevant_transcripts(var_g)          # → a LIST
am.g_to_c(var_g, "NM_000546.6")         # you must choose

There is no "the" coding description of a genomic variant. Choosing a transcript is a policy decision. See AssemblyMapper.

The full pipeline

import hgvs.parser, hgvs.dataproviders.uta, hgvs.assemblymapper
import hgvs.normalizer, hgvs.validator
from hgvs.exceptions import HGVSError

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

def annotate(hgvs_c):
    var_c = hp.parse_hgvs_variant(hgvs_c)     # 1. parse
    hv.validate(var_c)                         # 2. validate
    var_c = hn.normalize(var_c)                # 3. normalize
    var_g = am.c_to_g(var_c)                   # 4. project to genome
    var_p = am.c_to_p(var_c)                   # 5. predict protein
    return {"c": str(var_c), "g": str(var_g), "p": str(var_p),
            "transcript": var_c.ac, "assembly": am.assembly_name}

From genomic (e.g. a VCF)

def annotate_genomic(chrom_ac, pos, ref, alt, gene=None):
    var_g = hp.parse_hgvs_variant(f"{chrom_ac}:g.{pos}{ref}>{alt}")
    hv.validate(var_g)
    var_g = hn.normalize(var_g)

    txs = am.relevant_transcripts(var_g)
    if not txs:
        return {"g": str(var_g), "note": "intergenic"}

    tx = preferred_transcript(gene, txs)       # your MANE-based policy
    var_c = am.g_to_c(var_g, tx)
    var_p = am.c_to_p(var_c)
    return {"g": str(var_g), "c": str(var_c), "p": str(var_p), "transcript": tx}

Chromosome accession, not chr17. HGVS needs NC_000017.11 (GRCh38) or NC_000017.10 (GRCh37). The accession encodes the assembly, which is exactly the property chr17 lacks.

from bioutils.assemblies import make_ac_name_map
AC_BY_NAME = {v: k for k, v in make_ac_name_map("GRCh38").items()}
AC_BY_NAME["17"]        # 'NC_000017.11'

bioutils (a biocommons package, already a dependency) provides these maps so you do not hardcode them.

What can go wrong at each step

Step	Failure	Meaning
parse	HGVSParseError	malformed string
validate	HGVSInvalidVariantError	ref base mismatch → wrong assembly/version/strand
validate	HGVSInvalidIntervalError	position outside the sequence
normalize	HGVSDataNotAvailableError	sequence missing from SeqRepo
c_to_g	HGVSDataNotAvailableError	transcript not in UTA
c_to_g	HGVSInvalidIntervalError	position beyond the transcript
c_to_p	HGVSUnsupportedOperationError	protein prediction not possible for this edit type
any	HGVSUsageError	wrong method for the variant type (e.g. c_to_g on a g. variant)

Intronic variants

var_c = hp.parse_hgvs_variant("NM_000546.6:c.-28-2A>G")   # splice acceptor
var_g = am.c_to_g(var_c)                                    # works fine
var_p = am.c_to_p(var_c)                                    # p.? — unpredictable

Intronic variants project to genomic coordinates without difficulty. Protein prediction cannot work, because the effect depends on splicing, which the codon table knows nothing about. The result is p.?.

For splice-effect prediction you need a dedicated tool — SpliceAI, MaxEntScan, or Pangolin. hgvs will tell you where the variant is relative to the splice site (via the offset — see SequenceVariant), which is useful input to that decision, but it will not predict the consequence.

Frameshift prediction

var_c = hp.parse_hgvs_variant("NM_000546.6:c.215del")
am.c_to_p(var_c)          # NP_000537.3:p.(Pro72ArgfsTer13)

fsTer13 means "the new reading frame hits a stop 13 residues later". This is a mechanical calculation over the reference sequence — correct as arithmetic, and it does not tell you whether the transcript is actually degraded by nonsense-mediated decay (which depends on where the new stop sits relative to the last exon junction) or whether a downstream reinitiation produces a partial protein.

Round-tripping

var_c  = hp.parse_hgvs_variant("NM_000546.6:c.215C>G")
var_g  = am.c_to_g(var_c)
var_c2 = am.g_to_c(var_g, var_c.ac)
str(var_c) == str(var_c2)          # True — IF both were normalized

Round-tripping is a good self-test in a pipeline. It should be lossless for a normalized variant on a transcript with a clean alignment. Where it is not lossless, you have found something interesting: usually a transcript–genome alignment gap, which is a region where the transcript and reference genome genuinely disagree.

Multiple transcripts, reported honestly

def all_consequences(var_g, am):
    rows = []
    for tx in am.relevant_transcripts(var_g):
        try:
            var_c = am.g_to_c(var_g, tx)
            var_p = am.c_to_p(var_c)
            rows.append({"transcript": tx, "c": str(var_c), "p": str(var_p),
                         "mane": tx in MANE_TRANSCRIPTS})
        except HGVSError as e:
            rows.append({"transcript": tx, "error": str(e)})
    return pd.DataFrame(rows).sort_values("mane", ascending=False)

Reporting all consequences with the MANE transcript flagged is more honest than silently picking one, and it is what a clinical reviewer will want to see when the MANE consequence is benign but another transcript's is not.

Common mistakes

• Expecting g_to_c to pick a transcript for you.
• chr17 instead of NC_000017.11.
• Mixing assemblies.
• Treating p.(...) as observed.
• Expecting protein prediction for intronic variants.
• Not normalizing before projecting.
• Constructing AssemblyMapper per variant.
• Reporting one transcript's consequence without noting that others differ.
• Not recording the transcript, assembly and data versions in the output.

See also

AssemblyMapper · SequenceVariant · hgvs Normalizer · hgvs Validator · HGVS Nomenclature · Applied - Annotating a Variant Table · hgvs Pitfalls