ADD --pam-only and --bgzip, for libraries short enough to match everywhere

CLAUDE.md

@@ -331,7 +331,10 @@ and narrowing it would buy nothing.
       has the comparison);
     - classes count sites whose genomic PAM matches. A site whose PAM does
       not stays in the table with `pam_matches` `FALSE` and counts nowhere,
-      so a spacer matching twice with one PAM is `unique`;
+      so a spacer matching twice with one PAM is `unique`. `--pam-only`
+      drops those rows instead, for a library short enough that they are
+      the bulk of the table: it changes what is written, never what is
+      searched or counted, so `n_matches` and `qc.tsv` are unmoved;
     - `run` writes only inside `outdir`, and only three folders: the results
       in `alignments/`, this invocation's record in `run_metadata/` and the
       log in `logs/`. Nothing of the run's own sits beside an input;

src/crispr_analysis_utils/cli/guide_alignment.py

@@ -109,6 +109,23 @@ def add_arguments(parser: argparse.ArgumentParser) -> None:
         action="store_true",
         help="The guide table's first line holds column names.",
     )
+    search.add_argument(
+        "--pam-only",
+        action="store_true",
+        help=(
+            "Write only the sites whose genomic PAM matches, which is every "
+            "site the nuclease could act at. A short spacer matches far more "
+            "places than it can cut; n_matches still counts them all."
+        ),
+    )
+    search.add_argument(
+        "--bgzip",
+        action="store_true",
+        help=(
+            "Write guide_mapping.bed.gz and its tabix index instead of a "
+            "plain BED, so a browser or tabix can seek into it."
+        ),
+    )
 
 
 def run(args: argparse.Namespace) -> int:
@@ -128,6 +145,8 @@ def run(args: argparse.Namespace) -> int:
             orientation_check=args.orientation_check,
             fix_orientation=args.fix_orientation,
             contigs=_contigs(args.contigs),
+            pam_only=args.pam_only,
+            bgzip=args.bgzip,
             id_col=_column(args.id_col),
             spacer_col=_column(args.spacer_col),
             sep=args.sep,

src/crispr_analysis_utils/guide_alignment/pipeline.py

@@ -7,7 +7,8 @@
 that.
 
 `outdir` gets three folders and nothing else: ``alignments/`` with the
-results (``guide_genome_annotation.tsv``, ``guide_mapping.bed``,
+results (``guide_genome_annotation.tsv``, ``guide_mapping.bed``, or
+``guide_mapping.bed.gz`` and its ``.tbi`` under ``bgzip``,
 ``guide_genome_annotation_columns.yaml``, ``qc.tsv``), ``run_metadata/`` with
 this invocation's record, and ``logs/``.
 """
@@ -97,6 +98,8 @@ def run(
     orientation_check: str = "error",
     fix_orientation: bool = False,
     contigs: Collection[str] | None = None,
+    pam_only: bool = False,
+    bgzip: bool = False,
     id_col: int | str = 0,
     spacer_col: int | str = 1,
     sep: str = "\t",
@@ -130,6 +133,13 @@ def run(
         instead of leaving them as given.
     contigs
         If given, read only these contigs.
+    pam_only
+        Write only the sites whose genomic PAM matches, which is every site
+        a nuclease could act at. A short spacer matches far more places than
+        it can cut, and the rest are what makes the table unusable.
+    bgzip
+        Write ``guide_mapping.bed.gz`` and its tabix index instead of a
+        plain BED, so a browser or `tabix` can seek into it.
     id_col, spacer_col, sep, header
         See `read_guides`.
 
@@ -173,6 +183,8 @@ def run(
         "orientation_check": orientation_check,
         "fix_orientation": fix_orientation,
         "contigs": None if contigs is None else list(contigs),
+        "pam_only": pam_only,
+        "bgzip": bgzip,
         "id_col": id_col,
         "spacer_col": spacer_col,
         "sep": sep,
@@ -252,6 +264,7 @@ def _run(
             reference,
             pam=parameters["pam"],
             contigs=parameters["contigs"],
+            pam_only=parameters["pam_only"],
         )
     state["counts"]["scan"] = scan.counts()
     write_json(record, state)
@@ -266,12 +279,18 @@ def _run(
             evidence=scan.evidence,
             pam=parameters["pam"],
             max_mismatches=0,
+            matches=scan.matches,
         )
         tables = outdir / RESULTS_DIR
         write_guides_tsv(
             tables / "guide_genome_annotation.tsv", library, summaries, scan.sites
         )
-        write_guides_bed(tables / "guide_mapping.bed", library, scan.sites)
+        write_guides_bed(
+            tables / "guide_mapping.bed",
+            library,
+            scan.sites,
+            bgzip=parameters["bgzip"],
+        )
         write_columns_yaml(tables / "guide_genome_annotation_columns.yaml")
         write_qc_tsv(tables / "qc.tsv", summaries, counts=state["counts"])
 
@@ -351,6 +370,7 @@ def _check_orientation(
             reference,
             pam=parameters["pam"],
             contigs=parameters["contigs"],
+            pam_only=parameters["pam_only"],
         )
     state["orientation_fixed"] = len(backwards)
     state["counts"]["scan"] = scan.counts()

src/crispr_analysis_utils/guide_alignment/scan.py

@@ -51,6 +51,10 @@ class ScanResult:
         For a ``reverse_complement`` row, whether the reversed reading was a
         distinct sequence or only the spacer's reverse complement; None for
         every other row.
+    matches
+        How many places each guide row's protospacer occurs, counted whether
+        or not the site was kept. Under `pam_only` it is the only record of
+        the sites that were not.
     n_matches
         Exact spacer matches found, before any were set aside.
     n_outside_contigs
@@ -59,14 +63,17 @@ class ScanResult:
     n_off_contig_end
         Matches whose PAM would run off the end of their contig.
     n_no_pam
-        Matches kept with a genomic PAM that matches no pattern.
+        Matches kept with a genomic PAM that matches no pattern. Under
+        `pam_only` these are counted and not written, so the table holds
+        ``n_sites - n_no_pam`` rows.
     n_sites
         Matches kept.
     """
 
     sites: dict[int, list[Site]] = field(default_factory=dict)
     verdicts: dict[int, str] = field(default_factory=dict)
     evidence: dict[int, str | None] = field(default_factory=dict)
+    matches: dict[int, int] = field(default_factory=dict)
     n_matches: int = 0
     n_outside_contigs: int = 0
     n_off_contig_end: int = 0
@@ -90,6 +97,7 @@ def scan_guides(
     *,
     pam: str = "NGG",
     contigs: Collection[str] | None = None,
+    pam_only: bool = False,
 ) -> ScanResult:
     """Find every exact site of every guide in `reference`.
 
@@ -103,6 +111,11 @@ def scan_guides(
         The PAM pattern, in IUPAC codes, 3' of the protospacer.
     contigs
         If given, read only these contigs.
+    pam_only
+        Keep only the matches whose genomic PAM matches `pam`. Every match
+        is still counted, in `ScanResult.matches` and `n_matches`; what
+        changes is which of them become a `Site`. A short spacer matches
+        everywhere, and the sites it cannot be cut at are the bulk of them.
 
     Returns
     -------
@@ -121,6 +134,9 @@ def scan_guides(
     Each distinct spacer is scanned once and its matches are given to every
     guide that uses it, so two entries differing only in an added G cost one
     search between them.
+
+    `pam_only` shrinks what is held and what is written, not what is
+    searched: every match is found and counted either way.
     """
     from .. import _core
 
@@ -156,27 +172,43 @@ def scan_guides(
     )
 
     result = ScanResult()
+    result.matches = dict.fromkeys(range(len(guides)), 0)
     reversed_has_pam: set[int] = set()
     for index, chrom, start, strand, genomic_pam, upstream, matched in matches:
         if index >= len(spacers):
+            # The core's verdict is enough here: this is only evidence that
+            # the reversed reading could cut somewhere, and no column of
+            # any table reports it.
             if matched:
                 reversed_has_pam.update(flipped_rows[queries[index]])
             continue
         result.n_matches += 1
         if len(genomic_pam) != len(primary):
             result.n_off_contig_end += 1
             continue
-        if not matched:
+        # The core's own `matched` is stricter than this over a genomic N,
+        # which `base_matches` reads as a PAM at an N of the pattern and
+        # the core reads as no PAM at all. The `pam_matches` column has
+        # always followed this one, so the counter and the flag do too:
+        # otherwise qc.tsv contradicts the table beside it.
+        has_pam = pattern_matches(genomic_pam, primary)
+        if not has_pam:
             result.n_no_pam += 1
         result.n_sites += 1
+        keep = has_pam or not pam_only
         for row in rows_of[spacers[index]]:
+            result.matches[row] += 1
+            if not keep:
+                continue
             guide = guides[row]
             result.sites.setdefault(row, []).append(
                 _site(guide, chrom, start, strand, genomic_pam, upstream, primary)
             )
 
     for row, guide in enumerate(guides):
-        verdict = _verdict(result.sites.get(row, ()), row in reversed_has_pam)
+        verdict = _verdict(
+            result.sites.get(row, ()), row in reversed_has_pam, result.matches[row]
+        )
         result.verdicts[row] = verdict
         result.evidence[row] = (
             None
@@ -230,10 +262,15 @@ def _kept(extra: str, genomic: str | None) -> str:
     )
 
 
-def _verdict(sites: Sequence[Site], reversed_has_pam: bool) -> str:
-    """One guide's orientation, from both directions' matches."""
+def _verdict(sites: Sequence[Site], reversed_has_pam: bool, n_matches: int) -> str:
+    """One guide's orientation, from both directions' matches.
+
+    The last branch counts matches rather than kept sites: under `pam_only`
+    a guide that matched everywhere without a PAM keeps none of them, and
+    reading the empty list would call it unmatched.
+    """
     if any(site.pam_class != NO_PAM_CLASS for site in sites):
         return FORWARD
     if reversed_has_pam:
         return REVERSE_COMPLEMENT
-    return NO_PAM if sites else NO_MATCH
+    return NO_PAM if n_matches else NO_MATCH

src/crispr_analysis_utils/guide_alignment/summary.py

@@ -15,10 +15,12 @@
 
 from __future__ import annotations
 
+import io
 import os
 from collections import Counter
 from collections.abc import Mapping, Sequence
 from dataclasses import dataclass, replace
+from pathlib import Path
 from typing import TYPE_CHECKING, TextIO
 
 from ..accession import guide_accession
@@ -139,6 +141,7 @@ def summarize(
     evidence: Mapping[int, str | None] | None = None,
     pam: str = "NGG",
     max_mismatches: int = 3,
+    matches: Mapping[int, int] | None = None,
 ) -> list[GuideSummary]:
     """Classify every guide from its sites.
 
@@ -155,6 +158,11 @@ def summarize(
         The PAM pattern the sites were classified with.
     max_mismatches
         The mismatch limit the sites were kept with.
+    matches
+        How many places each guide's protospacer occurs, as in
+        `ScanResult.matches`. None counts the sites given instead,
+        which is right only when every match became one: after a `pam_only`
+        scan it would count the PAM-bearing sites and call them every match.
 
     Returns
     -------
@@ -179,11 +187,11 @@ class that is neither one of the patterns nor ``none``, a guide is a
     per_spacer = Counter(guide.spacer for guide in guides)
     summaries = []
     for guide, source in zip(guides, _source_rows(guides), strict=True):
-        n_matches = 0
+        n_kept = 0
         n_sites = [0] * (max_mismatches + 1)
         perfect = []
         for site in sites.get(source, ()):
-            n_matches += 1
+            n_kept += 1
             if site.n_mismatches > max_mismatches:
                 raise ValueError(
                     f"A site of guide {site.guide_id!r} has {site.n_mismatches} "
@@ -200,6 +208,8 @@ class that is neither one of the patterns nor ``none``, a guide is a
                     f"{', '.join(patterns)}."
                 )
 
+        n_matches = n_kept if matches is None else matches.get(source, 0)
+
         if orientation is None:
             verdict = NOT_CHECKED
         elif source in orientation:
@@ -320,8 +330,9 @@ def fan_out_sites(
         "unique, multi or no_site: how many PAM-bearing places this guide matches."
     ),
     "n_matches": (
-        "How many places the protospacer matches, which is how many rows"
-        " this guide has below."
+        "How many places the protospacer matches, whether or not a PAM sits"
+        " beside it. That is one row each below, unless the run was given"
+        " --pam-only, which writes only the rows with a PAM."
     ),
     "n_target_sites": "How many of those have a PAM, so Cas9 can act there.",
     "orientation": (
@@ -410,7 +421,9 @@ def write_guides_tsv(
 
     ``n_matches`` is not ``site_class``: it counts every place the protospacer
     occurs, while ``multi`` counts only the matches with a PAM beside them, so
-    a guide can match many places and still bear a PAM at one.
+    a guide can match many places and still bear a PAM at one. It counts them
+    all after a `pam_only` scan too, which is the only place that count
+    survives.
     """
     if len(summaries) != len(guides):
         raise ValueError(f"Got {len(summaries)} summaries for {len(guides)} guides.")
@@ -452,7 +465,8 @@ def write_guides_bed(
     path: str | os.PathLike[str],
     guides: Sequence[Guide],
     sites: Mapping[int, Sequence[Site]],
-) -> None:
+    bgzip: bool = False,
+) -> Path:
     """Write guide_mapping.bed: every site a guide could act at.
 
     Parameters
@@ -463,6 +477,15 @@ def write_guides_bed(
         The library, as returned by `read_guides`.
     sites
         The sites by guide row, as in `ScanResult.sites`.
+    bgzip
+        Write `path` with ``.gz`` appended, block-compressed, and build its
+        tabix index beside it. The rows are already in coordinate order, so
+        nothing is sorted twice.
+
+    Returns
+    -------
+    pathlib.Path
+        The file written, which carries ``.gz`` under `bgzip`.
 
     Notes
     -----
@@ -488,9 +511,25 @@ def write_guides_bed(
             if site.pam_class != NO_PAM_CLASS
         )
     rows.sort()
-    with _open_table(path, None) as handle:
+    if not bgzip:
+        with _open_table(path, None) as handle:
+            for chrom, start, end, name, strand in rows:
+                _write_row(handle, [chrom, start, end, name, 0, strand])
+        return Path(path)
+
+    import pysam
+
+    # Straight into BGZF: a library short enough to need an index is one
+    # whose plain BED nobody wants on disk even for a moment.
+    path = Path(f"{os.fspath(path)}.gz")
+    with pysam.BGZFile(os.fspath(path), "wb") as raw:
+        handle = io.TextIOWrapper(raw, encoding="utf-8", newline="\n")
         for chrom, start, end, name, strand in rows:
             _write_row(handle, [chrom, start, end, name, 0, strand])
+        handle.flush()
+        handle.detach()
+    pysam.tabix_index(os.fspath(path), preset="bed", force=True)
+    return path
 
 
 def write_columns_yaml(path: str | os.PathLike[str]) -> None: