stash

Author: jsstevenson

src/cool_seq_tool/mappers/exon_genomic_coords.py

@@ -944,10 +944,10 @@ async def _genomic_to_tx_segment(
                         )
         # gene is not required to liftover coordinates if tx_ac and genomic_ac are given, but we should set the associated gene
         if not gene:
-            _gene, err_msg = await self._get_tx_ac_gene(transcript)
-            if err_msg:
-                return GenomicTxSeg(errors=[err_msg])
-            gene = _gene
+            async with self.uta_db.repository() as uta:
+                gene = await uta.get_gene_from_tx_ac(transcript)
+            if not gene:
+                return GenomicTxSeg(errors=[f"No gene(s) found given {transcript}"])
 
         tx_exons = await self._get_all_exon_coords(
             tx_ac=transcript, genomic_ac=genomic_ac
@@ -967,9 +967,10 @@ async def _genomic_to_tx_segment(
 
         # Validate that the breakpoint occurs within 150 bp of the first and last exon for the selected transcript.
         # A breakpoint beyond this range is likely erroneous.
-        coordinate_check = await self._validate_genomic_breakpoint(
-            pos=genomic_pos, genomic_ac=genomic_ac, tx_ac=transcript
-        )
+        async with self.uta_db.repository() as uta:
+            coordinate_check = await uta.validate_genomic_breakpoint(
+                pos=genomic_pos, genomic_ac=genomic_ac, tx_ac=transcript
+            )
         if not coordinate_check:
             msg = (
                 f"{genomic_pos} on {genomic_ac} occurs more than 150 bp outside the "
@@ -1053,69 +1054,6 @@ async def _get_grch38_pos(
         )
         return liftover_data[1] if liftover_data else None
 
-    async def _validate_genomic_breakpoint(
-        self,
-        pos: int,
-        genomic_ac: str,
-        tx_ac: str,
-    ) -> bool:
-        """Validate that a genomic coordinate falls within the first and last exon
-            for a transcript on a given accession
-
-        :param pos: Genomic position on ``genomic_ac``
-        :param genomic_ac: RefSeq genomic accession, e.g. ``"NC_000007.14"``
-        :param transcript: A transcript accession
-        :return: ``True`` if the coordinate falls within 150bp of the first and last exon
-            for the transcript, ``False`` if not. Breakpoints past this threshold
-            are likely erroneous.
-        """
-        query = """
-            WITH tx_boundaries AS (
-                SELECT
-                MIN(alt_start_i) AS min_start,
-                MAX(alt_end_i) AS max_end
-                FROM tx_exon_aln_mv
-                WHERE tx_ac = %(tx_ac)s
-                AND alt_ac = %(genomic_ac)s
-            )
-            SELECT * FROM tx_boundaries
-            WHERE %(pos)s between (tx_boundaries.min_start - 150) and (tx_boundaries.max_end + 150);
-            """
-        async with self.uta_db.repository() as uta:
-            cursor = await uta.execute_query(
-                query, {"tx_ac": tx_ac, "genomic_ac": genomic_ac, "pos": pos}
-            )
-            result = await cursor.fetchone()
-        return bool(result)
-
-    async def _get_tx_ac_gene(
-        self,
-        tx_ac: str,
-    ) -> tuple[str | None, str | None]:
-        """Get gene given a transcript.
-
-        If multiple genes are found for a given ``tx_ac``, only one
-        gene will be returned.
-
-        :param tx_ac: RefSeq transcript, e.g. ``"NM_004333.6"``
-        :return: HGNC gene symbol associated to transcript and
-            warning
-        """
-        query = """
-            SELECT DISTINCT hgnc
-            FROM tx_exon_aln_mv
-            WHERE tx_ac = %(tx_ac)s
-            ORDER BY hgnc
-            LIMIT 1;
-            """
-        async with self.uta_db.repository() as uta:
-            cursor = await uta.execute_query(query, {"tx_ac": tx_ac})
-            result = await cursor.fetchone()
-        if not result:
-            return None, f"No gene(s) found given {tx_ac}"
-
-        return result[0], None
-
     @staticmethod
     def _is_exonic_breakpoint(pos: int, tx_genomic_coords: list[_ExonCoord]) -> bool:
         """Check if a breakpoint occurs on an exon

src/cool_seq_tool/mappers/mane_transcript.py

@@ -236,15 +236,7 @@ async def _liftover_to_38(self, genomic_tx_data: GenomicTxMetadata) -> None:
         chromosome, _ = descr
 
         async with self.uta_db.repository() as uta:
-            cursor = await uta.execute_query(
-                """SELECT DISTINCT alt_ac
-                FROM tx_exon_aln_mv
-                WHERE tx_ac = %(tx_ac)s;
-                """,
-                {"tx_ac": genomic_tx_data.tx_ac},
-            )
-            nc_acs = await cursor.fetchall()
-        nc_acs = [nc_ac[0] for nc_ac in nc_acs]
+            nc_acs = await uta.get_alt_acs_for_tx(genomic_tx_data.tx_ac)
         if nc_acs == [genomic_tx_data.alt_ac]:
             _logger.warning(
                 "UTA does not have GRCh38 assembly for %s",
@@ -263,26 +255,9 @@ async def _liftover_to_38(self, genomic_tx_data: GenomicTxMetadata) -> None:
             genomic_tx_data, "alt_pos_change_range", chromosome, Assembly.GRCH38
         )
 
-        # Change alt_ac to most recent
-        if genomic_tx_data.alt_ac.startswith("EN"):
-            order_by_cond = "ORDER BY alt_ac DESC;"
-        else:
-            order_by_cond = """
-            ORDER BY CAST(SUBSTR(alt_ac, position('.' in alt_ac) + 1,
-            LENGTH(alt_ac)) AS INT) DESC;
-            """
-        query = f"""
-            SELECT alt_ac
-            FROM genomic
-            WHERE alt_ac LIKE %(ac_pattern)s
-            {order_by_cond}
-            """  # noqa: S608
         async with self.uta_db.repository() as uta:
-            cursor = await uta.execute_query(
-                query, {"ac_pattern": f"{genomic_tx_data.alt_ac.split('.')[0]}%"}
-            )
-            nc_acs = await cursor.fetchall()
-        genomic_tx_data.alt_ac = nc_acs[0][0]
+            alt_ac = await uta.get_newest_assembly_ac(genomic_tx_data.alt_ac)
+            genomic_tx_data.alt_ac = alt_ac[0]
 
     def _set_liftover(
         self,
@@ -1280,7 +1255,7 @@ async def g_to_mane_c(
         coordinate_type = CoordinateType.INTER_RESIDUE
 
         async with self.uta_db.repository() as uta:
-            validation_result = uta.validate_genomic_ac(ac)
+            validation_result = await uta.validate_genomic_ac(ac)
         if not validation_result:
             _logger.warning("Genomic accession does not exist: %s", ac)
             return None

src/cool_seq_tool/schemas.py

@@ -1,7 +1,7 @@
 """Defines attribute constants, useful object structures, and API response schemas."""
 
 import datetime
-from enum import Enum, IntEnum
+from enum import IntEnum, StrEnum
 from typing import Literal
 
 from ga4gh.vrs.models import SequenceLocation
@@ -17,7 +17,7 @@
 _now = str(datetime.datetime.now(tz=datetime.UTC))
 
 
-class AnnotationLayer(str, Enum):
+class AnnotationLayer(StrEnum):
     """Create enum for supported annotation layers"""
 
     PROTEIN = "p"
@@ -32,7 +32,7 @@ class Strand(IntEnum):
     NEGATIVE = -1
 
 
-class Assembly(str, Enum):
+class Assembly(StrEnum):
     """Define supported genomic assemblies. Must be defined in ascending order"""
 
     GRCH37 = "GRCh37"
@@ -44,14 +44,14 @@ def values(cls) -> list[str]:
         return [item.value for item in cls]
 
 
-class ManeStatus(str, Enum):
+class ManeStatus(StrEnum):
     """Define constraints for mane status"""
 
     SELECT = "mane_select"
     PLUS_CLINICAL = "mane_plus_clinical"
 
 
-class TranscriptPriority(str, Enum):
+class TranscriptPriority(StrEnum):
     """Create Enum for Transcript Priority labels"""
 
     MANE_SELECT = ManeStatus.SELECT.value
@@ -60,7 +60,7 @@ class TranscriptPriority(str, Enum):
     GRCH38 = "grch38"
 
 
-class CoordinateType(str, Enum):
+class CoordinateType(StrEnum):
     """Create Enum for coordinate types.
 
     It is preferred to operate in inter-residue coordinates, but users should be

src/cool_seq_tool/sources/uta_database.py

@@ -113,6 +113,10 @@ async def execute_query(
         queries using the same DB connection. However, that means they are responsible
         for managing the cursor themselves.
 
+        For the sake of compactness and separation of concerns, other modules in CoolSeqTool
+        should avoid use of this method, and should instead add new methods to the repository
+        class itself.
+
         :param q: raw query. May need to specify schema depending on connection context.
         :param params: query variables, if needed. These should not be hard-coded into the query.
         :return: query result cursor
@@ -845,6 +849,73 @@ async def get_transcripts_from_genomic_pos(
         results = await cursor.fetchall()
         return [item for sublist in results for item in sublist]
 
+    async def get_alt_acs_for_tx(self, tx_ac: str) -> list[str]:
+        """Return genomic reference sequences associated with transcript accession
+
+        :param tx_ac: transcript accession
+        :return: list of genomic accessions for which alignments exist to ``tx_ac``
+        """
+        query = """
+            SELECT DISTINCT alt_ac
+            FROM exon_set
+            WHERE tx_ac = %(tx_ac)s
+              AND alt_aln_method != 'transcript';
+        """
+        cursor = await self.execute_query(query, {"tx_ac": tx_ac})
+        results = await cursor.fetchall()
+        return [r[0] for r in results]
+
+    async def get_gene_from_tx_ac(self, tx_ac: str) -> str | None:
+        """Return HGNC gene name for a transcript accession
+
+        :param tx_ac: transcript accession
+        :return: gene name, if found
+        """
+        query = """
+            SELECT hgnc
+            FROM tx_exon_aln_mv
+            WHERE tx_ac = %(tx_ac)s;
+        """
+        cursor = await self.execute_query(query, {"tx_ac": tx_ac})
+        result = await cursor.fetchone()
+        if result:
+            return result[0]
+        return None
+
+    async def validate_genomic_breakpoint(
+        self,
+        pos: int,
+        genomic_ac: str,
+        tx_ac: str,
+    ) -> bool:
+        """Validate that a genomic coordinate falls within the first and last exon
+            for a transcript on a given accession
+
+        :param pos: Genomic position on ``genomic_ac``
+        :param genomic_ac: RefSeq genomic accession, e.g. ``"NC_000007.14"``
+        :param transcript: A transcript accession
+        :return: ``True`` if the coordinate falls within 150bp of the first and last exon
+            for the transcript, ``False`` if not. Breakpoints past this threshold
+            are likely erroneous.
+        """
+        query = """
+            WITH tx_boundaries AS (
+                SELECT
+                MIN(alt_start_i) AS min_start,
+                MAX(alt_end_i) AS max_end
+                FROM tx_exon_aln_mv
+                WHERE tx_ac = %(tx_ac)s
+                AND alt_ac = %(genomic_ac)s
+            )
+            SELECT * FROM tx_boundaries
+            WHERE %(pos)s between (tx_boundaries.min_start - 150) and (tx_boundaries.max_end + 150);
+            """
+        cursor = await self.execute_query(
+            query, {"tx_ac": tx_ac, "genomic_ac": genomic_ac, "pos": pos}
+        )
+        result = await cursor.fetchone()
+        return bool(result)
+
 
 class ParseResult(UrlLibParseResult):
     """Subclass of url.ParseResult that adds database and schema methods, and provides stringification.

tests/sources/test_uta_database.py

@@ -366,6 +366,15 @@ async def test_get_alt_ac_start_or_end(
         await uta_repo.get_alt_ac_start_or_end("NM_152263.63", 822, 892, None)
 
 
+@pytest.mark.asyncio
+async def test_get_alt_acs_for_tx(uta_repo: UtaRepository):
+    resp = await uta_repo.get_alt_acs_for_tx("NM_004333.6")
+    assert set(resp) == {"NC_000007.13", "NC_000007.14", "NC_060931.1"}
+
+    resp = await uta_repo.get_alt_acs_for_tx("NM_99999999.99")
+    assert resp == []
+
+
 @pytest.mark.asyncio
 async def test_get_mane_transcripts_from_genomic_pos(uta_repo: UtaRepository):
     """Test that get_mane_transcripts_from_genomic_pos works correctly"""
@@ -398,6 +407,20 @@ async def test_get_mane_transcripts_from_genomic_pos(uta_repo: UtaRepository):
     assert resp == []
 
 
+@pytest.mark.asyncio
+async def test_get_gene_from_tx_ac(uta_repo: UtaRepository):
+    resp = await uta_repo.get_gene_from_tx_ac("NM_002529.3")
+    assert resp == "NTRK1"
+
+    resp = await uta_repo.get_gene_from_tx_ac("NM_99999999.99")
+    assert resp is None
+
+
+@pytest.mark.asyncio
+async def test_validate_genomic_breakpoint(uta_repo: UtaRepository):
+    raise NotImplementedError
+
+
 @pytest.mark.parametrize(
     ("raw_url", "expected"),
     [