Fix allele translator to_spdi issue for RLE deletion. Add RLE reconstruction when longer than rle_seq_limit. (#542)

This fixes a bug for `_to_spdi` for ReferenceLengthExpressions. 

It also adds an optional argument to `translate_to` called
`ref_seq_limit` which works similarly to `rle_seq_limit`, but for SPDI
expressions. If `ref_seq_limit` is nonzero, the 3rd term of the SPDI
expression will be the actual reference sequence, not the length.

To round-trip a clinvar SPDI expression which contains the reference
sequence:

```
spdi = "NC_000019.10:44908821:C:T"
vrs = translator.translate_from(spdi, "spdi")
to_spdi = translator.translate_to(vrs, "spdi", ref_seq_limit=None) # or some number >=1 since len("C")==1
assert spdi == to_spdi
```

To continue the existing behavior of using the reference sequence
length:
```
spdi = "NC_000019.10:44908821:1:T"
vrs = translator.translate_from(spdi, "spdi")
to_spdi = translator.translate_to(vrs, "spdi") # ref_seq_limit defaults to 0
assert spdi == to_spdi
```

---------

Co-authored-by: Kori Kuzma <korikuzma@gmail.com>

Author: theferrit32

src/ga4gh/vrs/extras/translator.py

@@ -9,16 +9,40 @@
 import re
 from abc import ABC
 from collections.abc import Mapping
+from typing import Protocol
 
 from ga4gh.core import ga4gh_identify
 from ga4gh.vrs import models, normalize
-from ga4gh.vrs.dataproxy import _DataProxy
+from ga4gh.vrs.dataproxy import SequenceProxy, _DataProxy
 from ga4gh.vrs.extras.decorators import lazy_property
+from ga4gh.vrs.normalize import denormalize_reference_length_expression
 from ga4gh.vrs.utils.hgvs_tools import HgvsTools
 
 _logger = logging.getLogger(__name__)
 
 
+class VariationToStrProtocol(Protocol):
+    """Protocol for translating VRS objects to other string expressions.
+
+    This protocol defines a callable interface for translating a VRS object
+    into variation strings, with optional keyword arguments for customization.
+    """
+
+    def __call__(self, vo: models._VariationBase, **kwargs) -> list[str]:
+        """Translate vrs object `vo` to variation string expressions"""
+
+
+class VariationFromStrProtocol(Protocol):
+    """Protocol for translating variation strings to VRS objects.
+
+    This protocol defines a callable interface for translating a variation
+    string into a VRS object, with optional keyword arguments for customization.
+    """
+
+    def __call__(self, expr: str, **kwargs) -> models._VariationBase | None:
+        """Translate variation string `expr` to a VRS object"""
+
+
 class _Translator(ABC):  # noqa: B024
     """abstract class / interface for VRS to/from translation needs
 
@@ -54,8 +78,8 @@ def __init__(
         self.data_proxy = data_proxy
         self.identify = identify
         self.rle_seq_limit = rle_seq_limit
-        self.from_translators = {}
-        self.to_translators = {}
+        self.from_translators: dict[str, VariationFromStrProtocol] = {}
+        self.to_translators: dict[str, VariationToStrProtocol] = {}
 
     def translate_from(
         self, var: str, fmt: str | None = None, **kwargs
@@ -110,10 +134,15 @@ def translate_from(
         msg = f"Unable to parse data as {', '.join(formats)}"
         raise ValueError(msg)
 
-    def translate_to(self, vo: models._VariationBase, fmt: str) -> str:
-        """Translate vrs object `vo` to named format `fmt`"""
+    def translate_to(self, vo: models._VariationBase, fmt: str, **kwargs) -> list[str]:
+        """Translate vrs object `vo` to named format `fmt`
+
+        kwargs:
+            ref_seq_limit Optional(int):
+                If vo.state is a ReferenceLengthExpression, and `ref_seq_limit` is specified, and `fmt` is `spdi`, the reference sequence is included in the SPDI expression if it is below the limit Otherwise only the length of the reference sequence is included. If the limit is None, the reference sequence is always included. In all cases, the alt sequence is included. Default is 0 (never include reference sequence).
+        """
         t = self.to_translators[fmt]
-        return t(vo)
+        return t(vo, **kwargs)
 
     ############################################################################
     # INTERNAL
@@ -385,12 +414,15 @@ def _from_spdi(self, spdi_expr: str, **kwargs) -> models.Allele | None:
         return self._create_allele(values, **kwargs)
 
     def _to_hgvs(
-        self, vo: models.Allele, namespace: str | None = "refseq"
+        self,
+        vo: models.Allele,
+        namespace: str | None = "refseq",
+        **kwargs,  # noqa: ARG002
     ) -> list[str]:
         return self.hgvs_tools.from_allele(vo, namespace)
 
     def _to_spdi(
-        self, vo: models.Allele, namespace: str | None = "refseq"
+        self, vo: models.Allele, namespace: str | None = "refseq", **kwargs
     ) -> list[str]:
         """Generate a *list* of SPDI expressions for VRS Allele.
 
@@ -400,6 +432,13 @@ def _to_spdi(
         If `namespace` is None, returns SPDI strings for all alias
         translations.
 
+        If `ref_seq_limit` is specified, the reference sequence is
+        included in the SPDI expression only if it is below the limit.
+        Otherwise only the length of the reference sequence is
+        included. If the limit is None, the reference sequence is
+        always included. In all cases, the alt sequence is included.
+        Default is 0 (never include reference sequence).
+
         If no alias translations are available, an empty list is
         returned.
 
@@ -411,9 +450,43 @@ def _to_spdi(
         sequence = f"ga4gh:{vo.location.get_refget_accession()}"
         aliases = self.data_proxy.translate_sequence_identifier(sequence, namespace)
         aliases = [a.split(":")[1] for a in aliases]
+        seq_proxies = {a: SequenceProxy(self.data_proxy, a) for a in aliases}
         start, end = vo.location.start, vo.location.end
-        spdi_tail = f":{start}:{end - start}:{vo.state.sequence.root}"
-        return [a + spdi_tail for a in aliases]
+        spdi_exprs = []
+
+        for alias in aliases:
+            # Get the reference sequence
+            seq_proxy = seq_proxies[alias]
+            ref_seq = seq_proxy[start:end]
+
+            if vo.state.type == models.VrsType.REF_LEN_EXPR.value:
+                # Derived from reference. sequence included if under limit, but
+                # we can derive it again from the reference.
+                alt_seq = denormalize_reference_length_expression(
+                    ref_seq=ref_seq,
+                    repeat_subunit_length=vo.state.repeatSubunitLength,
+                    alt_length=vo.state.length,
+                )
+                # Warn if the derived sequence is different from the one in the object
+                if vo.state.sequence and vo.state.sequence != alt_seq:
+                    _logger.warning(
+                        "Derived sequence '%s' is different from provided state.sequence '%s'",
+                        alt_seq,
+                        vo.state.sequence,
+                    )
+            else:
+                alt_seq = vo.state.sequence.root
+
+            # Optionally allow using the length of the reference sequence
+            # instead of the sequence itself.
+            ref_seq_limit = kwargs.get("ref_seq_limit", 0)
+            if ref_seq_limit is not None and len(ref_seq) > int(ref_seq_limit):
+                ref_seq = len(ref_seq)
+
+            spdi_expr = f"{alias}:{start}:{ref_seq}:{alt_seq}"
+            spdi_exprs.append(spdi_expr)
+
+        return spdi_exprs
 
     def _post_process_imported_allele(
         self, allele: models.Allele, **kwargs

src/ga4gh/vrs/normalize.py

@@ -221,6 +221,33 @@ def _normalize_allele(input_allele, data_proxy, rle_seq_limit=50):
     return new_allele
 
 
+def denormalize_reference_length_expression(
+    ref_seq: str,
+    repeat_subunit_length: int,
+    alt_length: int,
+) -> str:
+    """Reverse the process of repeat subunit compaction for a ReferenceLengthExpression.
+    The alt is reference-derived so repeat the repeat subunit until it is it to alt_length
+    characters long. May include a trailing partial copy of the repeat subunit.
+
+    e.g. "ACGT" (repeat_subunit_length=4, length=10) -> "ACGTACGTAC" (ACGT ACGT AC)
+
+    :param ref_seq: The reference sequence to be denormalized.
+    :param repeat_subunit_length: The length of the repeat subunit.
+    :param alt_length: The length of the alternate sequence that was compacted during normalization.
+    """
+    repeat_subunit = ref_seq[:repeat_subunit_length]
+    if len(repeat_subunit) != repeat_subunit_length:
+        raise ValueError(  # noqa: TRY003
+            f"Repeat subunit length {repeat_subunit_length} is not equal to the length of the repeat subunit {len(repeat_subunit)}"  # noqa: EM102
+        )
+    repeat_count = alt_length // repeat_subunit_length
+    remainder = alt_length % repeat_subunit_length
+    alt = repeat_subunit * repeat_count
+    alt += repeat_subunit[:remainder]
+    return alt
+
+
 def _factor_gen(n):
     """Yield all factors of an integer `n`, in descending order"""
     lower_factors = []
@@ -259,7 +286,10 @@ def _is_valid_cycle(template_start, template, target):
 # TODO _normalize_genotype?
 
 
-def _normalize_cis_phased_block(o, data_proxy: _DataProxy | None = None):  # noqa: ARG001
+def _normalize_cis_phased_block(
+    o,
+    data_proxy: _DataProxy | None = None,  # noqa: ARG001
+):
     o.members = sorted(o.members, key=ga4gh_digest)
     return o
 

tests/extras/cassettes/test_from_hgvs_cx[NC_000013.11:g.32379315_32379819del-None-expected1].yaml

@@ -0,0 +1,44 @@
+interactions:
+- request:
+    body: null
+    headers:
+      Accept:
+      - '*/*'
+      Accept-Encoding:
+      - gzip, deflate
+      Connection:
+      - keep-alive
+      User-Agent:
+      - python-requests/2.32.3
+    method: GET
+    uri: http://localhost:5001/seqrepo/1/metadata/refseq:NC_000013.11
+  response:
+    body:
+      string: "{\n  \"added\": \"2016-08-27T23:50:14Z\",\n  \"aliases\": [\n    \"GRCh38:13\",\n
+        \   \"GRCh38:chr13\",\n    \"GRCh38.p1:13\",\n    \"GRCh38.p1:chr13\",\n    \"GRCh38.p10:13\",\n
+        \   \"GRCh38.p10:chr13\",\n    \"GRCh38.p11:13\",\n    \"GRCh38.p11:chr13\",\n
+        \   \"GRCh38.p12:13\",\n    \"GRCh38.p12:chr13\",\n    \"GRCh38.p2:13\",\n
+        \   \"GRCh38.p2:chr13\",\n    \"GRCh38.p3:13\",\n    \"GRCh38.p3:chr13\",\n
+        \   \"GRCh38.p4:13\",\n    \"GRCh38.p4:chr13\",\n    \"GRCh38.p5:13\",\n    \"GRCh38.p5:chr13\",\n
+        \   \"GRCh38.p6:13\",\n    \"GRCh38.p6:chr13\",\n    \"GRCh38.p7:13\",\n    \"GRCh38.p7:chr13\",\n
+        \   \"GRCh38.p8:13\",\n    \"GRCh38.p8:chr13\",\n    \"GRCh38.p9:13\",\n    \"GRCh38.p9:chr13\",\n
+        \   \"MD5:a5437debe2ef9c9ef8f3ea2874ae1d82\",\n    \"NCBI:NC_000013.11\",\n
+        \   \"refseq:NC_000013.11\",\n    \"SEGUID:2oDBty0yKV9wHo7gg+Bt+fPgi5o\",\n
+        \   \"SHA1:da80c1b72d32295f701e8ee083e06df9f3e08b9a\",\n    \"VMC:GS__0wi-qoDrvram155UmcSC-zA5ZK4fpLT\",\n
+        \   \"sha512t24u:_0wi-qoDrvram155UmcSC-zA5ZK4fpLT\",\n    \"ga4gh:SQ._0wi-qoDrvram155UmcSC-zA5ZK4fpLT\"\n
+        \ ],\n  \"alphabet\": \"ACGKNTY\",\n  \"length\": 114364328\n}\n"
+    headers:
+      Connection:
+      - close
+      Content-Length:
+      - '1002'
+      Content-Type:
+      - application/json
+      Date:
+      - Wed, 09 Apr 2025 19:40:39 GMT
+      Server:
+      - Werkzeug/2.2.3 Python/3.10.12
+    status:
+      code: 200
+      message: OK
+version: 1

tests/extras/cassettes/test_hgvs[NC_000013.11:g.32936732=-expected0].yaml

@@ -23,7 +23,37 @@ interactions:
       Content-Type:
       - text/plain; charset=utf-8
       Date:
-      - Mon, 10 Mar 2025 16:22:49 GMT
+      - Wed, 16 Apr 2025 22:49:30 GMT
+      Server:
+      - Werkzeug/2.2.3 Python/3.10.12
+    status:
+      code: 200
+      message: OK
+- request:
+    body: null
+    headers:
+      Accept:
+      - '*/*'
+      Accept-Encoding:
+      - gzip, deflate
+      Connection:
+      - keep-alive
+      User-Agent:
+      - python-requests/2.32.3
+    method: GET
+    uri: http://localhost:5000/seqrepo/1/sequence/ga4gh:SQ._0wi-qoDrvram155UmcSC-zA5ZK4fpLT?start=32936731&end=32936732
+  response:
+    body:
+      string: C
+    headers:
+      Connection:
+      - close
+      Content-Length:
+      - '1'
+      Content-Type:
+      - text/plain; charset=utf-8
+      Date:
+      - Wed, 16 Apr 2025 22:49:30 GMT
       Server:
       - Werkzeug/2.2.3 Python/3.10.12
     status:
@@ -67,20 +97,20 @@ interactions:
       Content-Type:
       - text/plain
       Date:
-      - Mon, 10 Mar 2025 16:22:50 GMT
+      - Wed, 16 Apr 2025 22:49:30 GMT
       Keep-Alive:
       - timeout=4, max=40
       NCBI-PHID:
-      - D0BDD793D08B96450000396AA7F60D6E.1.1.m_5
+      - 322CF8979F5E85B5000050CD47341D44.1.1.m_5
       NCBI-SID:
-      - A10F065EB442C5FB_D344SID
+      - B4DD900AE926F06D_8D6DSID
       Referrer-Policy:
       - origin-when-cross-origin
       Server:
       - Finatra
       Set-Cookie:
-      - ncbi_sid=A10F065EB442C5FB_D344SID; domain=.nih.gov; path=/; expires=Tue, 10
-        Mar 2026 16:22:50 GMT
+      - ncbi_sid=B4DD900AE926F06D_8D6DSID; domain=.nih.gov; path=/; expires=Thu, 16
+        Apr 2026 22:49:30 GMT
       Strict-Transport-Security:
       - max-age=31536000; includeSubDomains; preload
       Transfer-Encoding: