docs: find a gene's exons with bedspec and tidy the README's files and names

README.md

@@ -34,15 +34,15 @@ from captainslog.riker import ReadCategory
 from captainslog.riker import RikerReader
 
 with (
-    RikerReader.from_path[AlignmentMetrics]("tests/data/riker/s.alignment-metrics.txt") as alignments,
-    RikerReader.from_path[IsizeMetrics]("tests/data/riker/s.isize-metrics.txt") as isizes,
+    RikerReader.from_path[AlignmentMetrics]("riker/s.alignment-metrics.txt") as alignments,
+    RikerReader.from_path[IsizeMetrics]("riker/s.isize-metrics.txt") as isizes,
 ):
     alignment = {row.category: row for row in alignments}
-    [isize] = isizes
+    [isize_metrics] = isizes
 
-frac_reads_in_fr_pairs = 2 * isize.read_pairs / alignment[ReadCategory.PAIR].total_reads
+frac_fr_reads = 2 * isize_metrics.read_pairs / alignment[ReadCategory.PAIR].total_reads
 
-assert round(frac_reads_in_fr_pairs, 3) == 0.972
+assert round(frac_fr_reads, 3) == 0.972
 ```
 
 ## Quickstart
@@ -56,19 +56,17 @@ from captainslog import riker
 from captainslog.riker import HybcapMetrics
 from captainslog.riker import IsizeHistogram
 
-assert riker.metrics_class("tests/data/riker/s.hybcap-metrics.txt") is HybcapMetrics
-
-with RikerReader.from_path[HybcapMetrics]("tests/data/riker/s.hybcap-metrics.txt") as reader:
+with RikerReader.from_path[HybcapMetrics]("riker/s.hybcap-metrics.txt") as reader:
     [hybcap] = reader
 
 assert hybcap.sample == "NA12891"
 assert hybcap.mean_target_coverage == 70.54
 
-with RikerReader.from_path[IsizeHistogram]("tests/data/riker/s.isize-histogram.txt") as reader:
-    isize = riker.columns(list(reader))
+with RikerReader.from_path[IsizeHistogram]("riker/s.isize-histogram.txt") as reader:
+    isize_histogram = riker.columns(list(reader))
 
-assert isize.insert_size[:3] == (112, 128, 134)
-assert isize.fr_count[:3] == (1, 1, 1)
+assert isize_histogram.insert_size[:3] == (112, 128, 134)
+assert isize_histogram.fr_count[:3] == (1, 1, 1)
 ```
 
 ### Reading Picard Metrics
@@ -79,7 +77,7 @@ A `PicardReader` reads a Picard metrics file into dataclasses, one per row.
 from captainslog.picard import HsMetrics
 from captainslog.picard import PicardReader
 
-with PicardReader.from_path[HsMetrics]("tests/data/picard/hs.hs_metrics") as reader:
+with PicardReader.from_path[HsMetrics]("picard/hs.hs_metrics.txt") as reader:
     [hs_metrics] = reader
 
 assert hs_metrics.mean_target_coverage == 70.54
@@ -108,7 +106,7 @@ from tempfile import TemporaryDirectory
 from captainslog.picard import PicardWriter
 
 temp_dir = TemporaryDirectory()
-path = Path(temp_dir.name) / "s.hs_metrics"
+path = Path(temp_dir.name) / "s.hs_metrics.txt"
 
 with PicardWriter.from_path[HsMetrics](path) as writer:
     writer.write(hs_metrics)
@@ -121,16 +119,17 @@ with PicardReader.from_path[HsMetrics](path) as reader:
 
 #### Finding the Class of a File
 
-The `metrics_class` and `histogram_class` functions return the metric type in a data file.
+The `metrics_class` function returns the metric type in a data file, and Picard's `histogram_class` returns the type of the histogram inside one.
 
 ```python
 from captainslog import picard
 from captainslog.picard import HsHistogram
-from captainslog.picard import MeanQualityByCycle
 
-assert picard.metrics_class("tests/data/picard/hs.hs_metrics") is HsMetrics
-assert picard.histogram_class("tests/data/picard/hs.hs_metrics") is HsHistogram
-assert picard.histogram_class("tests/data/picard/cmm.quality_by_cycle_metrics") is MeanQualityByCycle
+assert riker.metrics_class("riker/s.hybcap-metrics.txt") is HybcapMetrics
+assert riker.metrics_class("riker/s.isize-histogram.txt") is IsizeHistogram
+
+assert picard.metrics_class("picard/hs.hs_metrics.txt") is HsMetrics
+assert picard.histogram_class("picard/hs.hs_metrics.txt") is HsHistogram
 ```
 
 #### Histograms and Columns
@@ -147,25 +146,26 @@ from captainslog.picard import InsertSizeHistogram
 from captainslog.picard import PairOrientation
 from captainslog.picard import PicardHistogramReader
 
-with PicardHistogramReader.from_path[InsertSizeHistogram]("tests/data/picard/cmm.insert_size_metrics") as reader:
-    hist = picard.columns(list(reader))
+with PicardHistogramReader.from_path[InsertSizeHistogram]("picard/cmm.insert_size_metrics.txt") as reader:
+    insert_size_histogram = picard.columns(list(reader))
 
-assert hist.insert_size[:3] == (112, 128, 134)
+assert insert_size_histogram.insert_size[:3] == (112, 128, 134)
 
-series = picard.insert_size_series(hist)
+series = picard.insert_size_series(insert_size_histogram)
 all_reads = next(s for s in series if s.name is None and s.orientation is PairOrientation.FR)
 
 assert all_reads.counts[:3] == (1, 1, 1)
 assert sum(all_reads.counts) == 478
-assert [s.name for s in series] == [None, "NA12891", "Pond-492101", "NormPond-544932", "NA12891.r1", "NA12891.r2"]
+names = [s.name for s in series]
+assert names == [None, "NA12891", "Pond-492101", "NormPond-544932", "NA12891.r1", "NA12891.r2"]
 ```
 
 Histogram-only files read the same way:
 
 ```python
 from captainslog.picard import QualityScoreDistribution
 
-with PicardHistogramReader.from_path[QualityScoreDistribution]("tests/data/picard/cmm.quality_distribution_metrics") as reader:
+with PicardHistogramReader.from_path[QualityScoreDistribution]("picard/cmm.quality_distribution_metrics.txt") as reader:
     qualities = list(reader)
 
 assert qualities[0].quality == 1
@@ -180,7 +180,7 @@ Picard writes a row for each accumulation level: all reads, then each sample, li
 from captainslog.picard import AccumulationLevel
 from captainslog.picard import InsertSizeMetrics
 
-with PicardReader.from_path[InsertSizeMetrics]("tests/data/picard/cmm.insert_size_metrics") as reader:
+with PicardReader.from_path[InsertSizeMetrics]("picard/cmm.insert_size_metrics.txt") as reader:
     rows = list(reader)
 
 assert [picard.accumulation_level(row) for row in rows] == [
@@ -206,7 +206,7 @@ CollectHsMetrics and CollectTargetedPcrMetrics can also write plain tables of co
 from captainslog.picard import PerTargetCoverage
 from captainslog.picard import PicardCoverageReader
 
-with PicardCoverageReader.from_path[PerTargetCoverage]("tests/data/picard/hs.per_target_coverage.txt") as reader:
+with PicardCoverageReader.from_path[PerTargetCoverage]("picard/hs.per_target_coverage.txt") as reader:
     targets = list(reader)
 
 assert targets[0].name == "rs7555566"
@@ -223,12 +223,36 @@ class MyHsMetrics(HsMetrics):
     def frac_bases_on_target(self) -> float:
         return self.on_target_bases / self.pf_bases_aligned
 
-with PicardReader.from_path[MyHsMetrics]("tests/data/picard/hs.hs_metrics") as reader:
+
+with PicardReader.from_path[MyHsMetrics]("picard/hs.hs_metrics.txt") as reader:
     [hs_metrics] = reader
 
 assert round(hs_metrics.frac_bases_on_target, 4) == 0.2398
 ```
 
+## Intervals with bedspec
+
+With [bedspec](https://github.com/clintval/bedspec), find captainslog's targets, bases, SNPs, and CpGs by overlap:
+
+```python
+from bedspec import Bed3
+from bedspec import register_locatable
+from bedspec.ext.captainslog import as_interval
+from bedspec.overlap import TreeDetector
+
+gene = Bed3("chr1", start=500, end=2_000)
+
+targets = [
+    PerTargetCoverage(chrom="chr1", start=1_001, end=1_100, name="exon2"),
+    PerTargetCoverage(chrom="chr1", start=1_501, end=1_600, name="exon3"),
+    PerTargetCoverage(chrom="chr2", start=501, end=600, name="exon7"),
+]
+
+with register_locatable(PerTargetCoverage, as_interval):
+    exons = TreeDetector(targets).overlapping(gene)
+    assert {exon.name for exon in exons} == {"exon2", "exon3"}
+```
+
 ## Development and Testing
 
 See the [contributing guide](./CONTRIBUTING.md) for more information.

conftest.py

@@ -1,7 +1,17 @@
+from pathlib import Path
+
+import pytest
 from sybil import Sybil
 from sybil.parsers.markdown import PythonCodeBlockParser
 from sybil.parsers.markdown import SkipParser
 
+
+@pytest.fixture
+def in_test_data(monkeypatch: pytest.MonkeyPatch) -> None:
+    """Run each README example from the test data, so it names files as a user would."""
+    monkeypatch.chdir(Path(__file__).parent / "tests" / "data")
+
+
 pytest_collect_file = Sybil(
-    parsers=[PythonCodeBlockParser(), SkipParser()], patterns=["*.md"]
+    parsers=[PythonCodeBlockParser(), SkipParser()], patterns=["*.md"], fixtures=["in_test_data"]
 ).pytest()

pyproject.toml

@@ -37,6 +37,7 @@ classifiers = [
 
 [dependency-groups]
 test = [
+    "bedspec ~=4.3",
     "pytest ~=9.1",
     "pytest-cov ~=7.1",
     "sybil[pytest] ~=10.1",
@@ -81,8 +82,8 @@ only-include = [
 # caught before it can be locked. Must stay >= the Dependabot `cooldown` in
 # .github/dependabot.yml, or Dependabot will propose versions that can't lock.
 exclude-newer = "1 week"
-# typeline is released alongside captainslog, so its new releases are allowed right away.
-exclude-newer-package = { typeline = false }
+# bedspec and typeline are released alongside captainslog, so their new releases are allowed right away.
+exclude-newer-package = { bedspec = false, typeline = false }
 
 [tool.poe.tasks]
 fix-format = "ruff format"

scripts/fixtures/make.sh

@@ -71,9 +71,9 @@ levels=(
   METRIC_ACCUMULATION_LEVEL=READ_GROUP
 )
 
-picard MarkDuplicates I=in/merged.bam O="$bam" M=picard/dup.duplication_metrics
+picard MarkDuplicates I=in/merged.bam O="$bam" M=picard/dup.duplication_metrics.txt
 samtools index "$bam"
-picard MarkDuplicates I=in/r1.bam O=in/r1.marked.bam M=picard/one_library.duplication_metrics
+picard MarkDuplicates I=in/r1.bam O=in/r1.marked.bam M=picard/one_library.duplication_metrics.txt
 
 picard CollectMultipleMetrics I="$bam" O=picard/cmm R="$ref" "${levels[@]}" \
   PROGRAM=null \
@@ -86,27 +86,27 @@ picard CollectMultipleMetrics I="$bam" O=picard/cmm R="$ref" "${levels[@]}" \
   PROGRAM=MeanQualityByCycle \
   PROGRAM=QualityScoreDistribution
 # Picard 3.5.0 shares target depth histograms across levels and throws on any level but ALL_READS.
-picard CollectHsMetrics I="$bam" O=picard/hs.hs_metrics R="$ref" \
+picard CollectHsMetrics I="$bam" O=picard/hs.hs_metrics.txt R="$ref" \
   BAIT_INTERVALS=in/baits.interval_list \
   TARGET_INTERVALS=in/targets.interval_list \
   PER_TARGET_COVERAGE=picard/hs.per_target_coverage.txt \
   PER_BASE_COVERAGE=picard/hs.per_base_coverage.txt
-picard CollectTargetedPcrMetrics I="$bam" O=picard/pcr.targeted_pcr_metrics R="$ref" \
+picard CollectTargetedPcrMetrics I="$bam" O=picard/pcr.targeted_pcr_metrics.txt R="$ref" \
   AMPLICON_INTERVALS=in/baits.interval_list \
   TARGET_INTERVALS=in/targets.interval_list
-picard CollectWgsMetrics I="$bam" O=picard/wgs.wgs_metrics R="$ref"
-picard CollectRawWgsMetrics I="$bam" O=picard/raw.raw_wgs_metrics R="$ref" COVERAGE_CAP=250
-picard CollectWgsMetricsWithNonZeroCoverage I="$bam" O=picard/nonzero.wgs_metrics R="$ref" \
+picard CollectWgsMetrics I="$bam" O=picard/wgs.wgs_metrics.txt R="$ref"
+picard CollectRawWgsMetrics I="$bam" O=picard/raw.raw_wgs_metrics.txt R="$ref" COVERAGE_CAP=250
+picard CollectWgsMetricsWithNonZeroCoverage I="$bam" O=picard/nonzero.wgs_metrics.txt R="$ref" \
   CHART=picard/nonzero.wgs_metrics.pdf
-picard CollectRnaSeqMetrics I="$bam" O=picard/rna.rna_metrics "${levels[@]}" \
+picard CollectRnaSeqMetrics I="$bam" O=picard/rna.rna_metrics.txt "${levels[@]}" \
   REF_FLAT=in/genes.refflat \
   RIBOSOMAL_INTERVALS=in/ribosomal.interval_list \
   STRAND_SPECIFICITY=NONE \
   MINIMUM_LENGTH=100
-picard CollectOxoGMetrics I="$bam" O=picard/oxog.oxog_metrics R="$ref"
-picard EstimateLibraryComplexity I="$bam" O=picard/elc.duplication_metrics
-picard CollectJumpingLibraryMetrics I="$bam" O=picard/jump.jumping_metrics
-picard CrosscheckFingerprints O=picard/crosscheck.crosscheck_metrics HAPLOTYPE_MAP="$sites" \
+picard CollectOxoGMetrics I="$bam" O=picard/oxog.oxog_metrics.txt R="$ref"
+picard EstimateLibraryComplexity I="$bam" O=picard/elc.duplication_metrics.txt
+picard CollectJumpingLibraryMetrics I="$bam" O=picard/jump.jumping_metrics.txt
+picard CrosscheckFingerprints O=picard/crosscheck.crosscheck_metrics.txt HAPLOTYPE_MAP="$sites" \
   I="$bam" \
   I=in/NA12892.over.fingerprints.shifted.for.crams.r1.sam \
   CROSSCHECK_BY=LIBRARY
@@ -134,6 +134,10 @@ riker multi -i "$bam" -o riker/s -r "$ref" \
   --rna::ribosomal-intervals in/ribosomal.interval_list
 
 find picard riker \( -name '*.pdf' -o -name '*.bw' \) -delete
+# Some Picard tools take only an output prefix, so their outputs are given the .txt extension here.
+for path in picard/*; do
+  [[ "$path" == *.txt ]] || command mv "$path" "$path.txt"
+done
 command rm -rf "$repo/tests/data/picard" "$repo/tests/data/riker"
 mkdir -p "$repo/tests/data"
 command cp -Rf picard riker "$repo/tests/data/"

tests/picard/__init__.py

@@ -69,19 +69,19 @@ def _holds_a_histogram(path: Path) -> bool:
 PICARD_HISTOGRAM_FIXTURES = sorted(path for path in DATA.glob("*") if _holds_a_histogram(path))
 
 FIXTURE_CLASSES: dict[str, type[Any]] = {
-    "cmm.alignment_summary_metrics": ReadLengthHistogram,
-    "cmm.insert_size_metrics": InsertSizeHistogram,
-    "cmm.quality_by_cycle_metrics": MeanQualityByCycle,
-    "cmm.quality_distribution_metrics": QualityScoreDistribution,
-    "dup.duplication_metrics": DuplicateSetSizeHistogram,
-    "elc.duplication_metrics": LibraryComplexityHistogram,
-    "hs.hs_metrics": HsHistogram,
-    "nonzero.wgs_metrics": WgsNonZeroCoverageHistogram,
-    "one_library.duplication_metrics": DuplicationRoiHistogram,
-    "pcr.targeted_pcr_metrics": HsHistogram,
-    "raw.raw_wgs_metrics": WgsHistogram,
-    "rna.rna_metrics": RnaCoverageHistogram,
-    "wgs.wgs_metrics": WgsHistogram,
+    "cmm.alignment_summary_metrics.txt": ReadLengthHistogram,
+    "cmm.insert_size_metrics.txt": InsertSizeHistogram,
+    "cmm.quality_by_cycle_metrics.txt": MeanQualityByCycle,
+    "cmm.quality_distribution_metrics.txt": QualityScoreDistribution,
+    "dup.duplication_metrics.txt": DuplicateSetSizeHistogram,
+    "elc.duplication_metrics.txt": LibraryComplexityHistogram,
+    "hs.hs_metrics.txt": HsHistogram,
+    "nonzero.wgs_metrics.txt": WgsNonZeroCoverageHistogram,
+    "one_library.duplication_metrics.txt": DuplicationRoiHistogram,
+    "pcr.targeted_pcr_metrics.txt": HsHistogram,
+    "raw.raw_wgs_metrics.txt": WgsHistogram,
+    "rna.rna_metrics.txt": RnaCoverageHistogram,
+    "wgs.wgs_metrics.txt": WgsHistogram,
 }
 
 

tests/picard/test_classes.py

@@ -24,7 +24,7 @@ def test_metrics_class_names_the_dataclass_of_every_fixture(path: Path) -> None:
 
 def test_metrics_class_and_the_reader_read_a_gzipped_file(tmp_path: Path) -> None:
     """A gzipped file gives the same metrics class and the same rows as its plain file."""
-    source = DATA / "cmm.quality_yield_metrics"
+    source = DATA / "cmm.quality_yield_metrics.txt"
     path = tmp_path / "zipped.quality_yield_metrics.gz"
     _ = path.write_bytes(gzip.compress(source.read_bytes()))
     assert metrics_class(path) is QualityYieldMetrics
@@ -34,12 +34,12 @@ def test_metrics_class_and_the_reader_read_a_gzipped_file(tmp_path: Path) -> Non
 
 def test_metrics_class_accepts_a_path_given_as_text() -> None:
     """A path given as text gives the same metrics class as a `Path`."""
-    assert metrics_class(str(DATA / "cmm.quality_yield_metrics")) is QualityYieldMetrics
+    assert metrics_class(str(DATA / "cmm.quality_yield_metrics.txt")) is QualityYieldMetrics
 
 
 def test_metrics_class_points_a_histogram_only_file_to_histogram_class() -> None:
     """A file holding only a histogram is refused, pointing to `histogram_class`."""
-    path = DATA / "cmm.quality_by_cycle_metrics"
+    path = DATA / "cmm.quality_by_cycle_metrics.txt"
     message = (
         rf"^{re.escape(str(path))} holds no metrics table! "
         r"If it holds only a histogram, find its type with histogram_class\.$"
@@ -108,7 +108,7 @@ def test_metrics_class_refuses_a_file_whose_lines_end_with_a_lone_carriage_retur
     """A file whose lines end with a lone carriage return is refused by both."""
     path = tmp_path / "mac.metrics"
     _ = path.write_text(
-        (DATA / "cmm.quality_yield_metrics").read_text().replace("\n", "\r"), newline=""
+        (DATA / "cmm.quality_yield_metrics.txt").read_text().replace("\n", "\r"), newline=""
     )
     message = (
         rf"^{re.escape(str(path))} breaks its lines with lone carriage returns, not line feeds!$"

tests/picard/test_format.py

@@ -116,7 +116,7 @@ def test_picard_float_codec_reads_numbers_as_typeline_does(text: str, value: flo
 
 @pytest.mark.parametrize(
     ("name", "kind"),
-    [("hs.hs_metrics", HsMetrics), ("cmm.insert_size_metrics", InsertSizeMetrics)],
+    [("hs.hs_metrics.txt", HsMetrics), ("cmm.insert_size_metrics.txt", InsertSizeMetrics)],
     ids=["hs", "insert-size"],
 )
 def test_reader_reads_metrics_whose_floats_hold_a_decimal_comma(
@@ -140,8 +140,8 @@ def test_reader_reads_metrics_whose_floats_hold_a_decimal_comma(
 @pytest.mark.parametrize(
     ("name", "kind"),
     [
-        ("rna.rna_metrics", RnaCoverageHistogram),
-        ("cmm.quality_by_cycle_metrics", MeanQualityByCycle),
+        ("rna.rna_metrics.txt", RnaCoverageHistogram),
+        ("cmm.quality_by_cycle_metrics.txt", MeanQualityByCycle),
     ],
     ids=["rna", "mean-quality"],
 )