refactor(htsget): use htslurp as htsget client

pyproject.toml

@@ -36,6 +36,7 @@ dependencies = [
     "zarr>=3.0.0",
     "platformdirs>=4.4.0",
     "obstore>=0.8.2",
+    "htslurp>=0.1.0",
 ]
 
 [dependency-groups]

src/modos/genomics/htsget.py

@@ -45,8 +45,11 @@
 from typing import Any
 from urllib.parse import urlparse, parse_qs
 
+import htslurp
+
 from crypt4gh import CIPHER_SEGMENT_SIZE
 from crypt4gh.lib import decrypt
+
 from pydantic import HttpUrl, validate_call
 from pydantic.dataclasses import dataclass
 import pysam
@@ -55,7 +58,7 @@
 from modos.remote import get_session
 from modos.genomics.c4gh import derive_public_key, get_secret_key
 from modos.genomics.region import Region
-from modos.genomics.formats import GenomicFileSuffix, read_pysam
+from modos.genomics.formats import GenomicFileSuffix
 
 
 @validate_call
@@ -337,39 +340,58 @@ def to_file(self, path: Path):
                 for block in stream:
                     sink.write(block)
 
+    @property
+    def format(self) -> str:
+        return GenomicFileSuffix.from_path(self.path).name
+
     @classmethod
     def from_url(cls, url: str):
         """Open connection directly from an htsget URL."""
         host, path, region = parse_htsget_url(url)
         return cls(host, path, region=region)
 
+    def records(self, reference: Path | None = None) -> htslurp.RecordIter:
+        # NOTE: Does note support crypt4gh encryption (yet)
+        records = htslurp.stream_records(
+            base_url=self.url,
+            id=str(self.path),
+            format=self.format,
+            region=self.region,
+            reference=reference,
+        )
+        return records
+
     def to_pysam(
-        self, reference_filename: str | None = None
+        self, reference_filename: Path | None = None
     ) -> Iterator[pysam.AlignedSegment | pysam.VariantRecord]:
         """Convert the stream to a pysam object."""
 
-        # NOTE: pysam needs a path or file descriptor,
-        # we have to stream from drive until this is addressed:
+        # NOTE: we use a dedicated client because pysam does not support bytestreams
         # ref: https://github.com/pysam-developers/pysam/blob/0787ca9da997b5911c00fd12584dad9741c82fb4/pysam/libcalignmentfile.pyx#L855
-        # TODO: when above addressed, replace temporary file with
+        # TODO: if above addressed, replace temporary file with
         # self.open() to stream directly from in-memory buffer.
-        buffer = tempfile.NamedTemporaryFile(
-            "w+b", delete=False, suffix="".join(self.path.suffixes)
-        ).name
 
-        self.to_file(Path(buffer))
-        buffer = read_pysam(
-            Path(buffer), reference_filename=reference_filename
-        )
+        stream = self.records(reference_filename)
+
+        for record in stream:
+            match self.format:
+                case "CRAM" | "BAM" | "SAM":
+                    parsed = pysam.AlignedSegment.fromstring(
+                        record.decode(), stream.header
+                    )
+                case _:
+                    # NOTE: pysam does not support instantiating VariantRecord on the fly.
+                    raise ValueError(
+                        f"Cannot convert {self.format} records to pysam."
+                    )
 
-        for record in buffer:
             if self.region is None:
-                yield record
+                yield parsed
                 continue
 
             # htsget includes all returns in the bgzf block
             # we filter out records outside requested region
             record_region = Region.from_pysam(record)
             if not record_region.overlaps(self.region):
                 continue
-            yield record
+            yield parsed