From 6624c3d7c0fcf8a399fa63d876c8adc23c71c3d4 Mon Sep 17 00:00:00 2001 From: Rakesh Pai <41351936+developer-rpai@users.noreply.github.com> Date: Wed, 23 Sep 2026 22:15:37 -0700 Subject: [PATCH] fix: validate reference allele stated in HGVS expressions The HGVS parser preserves the input reference allele (e.g. the C in NM_006087.3:c.900C>A), but the translator never checked it against the reference sequence, so an incorrect reference allele silently produced a plausible-but-wrong VRS Allele (ga4gh/vrs-python#364). extract_allele_values now validates the stated reference allele against the data proxy (default require_validation=True), mirroring the existing gnomAD translator behavior. require_validation=False keeps the legacy behavior (mismatch only logged). Adds hermetic regression tests and the VCR cassette interactions the new lookups require. --- src/ga4gh/vrs/extras/translator.py | 19 ++- src/ga4gh/vrs/utils/hgvs_tools.py | 27 +++- tests/extras/cassettes/test_from_hgvs.yaml | 24 ++++ ...NC_000007.14:g.55181320A>T-expected1].yaml | 12 ++ ...vs[NM_001331029.1:c.722A>G-expected6].yaml | 12 ++ ...hgvs[NM_181798.1:c.1007G>T-expected7].yaml | 12 ++ tests/extras/test_allele_translator.py | 135 +++++++++++++++++- 7 files changed, 238 insertions(+), 3 deletions(-) diff --git a/src/ga4gh/vrs/extras/translator.py b/src/ga4gh/vrs/extras/translator.py index 91c45532..50cb4eee 100644 --- a/src/ga4gh/vrs/extras/translator.py +++ b/src/ga4gh/vrs/extras/translator.py @@ -357,7 +357,24 @@ def _from_gnomad(self, gnomad_expr: str, **kwargs) -> models.Allele | None: return self._create_allele(values, **kwargs) def _from_hgvs(self, hgvs_expr: str, **kwargs) -> models.Allele | None: - allele_values = self.hgvs_tools.extract_allele_values(hgvs_expr) + """Parse HGVS expression into VRS Allele + + kwargs: + require_validation (bool): If `True` then validation checks must pass in + order to return a VRS object. A `DataProxyValidationError` will be + raised if validation checks fail. If `False` then VRS object will be + returned even if validation checks fail. Defaults to `True`. + rle_seq_limit Optional(int): If RLE is set as the new state after + normalization, this sets the limit for the length of the `sequence`. + To exclude `sequence` from the response, set to 0. + For no limit, set to `None`. + Defaults value set in instance variable, `rle_seq_limit`. + do_normalize (bool): `True` if fully justified normalization should be + performed. `False` otherwise. Defaults to `True` + """ + allele_values = self.hgvs_tools.extract_allele_values( + hgvs_expr, require_validation=kwargs.get("require_validation", True) + ) if allele_values: return self._create_allele(allele_values, **kwargs) return None diff --git a/src/ga4gh/vrs/utils/hgvs_tools.py b/src/ga4gh/vrs/utils/hgvs_tools.py index 503d4022..09b872bb 100644 --- a/src/ga4gh/vrs/utils/hgvs_tools.py +++ b/src/ga4gh/vrs/utils/hgvs_tools.py @@ -131,9 +131,18 @@ def get_position_and_state(self, sv: HgvsSequenceVariant) -> tuple[int, int, str return start, end, state - def extract_allele_values(self, hgvs_expr: str) -> dict | None: + def extract_allele_values( + self, hgvs_expr: str, require_validation: bool = True + ) -> dict | None: """Parse hgvs into a VRS Allele Object + :param hgvs_expr: HGVS expression to parse + :param require_validation: If `True`, the reference allele stated in the + HGVS expression (when present, e.g. the `C` in `c.900C>A`) must match + the actual reference sequence, otherwise a + `DataProxyValidationError` is raised. If `False`, a mismatch is only + logged and the allele is still returned. Defaults to `True`. + kwargs: rle_seq_limit Optional(int): If RLE is set as the new state after normalization, this sets the limit for the length of the `sequence`. @@ -181,6 +190,22 @@ def extract_allele_values(self, hgvs_expr: str) -> dict | None: (start, end, state) = self.get_position_and_state(sv) + # The HGVS expression may state the expected reference allele (e.g. the + # `C` in `c.900C>A`, or the deleted bases in `g.44908822delC`). A + # mismatch means the input itself is invalid: silently emitting a VRS + # Allele would produce a plausible-but-wrong variant object + # (ga4gh/vrs-python#364), so validate it against the reference + # sequence just like the gnomAD translator does. + ref_allele = getattr(sv.posedit.edit, "ref", None) + if ref_allele: + self.data_proxy.validate_ref_seq( + sv.ac, + start, + end, + ref_allele.upper(), + require_validation=require_validation, + ) + return { "refget_accession": refget_accession, "start": start, diff --git a/tests/extras/cassettes/test_from_hgvs.yaml b/tests/extras/cassettes/test_from_hgvs.yaml index ed06582a..15da4b72 100644 --- a/tests/extras/cassettes/test_from_hgvs.yaml +++ b/tests/extras/cassettes/test_from_hgvs.yaml @@ -88,4 +88,28 @@ interactions: status: code: 200 message: OK +- request: + body: null + headers: {} + method: GET + uri: http://localhost:5000/seqrepo/1/sequence/NC_000019.10?start=44908821&end=44908822 + response: + body: + string: C + headers: {} + status: + code: 200 + message: OK +- request: + body: null + headers: {} + method: GET + uri: http://localhost:5000/seqrepo/1/sequence/NC_012920.1?start=10082&end=10083 + response: + body: + string: A + headers: {} + status: + code: 200 + message: OK version: 1 diff --git a/tests/extras/cassettes/test_hgvs[NC_000007.14:g.55181320A>T-expected1].yaml b/tests/extras/cassettes/test_hgvs[NC_000007.14:g.55181320A>T-expected1].yaml index 3c14c730..dc0c584f 100644 --- a/tests/extras/cassettes/test_hgvs[NC_000007.14:g.55181320A>T-expected1].yaml +++ b/tests/extras/cassettes/test_hgvs[NC_000007.14:g.55181320A>T-expected1].yaml @@ -70,4 +70,16 @@ interactions: status: code: 200 message: OK +- request: + body: null + headers: {} + method: GET + uri: http://localhost:5000/seqrepo/1/sequence/NC_000007.14?start=55181319&end=55181320 + response: + body: + string: A + headers: {} + status: + code: 200 + message: OK version: 1 diff --git a/tests/extras/cassettes/test_hgvs[NM_001331029.1:c.722A>G-expected6].yaml b/tests/extras/cassettes/test_hgvs[NM_001331029.1:c.722A>G-expected6].yaml index 44352fa1..0b5032f5 100644 --- a/tests/extras/cassettes/test_hgvs[NM_001331029.1:c.722A>G-expected6].yaml +++ b/tests/extras/cassettes/test_hgvs[NM_001331029.1:c.722A>G-expected6].yaml @@ -81,4 +81,16 @@ interactions: status: code: 200 message: OK +- request: + body: null + headers: {} + method: GET + uri: http://localhost:5000/seqrepo/1/sequence/NM_001331029.1?start=871&end=872 + response: + body: + string: A + headers: {} + status: + code: 200 + message: OK version: 1 diff --git a/tests/extras/cassettes/test_hgvs[NM_181798.1:c.1007G>T-expected7].yaml b/tests/extras/cassettes/test_hgvs[NM_181798.1:c.1007G>T-expected7].yaml index 20c3a261..b53c0b7a 100644 --- a/tests/extras/cassettes/test_hgvs[NM_181798.1:c.1007G>T-expected7].yaml +++ b/tests/extras/cassettes/test_hgvs[NM_181798.1:c.1007G>T-expected7].yaml @@ -99,4 +99,16 @@ interactions: status: code: 200 message: OK +- request: + body: null + headers: {} + method: GET + uri: http://localhost:5000/seqrepo/1/sequence/NM_181798.1?start=1262&end=1263 + response: + body: + string: G + headers: {} + status: + code: 200 + message: OK version: 1 diff --git a/tests/extras/test_allele_translator.py b/tests/extras/test_allele_translator.py index 15f87a89..79b6a80d 100644 --- a/tests/extras/test_allele_translator.py +++ b/tests/extras/test_allele_translator.py @@ -1,8 +1,12 @@ +from typing import ClassVar +from unittest.mock import MagicMock, patch + import pytest from ga4gh.vrs import models -from ga4gh.vrs.dataproxy import DataProxyValidationError +from ga4gh.vrs.dataproxy import DataProxyValidationError, _DataProxy from ga4gh.vrs.extras.translator import AlleleTranslator +from ga4gh.vrs.utils.hgvs_tools import HgvsTools @pytest.fixture(scope="module") @@ -982,3 +986,132 @@ def test_normalize_microsatellite_counts(tlr, case): def test_translate_to_invalid_fmt(tlr): with pytest.raises(NotImplementedError, match="gnomad is not supported"): tlr.translate_to(models.Allele.model_validate(snv_output), fmt="gnomad") + + +# --------------------------------------------------------------------------- +# Regression tests for https://github.com/ga4gh/vrs-python/issues/364 +# ("hgvs to vrs is returning valid results when hgvs has IncorrectReferenceAllele") +# +# These tests are hermetic: reference-sequence lookups are served from canned +# ground truth, so they run without seqrepo/UTA network access. + + +class _CannedDataProxy(_DataProxy): + """Minimal data proxy backed by canned ground-truth reference sequences. + + Only the sequence/metadata lookups are faked; reference validation uses + the real `_DataProxy.validate_ref_seq` logic. + """ + + # (accession, interbase start, interbase end) -> true reference sequence + TRUTH: ClassVar[dict] = { + # NM_006087.3 (TUBB4A): CDS is 373..1707, so c.900 == n.1272. The true + # base there is G, as reported by the ClinGen Allele Registry for the + # NM_006087.3:c.900C>A expression in issue #364 (IncorrectReferenceAllele: + # "given=C, found=G"), independently confirmed against NCBI RefSeq. + ("NM_006087.3", 1271, 1272): "G", + # GRCh38 chr19:44908822, true ref C (matches the C>T test expression) + ("NC_000019.10", 44908821, 44908822): "C", + } + + def get_sequence( + self, identifier: str, start: int | None = None, end: int | None = None + ) -> str: + return self.TRUTH[(identifier, start, end)] + + def get_metadata(self, _identifier: str) -> dict: + return {"aliases": ["ga4gh:SQ." + "A" * 32], "length": 10**6} + + +def _c_to_n_nm006087(_self, sv): + """Emulate the UTA c.->n. mapping for NM_006087.3. + + UTA is not reachable from every test environment, so apply the true + mapping directly: the RefSeq CDS annotation for NM_006087.3 is 373..1707, + hence c.900 maps to n.1272. + """ + assert sv.ac == "NM_006087.3" + offset = 372 # n. coordinate == c. coordinate + 372 on this transcript + sv.posedit.pos.start.base += offset + sv.posedit.pos.end.base += offset + sv.type = "n" + return sv + + +@pytest.fixture +def tlr_canned(): + """AlleleTranslator backed by canned reference data (no network/UTA).""" + with ( + patch("hgvs.dataproviders.uta.connect", return_value=MagicMock()), + patch.object(HgvsTools, "c_to_n", _c_to_n_nm006087), + ): + yield AlleleTranslator(data_proxy=_CannedDataProxy(), identify=False) + + +def test_from_hgvs_wrong_ref_allele_raises(tlr_canned): + """A mismatched reference allele in the HGVS expression must raise, not + silently produce a plausible-but-wrong VRS Allele (issue #364). + """ + error_msg = ( + "Reference mismatch at NM_006087.3 position 1271-1272 " + "(input gave 'C' but correct ref is 'G')" + ) + + with pytest.raises(DataProxyValidationError) as e: + tlr_canned._from_hgvs("NM_006087.3:c.900C>A", do_normalize=False) + assert str(e.value) == error_msg + + with pytest.raises(DataProxyValidationError) as e: + tlr_canned.translate_from( + "NM_006087.3:c.900C>A", fmt="hgvs", do_normalize=False + ) + assert str(e.value) == error_msg + + +def test_from_hgvs_wrong_ref_allele_no_validation(tlr_canned): + """require_validation=False keeps the legacy behavior: the allele is + returned and the mismatch is only logged. + """ + allele = tlr_canned._from_hgvs( + "NM_006087.3:c.900C>A", do_normalize=False, require_validation=False + ) + assert (allele.location.start, allele.location.end) == (1271, 1272) + assert allele.state.sequence.root == "A" + + +def test_from_hgvs_correct_ref_allele_passes(tlr_canned): + """HGVS expressions whose stated ref matches the reference sequence still + translate cleanly (substitution and deletion-with-ref forms). + """ + allele = tlr_canned.translate_from( + "NC_000019.10:g.44908822C>T", fmt="hgvs", do_normalize=False + ) + assert (allele.location.start, allele.location.end) == (44908821, 44908822) + assert allele.state.sequence.root == "T" + + allele = tlr_canned.translate_from( + "NC_000019.10:g.44908822delC", fmt="hgvs", do_normalize=False + ) + assert (allele.location.start, allele.location.end) == (44908821, 44908822) + assert allele.state.sequence.root == "" + + +def test_from_hgvs_wrong_ref_allele_del_raises(tlr_canned): + """The deletion-with-ref form is validated too.""" + with pytest.raises(DataProxyValidationError) as e: + tlr_canned.translate_from( + "NC_000019.10:g.44908822delG", fmt="hgvs", do_normalize=False + ) + assert "correct ref is 'C'" in str(e.value) + + +def test_from_hgvs_no_ref_allele_skips_validation(tlr_canned): + """Edits that state no reference allele (ins/dup/bare del) skip validation + entirely -- the canned proxy raises KeyError on any sequence lookup, so + this fails if validation is attempted. + """ + allele = tlr_canned.translate_from( + "NC_000019.10:g.44908822_44908823insT", fmt="hgvs", do_normalize=False + ) + assert (allele.location.start, allele.location.end) == (44908822, 44908822) + assert allele.state.sequence.root == "T"