feat: find fgumi and fgbio variant reviews with bedspec.ext.captainslog

README.md

@@ -442,7 +442,7 @@ Each module names them the same way: `as_interval` finds a record by its interva
 | --- | --- | --- |
 | `bedspec.ext.biopython` | `as_interval` | a `Bio.Align.Alignment`, like a BED12 line, by its first to last aligned base, on its strand |
 | | `as_blocks` | a `Bio.Align.Alignment` by its blocks, on its strand |
-| `bedspec.ext.captainslog` | `as_interval` | a Picard or Riker per-target or per-base coverage row, or a `FingerprintingDetailMetrics` SNP, converted from the 1-based coordinates it holds; a `RrbsCpgDetailMetrics` by the C of its CpG on the forward strand |
+| `bedspec.ext.captainslog` | `as_interval` | a Riker or Picard per-target or per-base coverage row, a `FingerprintingDetailMetrics` SNP, or an fgumi or fgbio `ConsensusVariantReviewInfo` variant, converted from the 1-based coordinates it holds; a `RrbsCpgDetailMetrics` by the C of its CpG on the forward strand |
 | `bedspec.ext.cyvcf2` | `as_interval` | a `Variant` from its start to its end, so a symbolic allele to its `END` |
 | `bedspec.ext.fgpyo` | `as_interval` | a `SupplementaryAlignment` from an `SA` tag, on its strand |
 | | `as_template_intervals` | a `Template` by the intervals of its primary and supplementary records, each on its strand |

bedspec/ext/_captainslog.py

@@ -1,6 +1,8 @@
 from typing import overload
 
 try:
+    from captainslog import fgbio
+    from captainslog import fgumi
     from captainslog.picard import FingerprintingDetailMetrics
     from captainslog.picard import PerBaseCoverage
     from captainslog.picard import PerTargetCoverage
@@ -28,17 +30,22 @@ def as_interval(record: HybcapPerBase) -> list[Bed4]: ...
 @overload
 def as_interval(record: FingerprintingDetailMetrics) -> list[Bed4]: ...
 @overload
+def as_interval(record: fgbio.ConsensusVariantReviewInfo) -> list[Bed3]: ...
+@overload
+def as_interval(record: fgumi.ConsensusVariantReviewInfo) -> list[Bed3]: ...
+@overload
 def as_interval(record: RrbsCpgDetailMetrics) -> list[Bed3]: ...
 # fmt: on
 def as_interval(record: object) -> list[Bed4] | list[Bed3]:
-    """Return the interval of a target, a base, a SNP, or a CpG, from the coordinates it holds.
+    """Return the interval of a target, a base, a SNP, a variant, or a CpG, from its coordinates.
 
     A target runs from its 1-based start to its end and a base or a SNP covers its 1-based position,
     each named for its target or SNP, while a CpG covers its C on the forward strand, at the 0-based
-    position Picard writes.
+    position Picard writes. A variant fgumi or fgbio reviewed covers its 1-based position, unnamed,
+    as its row names a consensus read and not the variant.
 
     Raises:
-        TypeError: If the record is not one of captainslog's targets, bases, SNPs, or CpGs.
+        TypeError: If the record is none of captainslog's targets, bases, SNPs, variants, or CpGs.
         ValueError: If the record is a CpG without a position.
     """
     if isinstance(record, PerTargetCoverage | HybcapPerTarget):
@@ -50,6 +57,8 @@ def as_interval(record: object) -> list[Bed4] | list[Bed3]:
     if isinstance(record, FingerprintingDetailMetrics):
         name = record.snp or None
         return [Bed4(record.chrom, start=record.position - 1, end=record.position, name=name)]
+    if isinstance(record, fgbio.ConsensusVariantReviewInfo | fgumi.ConsensusVariantReviewInfo):
+        return [Bed3(record.chrom, start=record.pos - 1, end=record.pos)]
     if isinstance(record, RrbsCpgDetailMetrics):
         if record.position is None:
             raise ValueError("A CpG needs a position to find its interval, not None!")

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "bedspec"
-version = "4.3.0"
+version = "4.4.0"
 description =  "An HTS-specs compliant BED toolkit."
 authors = [{ name = "Clint Valentine", email = "[email protected]" }]
 license = "MIT"
@@ -39,7 +39,7 @@ classifiers = [
 
 [project.optional-dependencies]
 biopython   = ["biopython >=1.88"]
-captainslog = ["captainslog-ngs >=0.1.0"]
+captainslog = ["captainslog-ngs >=0.2.0"]
 cyvcf2      = ["cyvcf2 >=0.34"]
 fgpyo       = ["fgpyo >=1.5.1"]
 mappy       = ["mappy >=2.31"]

tests/ext/test_captainslog.py

@@ -4,6 +4,11 @@
 from typing import Any
 
 import pytest
+from captainslog import fgbio
+from captainslog import fgumi
+from captainslog.fgbio import ErrorRateByReadPositionMetric
+from captainslog.fgbio import FgbioReader
+from captainslog.fgumi import FgumiReader
 from captainslog.picard import FingerprintingDetailMetrics
 from captainslog.picard import PerBaseCoverage
 from captainslog.picard import PerTargetCoverage
@@ -73,6 +78,14 @@
 CHR1 = "TTCATGCTGAAGCCCTCTTACGATCGTACAGATGCAAATATTAACAAACC"
 """The first bases of chr1 of the reference the RRBS rows were collected against."""
 
+REVIEW = """chrom\tpos\tref\tgenotype\tfilters\tA\tC\tG\tT\tN\tconsensus_read\tconsensus_insert\tconsensus_call\tconsensus_qual\ta\tc\tg\tt\tn
+chr1\t501\tA\tA/G\tPASS\t2\t0\t1\t0\t0\tS1-lib|S2-lib:4/1\tchr1:452-770 | F1R2\tG\t72\t0\t0\t1\t0\t0
+"""  # noqa: E501
+"""A consensus read carrying a variant, as fgumi and fgbio both write it."""
+
+CHR1_FROM_490 = "TCATAATCTCAGTGCTGGTG"
+"""The bases of chr1 from 0-based 490 of the reference the reviewed reads were aligned to."""
+
 
 def written(tmp_path: Path, name: str, text: str) -> Path:
     """Return the path of a file in a temporary directory holding text."""
@@ -147,6 +160,22 @@ def test_as_interval_of_a_cpg_covers_its_c_at_the_0_based_position_picard_writes
     assert CHR1[spans[0].start : spans[0].end + 1] == "CG"
 
 
+def test_as_interval_of_a_reviewed_variant_covers_its_1_based_position(tmp_path: Path) -> None:
+    """Test that a variant fgumi or fgbio reviewed covers its position, 0-based, and is unnamed."""
+    path = written(tmp_path, "review.txt", REVIEW)
+    with (
+        FgbioReader.from_path[fgbio.ConsensusVariantReviewInfo](path) as by_fgbio,
+        FgumiReader.from_path[fgumi.ConsensusVariantReviewInfo](path) as by_fgumi,
+    ):
+        variants: list[fgbio.ConsensusVariantReviewInfo | fgumi.ConsensusVariantReviewInfo] = [
+            *by_fgbio,
+            *by_fgumi,
+        ]
+    spans = [span for variant in variants for span in as_interval(variant)]
+    assert spans == 2 * [Bed3("chr1", start=500, end=501)]
+    assert CHR1_FROM_490[spans[0].start - 490 : spans[0].end - 490] == variants[0].ref
+
+
 def test_as_interval_of_a_target_or_base_without_a_name_is_unnamed() -> None:
     """Test that a target or a base whose name is empty has an interval without a name."""
     assert as_interval(PerTargetCoverage(chrom="1", start=201, end=300)) == [
@@ -196,13 +225,34 @@ def test_registered_targets_are_found_by_their_first_and_last_bases_only(tmp_pat
         assert [target.name for target in everywhere] == ["rs17272796"]
 
 
+def test_registered_reviewed_variants_are_found_by_their_base_only(tmp_path: Path) -> None:
+    """Test that registered reviewed variants are found by their base, not the bases beside it."""
+    path = written(tmp_path, "review.txt", REVIEW)
+    with (
+        FgumiReader.from_path[fgumi.ConsensusVariantReviewInfo](path) as reader,
+        register_locatable(fgumi.ConsensusVariantReviewInfo, as_interval),
+    ):
+        variants = TreeDetector(reader)
+        found = variants.overlapping(Bed3("chr1", start=500, end=501))
+        assert [variant.consensus_read for variant in found] == ["S1-lib|S2-lib:4/1"]
+        assert not variants.overlaps(Bed3("chr1", start=499, end=500))
+        assert not variants.overlaps(Bed3("chr1", start=501, end=502))
+
+
 def test_as_interval_of_a_record_of_another_class_raises_naming_it() -> None:
     """Test that a record of a class without an interval, like Riker's hybcap summary, raises."""
     record: Any = HybcapMetrics()
     with pytest.raises(TypeError, match=r"of type captainslog\.riker\..*HybcapMetrics!"):
         as_interval(record)
 
 
+def test_as_interval_of_an_error_rate_by_read_position_raises() -> None:
+    """Test that fgbio's error rate at a position within a read, not on a reference, raises."""
+    record: Any = ErrorRateByReadPositionMetric(position=501)
+    with pytest.raises(TypeError, match=r"captainslog\.fgbio\..*ErrorRateByReadPositionMetric!"):
+        as_interval(record)
+
+
 def test_as_interval_of_a_cpg_without_a_position_raises() -> None:
     """Test that a CpG made without a position has no interval, and raises saying so."""
     with pytest.raises(ValueError, match="needs a position"):

tests/test_ext_imports.py

@@ -13,7 +13,13 @@
     "pybedtools": ("pybedtools",),
     "pybwa": ("pybwa",),
     "pysam": ("pysam", "pysam.libctabixproxies"),
-    "captainslog": ("captainslog", "captainslog.picard", "captainslog.riker"),
+    "captainslog": (
+        "captainslog",
+        "captainslog.fgbio",
+        "captainslog.fgumi",
+        "captainslog.picard",
+        "captainslog.riker",
+    ),
 }
 """The modules each extra installs, blocked so importing them fails as if they were missing."""