HGVS Nomenclature09/11
Starting sandbox…
Beginner·4 graded tests

HGVS Pitfalls

hgvs Pitfalls

Core Principle: Clinical variant curation requires vigilance against common nomenclature pitfalls, including reference sequence mismatch, transcript version ambiguity, 5'/3' indel shifting divergence, and pseudoautosomal region mappings.

Read this once now, and again after you have built a pipeline. It lands differently the second time.

1. Not normalizing

The single biggest source of variant-matching failures.

An indel in a repeat has several valid descriptions. Without normalization, c.100del and c.103del are different strings for the same event. Your ClinVar lookup returns nothing, which reads as "novel variant" rather than "I asked the wrong question".

var = hn.normalize(hp.parse_hgvs_variant(s))    # always

Normalize both sides of any comparison, in the same coordinate system. See hgvs Normalizer.

2. HGVS and VCF shift in opposite directions

Standard Convention
HGVS 3'-most on the reference sequence
VCF leftmost (5'-most on the genome)

The same deletion is POS=100 in a VCF and c.103del in HGVS. Both correct. Converting between them without renormalizing at each step loses variants silently.

For VCF output, use a 5'-shifting normalizer:

hn5 = hgvs.normalizer.Normalizer(hdp, shuffle_direction=5)

3. Assembly confusion

NC_000017.10:g.7579472G>C     GRCh37
NC_000017.11:g.7676154G>C     GRCh38

Same variant, ~97,000 bp apart. Mixing assemblies within one analysis produces coordinates that are simply wrong, and nothing errors.

  • Track the assembly explicitly, in a column, in every table.
  • Use NC_ accessions, never chr17: the accession version encodes the assembly.
  • Validation catches only ~75% of assembly mixups, because a wrong position matches the claimed ref base by chance about a quarter of the time.
  • Transcript projection is not liftover. Use CrossMap or liftOver with official chain files.

4. Accession versions matter

NM_000546.5 and NM_000546.6 can number the same base differently: a UTR revision shifts every c. coordinate downstream of it.

Never strip a version to make a lookup succeed. If UTA lacks the version you have, that is information: your data and your UTA snapshot disagree. Record the substitution explicitly if you make one.

5. There is no "the" transcript

A gene has many transcripts. A genomic variant has a different c. and p. description in each, sometimes with different consequence classes.

  • Use MANE Select as the default.
  • Respect gene-specific clinical conventions where they exist.
  • Record the transcript in every output row.
  • Consider reporting all consequences with MANE flagged.

6. p.(...) is a prediction

The parentheses are HGVS's marker for inferred rather than observed. c_to_p applies the codon table and knows nothing about nonsense-mediated decay, alternative splicing, readthrough, or reinitiation.

Preserve the parentheses. Do not present a predicted protein change as an observation.

7. Minus-strand genes shift the other way

"3'-most" is relative to the reference sequence in use. For a minus-strand gene, shifting 3' on the transcript means moving to lower genomic coordinates.

Consequence: normalize after projecting, in the target coordinate system. Normalizing once at the start and assuming it carries through is wrong.

8. Intronic offsets are not coordinates

c.87+3 and c.87 are different positions. c.12 and c.*12 are different positions. Code that reads posedit.pos.start.base without also checking offset and datum is wrong.

pos.base; pos.offset; pos.datum       # all three, always

See SequenceVariant.

9. Constructing objects in a loop

# catastrophic
for s in variants:
    hp = hgvs.parser.Parser()                                    # compiles a grammar
    hdp = hgvs.dataproviders.uta.connect()                        # opens a DB connection
    am = hgvs.assemblymapper.AssemblyMapper(hdp, ...)             # caches transcripts

Build once, reuse. This alone is often a 100× difference.

10. The public UTA instance in production

It is a shared courtesy service. Batch work through it is slow, unreliable, and antisocial. Run UTA and SeqRepo locally for anything beyond exploration. See UTA and Data Providers and SeqRepo.

11. Regex-parsing HGVS

re.match(r"c\.(\d+)([ACGT])>([ACGT])", s)     # handles ~40% of real variants

Fails on indels, intronic offsets, UTR positions, repeats, delins, protein notation, and phased alleles. Worse, it partially matches malformed input rather than rejecting it.

String manipulation of HGVS is acceptable only for coarse filtering (which accession? is it coding?). Never for semantics. Use hgvs Parser.

12. Insertions are not a coordinate shift

HGVS c.76_77insT names the flanking bases. VCF anchors on the preceding base and includes it in both alleles (REF=A ALT=AT). BED uses a zero-width interval.

These are different representations, not different coordinate origins. Converting requires a library, not arithmetic.

13. ins where dup is required

If the inserted sequence duplicates the immediately preceding bases, HGVS requires dup. Normalization fixes this: another reason to normalize.

14. Losing phase

c.[215C>G];[847C>T] (in trans) and c.[215C>G;847C>T] (in cis) are clinically different: two LoF variants in trans in a recessive gene is a diagnosis; in cis it is carrier status. Most pipelines discard the brackets.

15. Silent failures in batch processing

[ ]python
out1ms

This code completed without producing output.

Classify errors by type, count them, and inspect the patterns. Failures cluster around systematic causes; the pattern is the diagnosis.

16. Not recording provenance

Three version numbers determine your output and all three must be in it:

print(f"hgvs={hgvs.__version__} "
      f"uta={hdp.data_version()} "
      f"seqrepo={os.environ.get('HGVS_SEQREPO_DIR')} "
      f"assembly={am.assembly_name} "
      f"aln_method=splign")

Pinning your Python packages is not enough: the data versions change the answers. See Applied - Reproducible Environment.

17. Confusing validation with interpretation

A valid variant is one that exists as described. It says nothing about pathogenicity, sequencing quality, or clinical relevance.

18. Pseudoautosomal regions

PAR loci exist identically on both X and Y. Without in_par_assume="X", projection raises. With it, you have made an assumption: record it.

The defensive checklist

Before trusting variant output, confirm:

  • Every variant parsed, or the failure was recorded
  • Every variant validated against the reference sequence
  • Every variant normalized, in its final coordinate system
  • The assembly is tracked in a column, not assumed
  • Accession versions are preserved
  • The transcript used is recorded per row
  • p. predictions are labelled as predictions
  • hgvs, UTA and SeqRepo versions are in the output
  • Error counts are reported by category
  • A known-answer set round-trips correctly

A regression test worth writing

[ ]python
out1ms

This code completed without producing output.

Run it whenever you change UTA snapshot, SeqRepo snapshot, hgvs version, or assembly. Two well-known variants take seconds to check and catch configuration errors that would otherwise reach a report.

See also

biocommons hgvs · hgvs Normalizer · hgvs Validator · AssemblyMapper · HGVS Nomenclature · Coordinate Systems · Applied - Variant Normalization Pipeline · Applied - Reproducible Environment

Exercise

The hgvs package is not available in this sandbox, and that suits this exercise: every problem below is visible in the string itself, before any database is consulted. Writing the audit yourself gives you the cheap pre-flight check worth running over a variant column before it reaches the parser.

Write audit(lines), taking the input lines and returning a list of output lines. Ignore blank lines. Emit <variant> <issue> for each remaining line, where <issue> is the first of these that applies, or ok if none does:

  • legacy_ivs: the string uses IVS intron numbering
  • missing_accession: there is no : at all
  • chr_reference: the reference is written chr17 rather than an accession
  • unversioned_accession: the accession carries no .<version>
  • mixed_assembly: the reference is an NC_ accession whose version differs from the majority NC_ version in this batch, breaking ties towards the higher version number
  • unmarked_prediction: the description is a p. change that is not wrapped in parentheses, so a prediction is being presented as an observation

Finish with checked=<n>, then one <issue>=<count> line per issue that occurred, in alphabetical order. checked=0 alone is the whole output when the input is empty. Note that mixed_assembly needs the whole batch before any one line can be judged, which is the point: assembly confusion is a property of a file, not of a variant.

exerciseEvery submission runs independently in a fresh namespace.

Test cases · 4

#viainputexpected stdout
1stdinNC_000017.11:g.7676154G>C NC_000017.11:g.7676155A>T NC_000013.11:g.32340300G>A NC_000017.10:g.7579472G>C chr17:g.7676154G>C IVS3+2T>C c.215C>G NM_000546:c.215C>G NP_000537.3:p.Pro72Arg NP_000537.3:p.(Pro72Arg) NC_000017.11:g.7676154G>C ok NC_000017.11:g.7676155A>T ok NC_000013.11:g.32340300G>A ok NC_000017.10:g.7579472G>C mixed_assembly chr17:g.7676154G>C chr_reference IVS3+2T>C legacy_ivs c.215C>G missing_accession NM_000546:c.215C>G unversioned_accession NP_000537.3:p.Pro72Arg unmarked_prediction NP_000537.3:p.(Pro72Arg) ok checked=10 chr_reference=1 legacy_ivs=1 missing_accession=1 mixed_assembly=1 ok=4 unmarked_prediction=1 unversioned_accession=1
2stdinNC_000017.10:g.7579472G>C NC_000013.10:g.32912299G>A NC_000017.11:g.7676154G>C NM_000546.6:c.215C>G NC_000017.10:g.7579472G>C ok NC_000013.10:g.32912299G>A ok NC_000017.11:g.7676154G>C mixed_assembly NM_000546.6:c.215C>G ok checked=4 mixed_assembly=1 ok=3
3stdinNM_000546.6:c.215C>G NM_000492.4:c.1521_1523delCTT NP_000483.3:p.(Phe508del) NM_000546.6:c.215C>G ok NM_000492.4:c.1521_1523delCTT ok NP_000483.3:p.(Phe508del) ok checked=3 ok=3
4stdin checked=0

Hints · 3

01Hint

Take two passes: tally the NC_ accession versions first, then classify each line against the majority.

02Hint

max(counts, key=lambda v: (counts[v], int(v))) resolves a tie towards the higher version number.

03Hint

A p. description is marked as predicted when it matches p\.\(.+\) in full.