Projection c to g to p

Core Principle: Coordinate projection maps sequence variants between genomic (g.), transcript (c.), and protein (p.) spaces, accounting for transcript alignments, strand orientation, and translation rules.

The coordinate systems

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.

'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

Exercise

The hgvs package is not available in this sandbox, so instead of calling am.c_to_g you will walk an exon table supplied on stdin. That table is the whole of what UTA provides for this operation, and stepping through it by hand is what makes the strand rule concrete rather than memorised.

Write project(text), taking the whole of stdin and returning a list of output lines. Above a line -- stdin gives strand + or strand -, genomic <accession>, cds_start <n> and cds_end <n> as transcript positions, and one exon <n_start> <n_end> <g_start> <g_end> line per exon, where g_start is the genomic coordinate of n_start and runs downwards on the minus strand. Below the -- come c. substitutions only. Convert the coding position to a transcript position: c.n is cds_start + n - 1, c.-n is cds_start - n, c.*n is cds_end + n: locate its exon, and step along the exon in the strand's direction. Apply any intronic offset in the same direction afterwards, and on the minus strand complement both the reference and the alternate base.

Emit <genomic_accession>:g.<pos><ref>><alt> for each variant, <variant> out_of_range when the transcript position falls in no exon, and <variant> unsupported for anything that is not a substitution. Finish with variants=<n>, which is the only output when no variants follow the --. The lesson's worked example is the check to keep in mind: NM_000546.6:c.215C>G projects to NC_000017.11:g.7676154G>C, the bases complemented because TP53 is on the minus strand.