feat: refactor hgvsTools and support c. translation

setup.cfg

@@ -92,7 +92,7 @@ extras =
     psycopg2-binary
     biocommons.seqrepo>=0.5.1
     bioutils>=0.5.2
-    hgvs@git+https://github.com/biocommons/hgvs@225-uncertain-ranges
+    hgvs>=1.4
     requests
     dill~=0.3.7
     click

src/ga4gh/vrs/_internal/models.py

@@ -414,7 +414,15 @@ class SequenceLocation(_Ga4ghIdentifiableObject):
     end: Union[Range, int] = Field(
         ...,
         description='The end coordinate or range of the SequenceLocation. The minimum value of this coordinate or range is 0. MUST represent a coordinate or range greater than the value of `start`.',
+
     )
+    def get_refget_accession(self):
+        if isinstance(self.sequenceReference, SequenceReference):
+            return self.sequenceReference.refgetAccession
+        elif isinstance(self.sequenceReference, IRI):
+            return self.sequenceReference.root
+        else:
+            return None
 
     class ga4gh(_Ga4ghIdentifiableObject.ga4gh):
         prefix = 'SL'
@@ -426,6 +434,7 @@ class ga4gh(_Ga4ghIdentifiableObject.ga4gh):
         ]
 
 
+
 class _VariationBase(_Ga4ghIdentifiableObject):
     """Base class for variation"""
 
@@ -526,27 +535,6 @@ class ga4gh(_Ga4ghIdentifiableObject.ga4gh):
         ]
 
 
-# class GenotypeMember(_ValueObject):
-#     """A class for expressing the count of a specific `MolecularVariation` present
-#     in-trans at a genomic locus represented by a `Genotype`.
-#     """
-#
-#     type: Literal['GenotypeMember'] = Field('GenotypeMember', description='MUST be "GenotypeMember".')
-#     count: Union[Range, int] = Field(
-#         ..., description='The number of copies of the `variation` at a Genotype locus.'
-#     )
-#     variation: Union[Allele, Haplotype] = Field(
-#         ..., description='A MolecularVariation at a Genotype locus.'
-#     )
-#
-#     class ga4gh(_Ga4ghIdentifiableObject.ga4gh):
-#         keys = [
-#             'type',
-#             'count',
-#             'variation'
-#         ]
-
-
 class MolecularVariation(RootModel):
     """A variation on a contiguous molecule."""
 
@@ -558,34 +546,6 @@ class MolecularVariation(RootModel):
         discriminator='type'
     )
 
-
-# class Genotype(_VariationBase):
-#     """A quantified set of _in-trans_ `MolecularVariation` at a genomic locus."""
-#
-#     type: Literal['Genotype'] = Field(
-#         'Genotype',
-#         description='MUST be "Genotype"'
-#     )
-#     # TODO members temporarily typed as List instead of Set + validate unique items
-#     members: List[GenotypeMember] = Field(
-#         ...,
-#         description='Each GenotypeMember in `members` describes a MolecularVariation and the count of that variation at the locus.',
-#         min_length=1,
-#     )
-#     count: Union[Range, int] = Field(
-#         ...,
-#         description='The total number of copies of all MolecularVariation at this locus, MUST be greater than or equal to the sum of GenotypeMember copy counts. If greater than the total counts, this implies additional MolecularVariation that are expected to exist but are not explicitly indicated.',
-#     )
-#
-#     class ga4gh(_Ga4ghIdentifiableObject.ga4gh):
-#         prefix = 'GT'
-#         keys = [
-#             'count',
-#             'members',
-#             'type'
-#         ]
-
-
 class SequenceExpression(RootModel):
     root: Union[LiteralSequenceExpression, ReferenceLengthExpression] = Field(
         ...,

src/ga4gh/vrs/dataproxy.py

@@ -6,6 +6,7 @@
 """
 
 from abc import ABC, abstractmethod
+from typing import Tuple
 from collections.abc import Sequence
 import datetime
 import functools
@@ -65,6 +66,36 @@ def get_metadata(self, identifier):
 
         """
 
+    @staticmethod
+    def extract_sequence_type(alias: str) -> str:
+        """
+        The purpose of this function is to provide a convenient way to extract the sequence type from an accession by matching its prefix to a known set of prefixes.
+
+        Args:
+        alias (str): The accession string.
+
+        Returns:
+        str or None: The sequence type associated with the accession string, or None if no matching prefix is found.
+        """
+
+        prefix_dict = {
+            "refseq:NM_": "c",
+            "refseq:NC_012920": "m",
+            "refseq:NG_": "g",
+            "refseq:NC_00": "g",
+            "refseq:NW_": "g",
+            "refseq:NR_": "n",
+            "refseq:NP_": "p",
+            "refseq:XM_": "c",
+            "refseq:XR_": "n",
+            "refseq:XP_": "p",
+            "GRCh": "g",
+        }
+
+        for prefix, seq_type in prefix_dict.items():
+            if alias.startswith(prefix):
+                return seq_type
+        return None
 
     @functools.lru_cache()
     def translate_sequence_identifier(self, identifier, namespace=None):
@@ -88,6 +119,56 @@ def translate_sequence_identifier(self, identifier, namespace=None):
             nsd = namespace + ":"
             aliases = [a for a in aliases if a.startswith(nsd)]
         return aliases
+
+    def derive_refget_accession(self, ac: str):
+        """Derive the refget accession from a public accession identifier
+
+        :param ac: public accession in simple or curie form from which to derive the refget accession
+        :return: Refget Accession if found
+        """
+
+        if ac is None:
+            return None
+
+        if ":" not in ac[1:]:
+            # always coerce the namespace if none provided 
+            ac = coerce_namespace(ac)
+
+        refget_accession = None
+        try:
+            aliases = self.translate_sequence_identifier(ac, namespace="ga4gh")
+        except KeyError:
+            _logger.error("KeyError when getting refget accession: %s", ac)
+        else:
+            if aliases:
+                refget_accession = aliases[0].split("ga4gh:")[-1]
+
+        return refget_accession
+
+    def is_valid_ref_seq(
+        self, sequence_id: str, start_pos: int, end_pos: int, ref: str
+    ) -> Tuple[bool, str]:
+        """Return wether or not the expected reference sequence matches the actual
+        reference sequence
+
+        :param sequence_id: Sequence ID to use
+        :param start_pos: Start pos (inter-residue) on the sequence_id
+        :param end_pos: End pos (inter-residue) on the sequence_id
+        :param ref: The expected reference sequence on the sequence_id given the
+            start_pos and end_pos
+        :return: Tuple containing whether or not actual reference sequence matches
+            the expected reference sequence and error message if mismatch
+        """
+        actual_ref = self.get_sequence(sequence_id, start_pos, end_pos)
+        is_valid = actual_ref == ref
+        err_msg = ""
+        if not is_valid:
+            err_msg = (
+                f"Expected reference sequence {ref} on {sequence_id} at positions "
+                f"({start_pos}, {end_pos}) but found {actual_ref}"
+            )
+            _logger.warning(err_msg)
+        return is_valid, err_msg
 
 
 class _SeqRepoDataProxyBase(_DataProxy):
@@ -224,16 +305,13 @@ def _isoformat(o):
     # eg: '2015-09-25T23:14:42.588601Z'
     return o.isoformat("T") + "Z"
 
-
 # Future implementations
 # * The RefGetDataProxy is waiting on support for sequence lookup by alias
 # class RefGetDataProxy(_DataProxy):
 #     def __init__(self, base_url):
 #         super().__init__()
 #         self.base_url = base_url
 
-
-
 def create_dataproxy(uri: str = None) -> _DataProxy:
     """Create a dataproxy from uri or GA4GH_VRS_DATAPROXY_URI