From 1a64d75b69a83b941c13f9bf4debfa1ffd74c7f2 Mon Sep 17 00:00:00 2001 From: Rakesh Pai <41351936+developer-rpai@users.noreply.github.com> Date: Sat, 19 Sep 2026 22:30:17 -0700 Subject: [PATCH] feat: support to_hgvs for single-residue protein variants (#633) Implements HGVS export for protein Alleles with a single-residue LiteralSequenceExpression state: missense/nonsense substitutions (e.g. NP_060204.1:p.Val261Ala) and single-residue deletions (e.g. NP_060204.1:p.Val261del). Multi-residue changes, insertions, and non-literal states raise ValueError. Adds tests/extras/test_to_hgvs_protein.py with 9 tests covering the issue example plus negative controls. --- src/ga4gh/vrs/utils/hgvs_tools.py | 65 ++++++++- tests/extras/test_to_hgvs_protein.py | 189 +++++++++++++++++++++++++++ 2 files changed, 250 insertions(+), 4 deletions(-) create mode 100644 tests/extras/test_to_hgvs_protein.py diff --git a/src/ga4gh/vrs/utils/hgvs_tools.py b/src/ga4gh/vrs/utils/hgvs_tools.py index 503d4022..81ee9883 100644 --- a/src/ga4gh/vrs/utils/hgvs_tools.py +++ b/src/ga4gh/vrs/utils/hgvs_tools.py @@ -8,6 +8,7 @@ import hgvs.normalizer import hgvs.parser import hgvs.variantmapper +from bioutils.sequences import aa1_to_aa3_lut from hgvs.sequencevariant import SequenceVariant as HgvsSequenceVariant from ga4gh.vrs import models @@ -266,10 +267,7 @@ def _to_sequence_variant( """Create a SequenceVariant object from an Allele object.""" # build interval and edit depending on sequence type if sequence_type == "p": - msg = "Only nucleic acid variation is currently supported" - raise ValueError(msg) - # ival = hgvs.location.Interval(start=start, end=end) - # edit = hgvs.edit.AARefAlt(ref=None, alt=vo.state.sequence) + return self._to_protein_sequence_variant(vo, sequence, accession) start, end = vo.location.start, vo.location.end # ib: 0 1 2 3 4 5 # h: 1 2 3 4 5 @@ -310,6 +308,65 @@ def _to_sequence_variant( return var + def _to_protein_sequence_variant( + self, vo: models.Allele, sequence: str, accession: str + ) -> HgvsSequenceVariant: + """Create a protein SequenceVariant object from an Allele object. + + Only single-residue changes with a LiteralSequenceExpression state are + currently supported: missense/nonsense substitutions (e.g. + ``NP_060204.1:p.Val261Ala``) and single-residue deletions (e.g. + ``NP_060204.1:p.Val261del``). + + :param vo: VRS Allele object on a protein sequence + :param sequence: ga4gh sequence identifier for the reference sequence + :param accession: sequence accession for the HGVS expression + :return: HGVS SequenceVariant object + :raises ValueError: if the allele is not a single-residue change with a + LiteralSequenceExpression state, or the state cannot be expressed + as an HGVS protein variant + """ + start, end = vo.location.start, vo.location.end + if end - start != 1: + msg = "Only single-residue protein changes are currently supported" + raise ValueError(msg) + if vo.state.type != models.VrsType.LIT_SEQ_EXPR.value: + msg = "Only LiteralSequenceExpression states are currently supported for protein variants" + raise ValueError(msg) + + ref = self.data_proxy.get_sequence(sequence, start, end) + alt = str(vo.state.sequence.root) or None + if alt == ref: + msg = "Reference alleles cannot be expressed as HGVS protein variants" + raise ValueError(msg) + if alt is not None and ( + len(alt) != 1 or (alt not in aa1_to_aa3_lut and alt != "?") + ): + msg = f"Unsupported protein alternate residue: {alt!r}" + raise ValueError(msg) + + # protein coordinates are 1-based; the reference residue is carried + # by the position, per HGVS protein notation (e.g. p.Val261Ala) + pos = hgvs.location.AAPosition(base=start + 1, aa=ref) + ival = hgvs.location.Interval(start=pos, end=pos) + if alt is None: + edit = hgvs.edit.AARefAlt(ref="", alt=None) + else: + edit = hgvs.edit.AASub(ref=ref, alt=alt) + + posedit = hgvs.posedit.PosEdit(pos=ival, edit=edit) + var = HgvsSequenceVariant(ac=accession, type="p", posedit=posedit) + + try: + # if the namespace is GRC, can't normalize, since hgvs can't deal with it + parsed = self.parse(str(var)) + var = self.normalize(parsed) + + except hgvs.exceptions.HGVSDataNotAvailableError: + _logger.warning("No data found for accession %s", accession) + + return var + def normalize(self, hgvs: HgvsSequenceVariant) -> HgvsSequenceVariant: """Perform hgvs library `normalize()` method""" return self.normalizer.normalize(hgvs) diff --git a/tests/extras/test_to_hgvs_protein.py b/tests/extras/test_to_hgvs_protein.py new file mode 100644 index 00000000..d69c4d17 --- /dev/null +++ b/tests/extras/test_to_hgvs_protein.py @@ -0,0 +1,189 @@ +"""Tests for to_hgvs support of protein variants (single-residue changes). + +See https://github.com/ga4gh/vrs-python/issues/633. + +These tests use a stub data proxy and a mocked UTA connection so they run +without SeqRepo data or external services. +""" + +from unittest import mock + +import hgvs.exceptions +import pytest + +from ga4gh.vrs import models +from ga4gh.vrs.dataproxy import _DataProxy +from ga4gh.vrs.utils.hgvs_tools import HgvsTools + + +class _StubProteinDataProxy(_DataProxy): + """Minimal data proxy stub exposing a single protein sequence. + + Residue 261 (1-based) of refseq:NP_060204.1 is Val. + """ + + REFGET_ACCESSION = "SQ.dvlWjX2CGulwfb2ehmkCFn02ah7tEEVB" + ACCESSION = "NP_060204.1" + + def get_sequence( + self, identifier: str, start: int | None = None, end: int | None = None + ) -> str: + assert identifier == f"ga4gh:{self.REFGET_ACCESSION}" + assert (start, end) == (260, 261) + return "V" + + def get_metadata(self, identifier: str) -> dict: + assert identifier == f"ga4gh:{self.REFGET_ACCESSION}" + return { + "aliases": [ + f"ga4gh:{self.REFGET_ACCESSION}", + f"refseq:{self.ACCESSION}", + ], + "length": 597, + "alphabet": "ACDEFGHIKLMNPQRSTVWXY", + } + + +@pytest.fixture(scope="module") +def hgvs_tools(): + """HgvsTools with stubbed data proxy and no UTA connection. + + ``normalize`` is stubbed to raise HGVSDataNotAvailableError so the + no-UTA-data fallback path is exercised deterministically. + """ + with ( + mock.patch("hgvs.dataproviders.uta.connect", return_value=None), + mock.patch.object( + HgvsTools, + "normalize", + side_effect=hgvs.exceptions.HGVSDataNotAvailableError("no UTA data"), + ), + ): + yield HgvsTools(_StubProteinDataProxy()) + + +def _protein_allele( + start: int = 260, + end: int = 261, + state: dict | None = None, +) -> models.Allele: + return models.Allele.model_validate( + { + "type": "Allele", + "location": { + "type": "SequenceLocation", + "sequenceReference": { + "type": "SequenceReference", + "refgetAccession": _StubProteinDataProxy.REFGET_ACCESSION, + }, + "start": start, + "end": end, + }, + "state": state or {"type": "LiteralSequenceExpression", "sequence": "A"}, + } + ) + + +def test_to_hgvs_protein_substitution(hgvs_tools): + """The exact example from issue #633.""" + allele = models.Allele.model_validate( + { + "id": "ga4gh:VA.AIm-GH_iqp_bpcIVzi431fl7z-cQimNN", + "type": "Allele", + "digest": "AIm-GH_iqp_bpcIVzi431fl7z-cQimNN", + "location": { + "id": "ga4gh:SL.wgFs8Z2Nk4uD7nq_gFbdYH1FhagXkMhL", + "type": "SequenceLocation", + "digest": "wgFs8Z2Nk4uD7nq_gFbdYH1FhagXkMhL", + "sequenceReference": { + "type": "SequenceReference", + "refgetAccession": _StubProteinDataProxy.REFGET_ACCESSION, + }, + "start": 260, + "end": 261, + }, + "state": {"type": "LiteralSequenceExpression", "sequence": "A"}, + } + ) + assert hgvs_tools.from_allele(allele, namespace="refseq") == [ + "NP_060204.1:p.Val261Ala" + ] + + +def test_to_hgvs_protein_nonsense(hgvs_tools): + """Single-residue nonsense substitution -> p.Val261Ter.""" + allele = _protein_allele( + state={"type": "LiteralSequenceExpression", "sequence": "*"} + ) + assert hgvs_tools.from_allele(allele, namespace="refseq") == [ + "NP_060204.1:p.Val261Ter" + ] + + +def test_to_hgvs_protein_deletion(hgvs_tools): + """Single-residue deletion (empty alt) -> p.Val261del.""" + allele = _protein_allele( + state={"type": "LiteralSequenceExpression", "sequence": ""} + ) + assert hgvs_tools.from_allele(allele, namespace="refseq") == [ + "NP_060204.1:p.Val261del" + ] + + +def test_to_hgvs_protein_multi_residue_rejected(hgvs_tools): + """Alleles spanning more than one residue are not supported (yet).""" + allele = _protein_allele(start=260, end=262) + with pytest.raises(ValueError, match="single-residue"): + hgvs_tools.from_allele(allele, namespace="refseq") + + +def test_to_hgvs_protein_insertion_rejected(hgvs_tools): + """Insertions (start == end) are not single-residue changes.""" + allele = _protein_allele(start=260, end=260) + with pytest.raises(ValueError, match="single-residue"): + hgvs_tools.from_allele(allele, namespace="refseq") + + +def test_to_hgvs_protein_rle_state_rejected(hgvs_tools): + """ReferenceLengthExpression states are not supported for proteins.""" + allele = _protein_allele( + state={ + "type": "ReferenceLengthExpression", + "length": 1, + "repeatSubunitLength": 1, + "sequence": "V", + } + ) + with pytest.raises(ValueError, match="LiteralSequenceExpression"): + hgvs_tools.from_allele(allele, namespace="refseq") + + +def test_to_hgvs_protein_reference_allele_rejected(hgvs_tools): + """An allele identical to reference cannot be expressed as a variant.""" + allele = _protein_allele( + state={"type": "LiteralSequenceExpression", "sequence": "V"} + ) + with pytest.raises(ValueError, match="Reference alleles"): + hgvs_tools.from_allele(allele, namespace="refseq") + + +@pytest.mark.parametrize("alt", ["AA"]) +def test_to_hgvs_protein_multi_residue_alt_rejected(hgvs_tools, alt): + """Alternate alleles longer than one residue are rejected.""" + allele = _protein_allele( + state={"type": "LiteralSequenceExpression", "sequence": alt} + ) + with pytest.raises(ValueError, match="Unsupported protein"): + hgvs_tools.from_allele(allele, namespace="refseq") + + +def test_to_hgvs_protein_invalid_alt_rejected(hgvs_tools): + """Alternate residues that are not valid amino acid codes are rejected. + + (Constructed via attribute assignment since the VRS model itself + rejects such sequences at validation time.) + """ + allele = _protein_allele() + allele.state.sequence.root = "!" + with pytest.raises(ValueError, match="Unsupported protein"): + hgvs_tools.from_allele(allele, namespace="refseq")