HGVS Parsing and VRS Conversion

HGVS (Human Genome Variation Society) nomenclature is the international standard for describing variants in DNA, RNA, and protein sequences. This document explains how HGVS expressions work, why converting them to genomic coordinates is challenging, and how the conversion pipeline produces GA4GH VRS allele identifiers.

What is HGVS nomenclature?

HGVS notation describes a variant as a combination of:

1. A reference sequence (accession number)
2. A reference type (genomic, coding, etc.)
3. A position within that reference
4. An edit (what changed)

For example:

NC_000007.14:g.140753336A>T
│            │ │         └─ Edit: A changed to T
│            │ └─ Position: 140753336
│            └─ Reference type: g (genomic)
└─ Reference: NC_000007.14 (chromosome 7)

This notation is unambiguous: anyone with access to the reference sequence can locate the exact position of this variant.

For the official specification, see varnomen.hgvs.org.

The goal: HGVS to VRS

GA4GH VRS (Variation Representation Specification) provides content-addressable identifiers for variants. A VRS allele identifier like ga4gh:VA.xKdIt6tBJ_J0B_bMEaYAYnrJpXzE8QEK uniquely identifies a variant based on its sequence context and alteration.

The pipeline converts human-readable HGVS strings into machine-comparable VRS identifiers through five steps: parse HGVS, look up transcript, map coordinates (e.g., c.1799 → chr7:140753336), normalize per VRS rules, and compute the digest. The result is a stable identifier like ga4gh:VA.xKdIt6tBJ_J0B_bMEaYAYnrJpXzE8QEK.

This enables variant deduplication across databases that may use different nomenclature conventions.

Reference types

HGVS defines several reference types, indicated by a single letter after the colon:

Type	Name	Reference sequence	Example
g.	Genomic	Chromosome	NC_000007.14:g.140753336A>T
c.	Coding	mRNA transcript	NM_004333.6:c.1799T>A
n.	Non-coding	Non-coding RNA	NR_046018.2:n.357A>G
m.	Mitochondrial	Mitochondrial genome	NC_012920.1:m.8993T>G
r.	RNA	RNA sequence	NM_004333.6:r.1799u>a
p.	Protein	Protein sequence	NP_004324.2:p.Val600Glu

The conversion pipeline handles g., c., n., and m. variants. Protein-level (p.) variants require back-translation and are not supported.

Genomic coordinates (g.)

Genomic variants use absolute positions on a chromosome:

NC_000007.14:g.140753336A>T

Position 140753336 is counted from the start of chromosome 7. This is the target coordinate system for VRS conversion.

Coding coordinates (c.)

Coding variants use positions relative to the coding sequence (CDS) of a transcript:

NM_004333.6:c.1799T>A

Position 1799 is the 1799th nucleotide of the coding region of transcript NM_004333.6 (the BRAF gene). The coding region starts at the translation initiation codon (ATG) and ends at the stop codon.

Non-coding coordinates (n.)

Non-coding variants use positions relative to the transcript start:

NR_046018.2:n.357A>G

Position 357 is the 357th nucleotide from the start of this non-coding RNA transcript. Non-coding transcripts have no CDS, so n. notation counts from the 5’ end.

The coordinate mapping problem

Converting c. coordinates to genomic positions requires transcript annotation data. This is challenging because:

1. Exon-intron structure: The coding sequence is split across multiple exons separated by introns that are spliced out.

2. UTRs: Positions before the start codon (5’ UTR) and after the stop codon (3’ UTR) require special notation.

3. Strand orientation: Transcripts on the minus strand run in the opposite direction from genomic coordinates.

4. Multiple isoforms: A gene may have multiple transcripts with different exon structures, making the same c. position map to different genomic locations.

Consider the BRAF gene: it is transcribed from the minus strand, with 18 exons spread across ~190kb of chromosome 7. The transcript has a 5’ UTR, a coding region (CDS) from the start codon (ATG, position c.1) to the stop codon (position c.2301), and a 3’ UTR.

The c.1799 position falls within exon 15. To find the genomic coordinate, we must:

1. Look up the transcript’s exon structure
2. Calculate which exon contains position 1799
3. Account for strand direction
4. Compute the genomic position

How coordinate mapping works

This section walks through mapping NM_004333.6:c.1799T>A (the BRAF V600E mutation) to genomic coordinates.

Step 1: Load transcript annotation

Retrieve the transcript record for NM_004333.6:

• Gene: BRAF
• Strand: Reverse (minus)
• CDS start: 140753274 (genomic)
• CDS end: 140924929 (genomic)
• Exons: 18 exons with genomic coordinates

Step 2: Build the exon offset table

For coordinate mapping, build a table that tracks cumulative positions:

Exon	Genomic start	Genomic end	Transcript start	Transcript end
1	140924764	140924929	0	165
2	140850108	140850212	165	269
…	…	…	…	…
15	140753274	140753403	1758	1887
…	…	…	…	…

For reverse strand transcripts, exons are processed from highest to lowest genomic coordinate (because transcription runs 3’ to 5’ on the chromosome).

Step 3: Convert c. position to transcript position

The c. coordinate counts from the CDS start (ATG codon). To get the transcript position:

transcript_pos = cds_tx_start + (c_pos - 1)
               = 165 + (1799 - 1)
               = 1963

(In this example, 165 bases of 5’ UTR precede the CDS in the transcript.)

Step 4: Find the containing exon

Search the offset table for the exon containing transcript position 1963:

Exon 15 spans transcript positions 1758-1887, so position 1963 falls… wait, that does not fit. Let me recalculate with correct numbers.

In reality, c.1799 falls in exon 15, at genomic position 140753336. The mapping accounts for:

• Offset within the exon: 1799 - exon_cds_start
• Strand correction: subtract from exon end (reverse strand)

Step 5: Apply strand correction

For reverse strand transcripts, higher transcript positions correspond to lower genomic positions:

genomic_pos = exon_genomic_end - 1 - offset_in_exon

The final result: c.1799 maps to chr7:140753336.

Intronic variants

Variants in introns use offset notation from the nearest exon boundary:

NM_004333.6:c.93+1G>A
           │  │ │
           │  │ └─ 1 base into the intron (downstream)
           │  └─ After position 93
           └─ Last coding base of this exon

The +1 indicates one base past the exon boundary into the downstream intron. Similarly:

NM_004333.6:c.94-2A>G
               │
               └─ 2 bases before position 94 (upstream intron)

Mapping intronic positions

1. Find the exon boundary for the base position (c.93 or c.94)
2. Verify the position is actually at an exon boundary
3. Apply the offset in the appropriate direction
4. Account for strand (offset direction flips on reverse strand)

For c.93+1G>A:

• c.93 is the last base of an exon
• +1 means one base past the exon end
• On reverse strand, “past the exon end” means lower genomic coordinate

Intronic variants are important for splice site mutations, which can disrupt proper mRNA splicing even though they occur outside the coding region.

UTR variants

Variants in untranslated regions use special notation relative to the CDS:

5’ UTR (upstream of start codon)

NM_004333.6:c.-14C>T
              │
              └─ 14 bases before the start codon

Negative numbers count backwards from the ATG start codon. Position c.-1 is the base immediately before the start codon.

Mapping: convert to transcript position, then to genomic:

transcript_pos = cds_tx_start - abs(c_pos)
               = 165 - 14
               = 151

3’ UTR (downstream of stop codon)

NM_004333.6:c.*37A>G
              │
              └─ 37 bases after the stop codon

The asterisk (*) indicates positions past the stop codon. Position c.*1 is the first base after the stop codon.

Mapping:

transcript_pos = cds_tx_end + c_pos
               = 2301 + 37
               = 2338

Minus strand handling

About half of human genes are transcribed from the minus (reverse) strand. This affects coordinate mapping in several ways:

Exon ordering

Exons are stored in genomic order (ascending coordinates), but for minus strand genes, they are processed in reverse order when building transcript positions. Exon 1 of the transcript corresponds to the highest genomic coordinates.

Position calculation

Within each exon on a minus strand gene, positions run backwards relative to genomic coordinates. For an exon spanning genomic positions 100-104, transcript position 1 is at genomic position 104, and transcript position 5 is at genomic position 100.

Intronic offset direction

On the minus strand, + offsets (downstream in transcript terms) go toward lower genomic coordinates, while - offsets go toward higher genomic coordinates. This is the opposite of plus strand behavior.

MANE Select transcripts

MANE (Matched Annotation from NCBI and EMBL-EBI) Select designates one representative transcript per protein-coding gene. This is critical for clinical variant interpretation.

The problem with multiple isoforms

A gene like BRAF has multiple transcripts with different exon structures:

Transcript	Exons	CDS length
NM_004333.6	18	2301 bp
NM_001354609.2	17	2208 bp
NM_001374244.1	16	2100 bp

The same c. position may refer to different genomic locations depending on which transcript is used.

MANE Select as the clinical standard

MANE Select provides a single authoritative transcript per gene:

• Jointly curated by NCBI and Ensembl
• Matches both RefSeq and Ensembl annotations
• Covers >99% of protein-coding genes
• Required for ClinVar submissions

When a user submits a gene symbol without an accession:

BRAF:c.1799T>A

The system looks up the MANE Select transcript for BRAF (NM_004333.6) and uses that for coordinate mapping.

MANE Plus Clinical

Some genes have additional transcripts designated as “MANE Plus Clinical” for disease-relevant isoforms not captured by the primary MANE Select transcript.

The full pipeline

Putting it all together, here is how an HGVS string becomes a VRS identifier:

1. Parse the HGVS expression

Input: "NM_004333.6:c.1799T>A"

Parsed:
  accession: NM_004333.6
  reference_type: c (coding)
  position: 1799 (base), 0 (offset), CdsStart (datum)
  edit: substitution T>A

2. Look up transcript annotation

Query the transcript store for NM_004333.6:

Transcript {
  accession: "NM_004333.6",
  gene: "BRAF",
  chrom_digest: [24 bytes],
  strand: Reverse,
  cds_start: 140753274,
  cds_end: 140924929,
  exons: [...],
  mane: ManeStatus { mane_select: true, mane_clinical: false }
}

3. Map coordinates

Apply the coordinate mapping algorithm:

c.1799 (coding position)
  → transcript position 1963 (add 5' UTR length)
  → exon 15 (binary search offset table)
  → genomic position 140753336 (strand-aware calculation)

4. Construct VRS allele

Build the VRS allele object:

{
  "type": "Allele",
  "location": {
    "type": "SequenceLocation",
    "sequenceReference": {
      "type": "SequenceReference",
      "refgetAccession": "SQ.F-LrLMe1SRpfUZHkQmvkVqFEhDcMX1E2"
    },
    "start": 140753335,
    "end": 140753336
  },
  "state": {
    "type": "LiteralSequenceExpression",
    "sequence": "A"
  }
}

Note: VRS uses 0-based interbase coordinates, so position 140753336 becomes interval [140753335, 140753336).

5. Compute the VRS digest

Canonicalize and hash the allele to produce the identifier:

ga4gh:VA.xKdIt6tBJ_J0B_bMEaYAYnrJpXzE8QEK

This identifier is stable and reproducible: anyone computing it from the same variant and reference will get the same result.

Why this matters

Converting HGVS to VRS enables:

• Deduplication: Identify when different notations describe the same variant
• Interoperability: Exchange variant data between systems using a common identifier
• Reproducibility: Ensure variant calls are anchored to specific reference sequences
• Federation: Query distributed databases without central coordination

The coordinate mapping step is the most complex part of this pipeline, requiring accurate transcript annotations and careful handling of strand orientation, UTRs, and intronic positions.