From 3453e9718c01aee9d17a2503ba17113a9b65aa4f Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 9 Nov 2026 11:06:55 -0500 Subject: [PATCH 1/4] ADD --pam-only, to write just the sites a nuclease could act at A spacer's matches grow as the genome over four to its length: a 9-mer matches tens of thousands of places and can be cut at about one in sixteen of them, so the table is mostly rows nothing can be done with. The sites are dropped as the scan reads them, so neither the table nor the memory holds them, and every match is still counted: ScanResult gains a per-row count that summarize takes, so n_matches keeps meaning every place the protospacer occurs. The orientation verdict reads that count too -- a guide that matched everywhere without a PAM keeps none of its sites, and the empty list would otherwise call it unmatched. Co-Authored-By: Claude Opus 5 --- CLAUDE.md | 5 +- .../cli/guide_alignment.py | 10 ++ .../guide_alignment/pipeline.py | 9 ++ .../guide_alignment/scan.py | 37 +++++- .../guide_alignment/summary.py | 21 +++- .../skills/cau-guide-alignment/SKILL.md | 25 ++++ .../cau-guide-alignment/references/outputs.md | 10 +- tests/test_cli_guide_alignment.py | 3 + tests/test_guide_alignment_pipeline.py | 112 ++++++++++++++++++ 9 files changed, 218 insertions(+), 14 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index b7e94d7..db5a5c8 100644 --- a/CLAUDE.md +++ b/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; diff --git a/src/crispr_analysis_utils/cli/guide_alignment.py b/src/crispr_analysis_utils/cli/guide_alignment.py index 933ad0a..98d1e05 100644 --- a/src/crispr_analysis_utils/cli/guide_alignment.py +++ b/src/crispr_analysis_utils/cli/guide_alignment.py @@ -109,6 +109,15 @@ 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." + ), + ) def run(args: argparse.Namespace) -> int: @@ -128,6 +137,7 @@ 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, id_col=_column(args.id_col), spacer_col=_column(args.spacer_col), sep=args.sep, diff --git a/src/crispr_analysis_utils/guide_alignment/pipeline.py b/src/crispr_analysis_utils/guide_alignment/pipeline.py index ab40c7f..936b71d 100644 --- a/src/crispr_analysis_utils/guide_alignment/pipeline.py +++ b/src/crispr_analysis_utils/guide_alignment/pipeline.py @@ -97,6 +97,7 @@ def run( orientation_check: str = "error", fix_orientation: bool = False, contigs: Collection[str] | None = None, + pam_only: bool = False, id_col: int | str = 0, spacer_col: int | str = 1, sep: str = "\t", @@ -130,6 +131,10 @@ 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. id_col, spacer_col, sep, header See `read_guides`. @@ -173,6 +178,7 @@ def run( "orientation_check": orientation_check, "fix_orientation": fix_orientation, "contigs": None if contigs is None else list(contigs), + "pam_only": pam_only, "id_col": id_col, "spacer_col": spacer_col, "sep": sep, @@ -252,6 +258,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,6 +273,7 @@ def _run( evidence=scan.evidence, pam=parameters["pam"], max_mismatches=0, + matches=scan.matches, ) tables = outdir / RESULTS_DIR write_guides_tsv( @@ -351,6 +359,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() diff --git a/src/crispr_analysis_utils/guide_alignment/scan.py b/src/crispr_analysis_utils/guide_alignment/scan.py index a51b09d..27f2547 100644 --- a/src/crispr_analysis_utils/guide_alignment/scan.py +++ b/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,7 +63,9 @@ 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. """ @@ -67,6 +73,7 @@ class ScanResult: 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,6 +172,7 @@ 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): @@ -170,13 +187,18 @@ def scan_guides( result.n_no_pam += 1 result.n_sites += 1 for row in rows_of[spacers[index]]: + result.matches[row] += 1 + if pam_only and not matched: + 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 +252,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 diff --git a/src/crispr_analysis_utils/guide_alignment/summary.py b/src/crispr_analysis_utils/guide_alignment/summary.py index ecca9ab..fc962ba 100644 --- a/src/crispr_analysis_utils/guide_alignment/summary.py +++ b/src/crispr_analysis_utils/guide_alignment/summary.py @@ -139,6 +139,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 +156,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 +185,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 +206,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 +328,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 +419,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.") diff --git a/src/crispr_analysis_utils/skills/cau-guide-alignment/SKILL.md b/src/crispr_analysis_utils/skills/cau-guide-alignment/SKILL.md index 483ea98..31dfd69 100644 --- a/src/crispr_analysis_utils/skills/cau-guide-alignment/SKILL.md +++ b/src/crispr_analysis_utils/skills/cau-guide-alignment/SKILL.md @@ -57,6 +57,7 @@ the PAM, is the library designer's answer, not something to infer. | `--orientation-check` | `error` | `error`, `warn` or `off` | | `--fix-orientation` | off | Align the reverse complement of flagged entries | | `--contigs` | all | Comma-separated contig names | +| `--pam-only` | off | Write only the sites that have a PAM | | `--id-col`, `--spacer-col` | `0`, `1` | Position or name | | `--sep`, `--header` | tab, off | `tab`, `comma`, or one character | @@ -87,6 +88,30 @@ Three things trip people up: - **`guide_id` and `mapped_sequence` line up base for base.** Diff them to see whether the 5' base paired and what the real PAM is. +## Short spacers, and `--pam-only` + +A spacer's matches grow as the genome divided by 4 to its length. A 20-mer +matches a handful of places; a 9-mer matches tens of thousands, and all but +about one in sixteen of them have no PAM, so the nuclease cannot act there. +Written out, that table is mostly rows nothing can be done with. + +`--pam-only` keeps the sites whose genomic PAM matches and drops the rest. + +- **It changes what is written, not what is searched or counted.** + `n_matches` still counts every place the protospacer occurs, which under + this flag is the only record that they exist, and `qc.tsv` still reports + `n_sites` and `n_no_pam`, so two runs stay comparable. The rows written + are `n_sites - n_no_pam`. +- **A guide whose every site is dropped keeps one row**, with the site + columns empty, exactly as a guide with no match at all does. It is in the + library, and the table says so. +- **The orientation verdict is unchanged.** A guide that matched everywhere + without a PAM is still `no_pam`, not `no_match`. +- **`guide_mapping.bed` already held only these sites**, so it does not + change. +- Leave it off for a 20-mer library: the non-PAM rows are what tell you a + spacer is promiscuous, and there are few enough of them to read. + ## Orientation Both readings of every guide go into the same scan, so an entry is diff --git a/src/crispr_analysis_utils/skills/cau-guide-alignment/references/outputs.md b/src/crispr_analysis_utils/skills/cau-guide-alignment/references/outputs.md index 265a8f2..d1ad2a5 100644 --- a/src/crispr_analysis_utils/skills/cau-guide-alignment/references/outputs.md +++ b/src/crispr_analysis_utils/skills/cau-guide-alignment/references/outputs.md @@ -29,7 +29,7 @@ empty. The guide's own columns repeat on each of its rows. | `oligo` | The sequence the library ordered, without the PAM: the spacer plus any base 5' of it. | | `n_oligos` | How many entries share this spacer. Above one the counts are a sum, because one feature is all the reads can resolve. | | `site_class` | `unique`, `multi` or `no_site`: how many places this guide can cut. | -| `n_matches` | How many places the protospacer occurs in the genome. | +| `n_matches` | How many places the protospacer occurs in the genome, whether or not a PAM sits beside it. `--pam-only` does not change it. | | `n_target_sites` | How many of those have a PAM. | | `orientation` | `forward`, `reverse_complement`, `no_pam`, `no_match` or `not_checked`. | | `orientation_evidence` | For `reverse_complement` only: `distinct` or `ambiguous`. | @@ -64,7 +64,9 @@ n_matches 89 n_target_sites 1 site_class unique The protospacer occurs 89 times; one of those has a PAM. Cas9 cuts there and nowhere else. Every one of the 89 is a row in the table, 88 of them with -`pam_matches` `FALSE`. +`pam_matches` `FALSE` -- unless the run was given `--pam-only`, which writes +only the one. `n_matches` still reads 89 either way, and after a `--pam-only` +run it is the only record that the other 88 places exist. ## `guide_mapping.bed` @@ -77,7 +79,9 @@ group by name. Matches with no PAM are not here; they are in the TSV. Three columns, `section`, `name`, `value`. The `class` and `orientation` sections each sum to the number of entries. `scan` holds what the pass saw: -`n_matches`, `n_off_contig_end`, `n_no_pam`, `n_sites`. +`n_matches`, `n_off_contig_end`, `n_no_pam`, `n_sites`. These count what was +found, not what was written, so a `--pam-only` run is comparable with a full +one: it wrote `n_sites - n_no_pam` site rows. ## `run_metadata/guide_alignment_.json` diff --git a/tests/test_cli_guide_alignment.py b/tests/test_cli_guide_alignment.py index c36318d..9355aba 100644 --- a/tests/test_cli_guide_alignment.py +++ b/tests/test_cli_guide_alignment.py @@ -59,6 +59,7 @@ def test_parse_run_defaults(): assert args.orientation_check == "error" assert args.fix_orientation is False assert args.contigs is None + assert args.pam_only is False # The columns stay strings here: `run` turns a digit into a position. assert (args.id_col, args.spacer_col, args.sep, args.header) == ( "0", @@ -89,6 +90,7 @@ def test_parse_run_every_flag(): "--sep", ",", "--header", + "--pam-only", ) assert args.pam == "NRG" assert args.spacer_length == 19 @@ -96,6 +98,7 @@ def test_parse_run_every_flag(): assert args.orientation_check == "warn" assert args.fix_orientation is True assert args.contigs == "chr1, chr2" + assert args.pam_only is True assert (args.id_col, args.spacer_col, args.sep, args.header) == ( "id", "2", diff --git a/tests/test_guide_alignment_pipeline.py b/tests/test_guide_alignment_pipeline.py index 1b4491b..eff3ccf 100644 --- a/tests/test_guide_alignment_pipeline.py +++ b/tests/test_guide_alignment_pipeline.py @@ -617,3 +617,115 @@ def test_a_5prime_base_at_a_contig_edge_does_not_crash(setup, tmp_path): assert row.site_class == "unique" assert result.scan.sites[0][0].genomic_5p_base is None assert result.scan.sites[0][0].mapped_sequence == SPACER_A + "TGG" + + +def build_blunt_reference(path): + """SPACER_A three times, a PAM beside one; SPACER_C once, with none. + + This is the shape a short spacer makes at scale: it occurs far more often + than a nuclease could act on it. Every site carries TTT on its 5' side, + so no reversed reading finds a PAM there. + """ + rng = random.Random(11) + + def filler(n): + return "".join(rng.choice("ACGT") for _ in range(n)) + + chr1 = Contig() + chr1.add(filler(30) + "TTT") + chr1.add(SPACER_A + "TGG") + chr1.add(filler(30) + "TTT") + chr1.add(SPACER_A + "TTT") + chr1.add(filler(30) + "TTT") + chr1.add(SPACER_A + "CAT") + chr1.add(filler(30) + "TTT") + chr1.add(SPACER_C + "TTT") + chr1.add(filler(30)) + path.parent.mkdir(parents=True, exist_ok=True) + write_fasta(path, {"chr1": chr1.sequence}) + return path + + +@pytest.fixture +def blunt(tmp_path): + """A library whose matches are mostly places it cannot be cut.""" + return { + "reference": build_blunt_reference(tmp_path / "ref" / "genome.fa"), + "guides": write_guides( + tmp_path / "guides.tsv", + [("g_a", SPACER_A), ("g_c", SPACER_C), ("g_nt", SPACER_NT)], + ), + "outdir": tmp_path / "out", + } + + +def both_ways(blunt): + """The same library aligned with the flag and without it.""" + full = run(blunt["guides"], blunt["reference"], blunt["outdir"] / "full") + only = run( + blunt["guides"], blunt["reference"], blunt["outdir"] / "only", pam_only=True + ) + return full, only + + +def annotation(result): + return table(result.outdir / "alignments" / "guide_genome_annotation.tsv") + + +def test_pam_only_writes_the_sites_the_nuclease_could_act_at(blunt): + full, only = both_ways(blunt) + + # Two of g_a's three matches cannot be cut, and g_c's one cannot either. + assert len(annotation(full)) == 5 + assert [row["pam_matches"] for row in annotation(only)] == ["TRUE", "", ""] + # Every site row left is one the full table also holds, PAM and all. + with_pam = [row for row in annotation(full) if row["pam_matches"] == "TRUE"] + assert [row for row in annotation(only) if row["pam_matches"]] == with_pam + # A guide whose every site went keeps the empty row a guide with no site + # gets: it is in the library, and the table still says so. + assert [row["alias"] for row in annotation(only)] == ["g_a", "g_c", "g_nt"] + + +def test_a_match_without_a_pam_is_counted_even_when_it_is_not_written(blunt): + full, only = both_ways(blunt) + + # n_matches is every place the protospacer occurs, and stays so: with + # the flag it is the only record that those places exist. + assert [row["n_matches"] for row in annotation(full)[:3]] == ["3", "3", "3"] + assert {row["alias"]: row["n_matches"] for row in annotation(only)} == { + "g_a": "3", + "g_c": "1", + "g_nt": "0", + } + + +def test_every_per_guide_column_is_the_same_either_way(blunt): + full, only = both_ways(blunt) + + assert [s.__dict__ for s in full.summaries] == [s.__dict__ for s in only.summaries] + + +def test_a_guide_that_only_matches_without_a_pam_is_not_called_unmatched(blunt): + _, only = both_ways(blunt) + + # Its sites are all dropped, so reading the empty list would call it + # unmatched. It matched; there is just nowhere to cut. + verdicts = {row["alias"]: row["orientation"] for row in annotation(only)} + assert verdicts["g_c"] == "no_pam" + assert verdicts["g_nt"] == "no_match" + + +def test_the_counts_do_not_change_with_what_is_written(blunt): + full, only = both_ways(blunt) + + # qc.tsv stays comparable between runs: n_sites is what was found and + # n_no_pam is how many of those the flag left out. + assert full.scan.counts() == only.scan.counts() + assert only.scan.counts()["n_no_pam"] == 3 + + +def test_the_run_records_whether_it_kept_only_the_pam_sites(blunt): + full, only = both_ways(blunt) + + assert full.metadata["parameters"]["pam_only"] is False + assert only.metadata["parameters"]["pam_only"] is True From d8e9b90ed7c54a6a8eb9c35aa2e2015fe035de9b Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 9 Nov 2026 11:12:29 -0500 Subject: [PATCH 2/4] FIX --pam-only, which asked the core and not the column The Rust core's `matched` and the table's `pam_matches` part company over a genomic N: the column reads one at an N of the pattern as a PAM, the core does not. The flag asked the core, so on a reference that writes its PAM as NGG every site was dropped though the table calls them all TRUE. It now asks `pattern_matches`, the same question the column answers. The counters still come from the core, so `n_no_pam` can disagree with the column on such a reference; that is older than this flag. Co-Authored-By: Claude Opus 5 --- .../guide_alignment/scan.py | 8 +++++++- tests/test_guide_alignment_pipeline.py | 19 +++++++++++++++++++ 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/src/crispr_analysis_utils/guide_alignment/scan.py b/src/crispr_analysis_utils/guide_alignment/scan.py index 27f2547..bcca9f6 100644 --- a/src/crispr_analysis_utils/guide_alignment/scan.py +++ b/src/crispr_analysis_utils/guide_alignment/scan.py @@ -186,9 +186,15 @@ def scan_guides( if not matched: result.n_no_pam += 1 result.n_sites += 1 + # Ask the question the `pam_matches` column answers, not the core's + # stricter `matched`: the two part company over a genomic N, which + # the column counts as a PAM at an N of the pattern and the core + # does not. Dropping a row the column would call TRUE is the one + # thing this flag must never do. + keep = not pam_only or pattern_matches(genomic_pam, primary) for row in rows_of[spacers[index]]: result.matches[row] += 1 - if pam_only and not matched: + if not keep: continue guide = guides[row] result.sites.setdefault(row, []).append( diff --git a/tests/test_guide_alignment_pipeline.py b/tests/test_guide_alignment_pipeline.py index eff3ccf..fc4e680 100644 --- a/tests/test_guide_alignment_pipeline.py +++ b/tests/test_guide_alignment_pipeline.py @@ -729,3 +729,22 @@ def test_the_run_records_whether_it_kept_only_the_pam_sites(blunt): assert full.metadata["parameters"]["pam_only"] is False assert only.metadata["parameters"]["pam_only"] is True + + +def test_a_genomic_n_in_the_pam_is_kept_because_the_column_calls_it_a_pam(tmp_path): + """`pam_matches` reads a genomic N at an N of the pattern as a PAM. + + The Rust core's own flag does not, so a run that trusted it would drop + the very rows the table marks TRUE. + """ + chr1 = Contig() + chr1.add("TTT") + chr1.add(SPACER_A + "NGG") + chr1.add("TTT") + reference = tmp_path / "ref.fa" + write_fasta(reference, {"chr1": chr1.sequence}) + guides = write_guides(tmp_path / "guides.tsv", [("g_a", SPACER_A)]) + + kept = run(guides, reference, tmp_path / "only", pam_only=True) + (row,) = annotation(kept) + assert (row["pam_matches"], row["site_class"]) == ("TRUE", "unique") From c298e888382887f5b5a67063a37749c19d580fc7 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 9 Nov 2026 11:15:31 -0500 Subject: [PATCH 3/4] FIX n_no_pam, which disagreed with the column beside it The counter came from the Rust core's PAM verdict and pam_matches from `base_matches`, and the two part company over a genomic N: the core's `code` has no value for one, so `matches_pam` rejects it even at an N of the pattern, while `base_matches` accepts it there by design. On a reference that writes its PAM as NGG, qc.tsv reported every site as having none while the table marked them all TRUE. Both now ask `pattern_matches`, so n_sites - n_no_pam is the number of site rows a --pam-only run writes. The reversed reading's evidence still reads the core's flag: no column reports it. Co-Authored-By: Claude Opus 5 --- .../guide_alignment/scan.py | 18 +++++++++++------- tests/test_guide_alignment_pipeline.py | 14 ++++++++++++++ 2 files changed, 25 insertions(+), 7 deletions(-) diff --git a/src/crispr_analysis_utils/guide_alignment/scan.py b/src/crispr_analysis_utils/guide_alignment/scan.py index bcca9f6..b17fdb3 100644 --- a/src/crispr_analysis_utils/guide_alignment/scan.py +++ b/src/crispr_analysis_utils/guide_alignment/scan.py @@ -176,6 +176,9 @@ def scan_guides( 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 @@ -183,15 +186,16 @@ def scan_guides( 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 - # Ask the question the `pam_matches` column answers, not the core's - # stricter `matched`: the two part company over a genomic N, which - # the column counts as a PAM at an N of the pattern and the core - # does not. Dropping a row the column would call TRUE is the one - # thing this flag must never do. - keep = not pam_only or pattern_matches(genomic_pam, primary) + keep = has_pam or not pam_only for row in rows_of[spacers[index]]: result.matches[row] += 1 if not keep: diff --git a/tests/test_guide_alignment_pipeline.py b/tests/test_guide_alignment_pipeline.py index fc4e680..f76f77c 100644 --- a/tests/test_guide_alignment_pipeline.py +++ b/tests/test_guide_alignment_pipeline.py @@ -748,3 +748,17 @@ def test_a_genomic_n_in_the_pam_is_kept_because_the_column_calls_it_a_pam(tmp_pa kept = run(guides, reference, tmp_path / "only", pam_only=True) (row,) = annotation(kept) assert (row["pam_matches"], row["site_class"]) == ("TRUE", "unique") + + +def test_qc_counts_the_matches_without_a_pam_the_way_the_table_marks_them(blunt): + full, only = both_ways(blunt) + + # qc.tsv sitting beside a table that contradicts it is worse than + # either being wrong alone. + rows = annotation(full) + assert full.scan.counts()["n_no_pam"] == sum( + row["pam_matches"] == "FALSE" for row in rows + ) + assert len(annotation(only)) - 2 == ( + only.scan.counts()["n_sites"] - only.scan.counts()["n_no_pam"] + ) From f9f52009d2f6357ec4a86b23031e09a3a5c38b0b Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 9 Nov 2026 11:18:08 -0500 Subject: [PATCH 4/4] ADD --bgzip, to write the BED block-compressed and indexed guide_mapping.bed.gz and its tabix index in place of the plain BED, so a browser or tabix can seek into a file too big to load. The rows are already in coordinate order, and they go straight into the compressed file through pysam's BGZFile, so the plain text never reaches the disk. Co-Authored-By: Claude Opus 5 --- .../cli/guide_alignment.py | 9 ++++ .../guide_alignment/pipeline.py | 15 +++++- .../guide_alignment/summary.py | 32 ++++++++++++- .../skills/cau-guide-alignment/SKILL.md | 6 +++ .../cau-guide-alignment/references/outputs.md | 4 ++ tests/test_cli_guide_alignment.py | 3 ++ tests/test_guide_alignment_pipeline.py | 47 +++++++++++++++++++ 7 files changed, 112 insertions(+), 4 deletions(-) diff --git a/src/crispr_analysis_utils/cli/guide_alignment.py b/src/crispr_analysis_utils/cli/guide_alignment.py index 98d1e05..d198f3d 100644 --- a/src/crispr_analysis_utils/cli/guide_alignment.py +++ b/src/crispr_analysis_utils/cli/guide_alignment.py @@ -118,6 +118,14 @@ def add_arguments(parser: argparse.ArgumentParser) -> None: "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: @@ -138,6 +146,7 @@ def run(args: argparse.Namespace) -> int: 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, diff --git a/src/crispr_analysis_utils/guide_alignment/pipeline.py b/src/crispr_analysis_utils/guide_alignment/pipeline.py index 936b71d..bd4b793 100644 --- a/src/crispr_analysis_utils/guide_alignment/pipeline.py +++ b/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/``. """ @@ -98,6 +99,7 @@ def run( 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", @@ -135,6 +137,9 @@ def run( 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`. @@ -179,6 +184,7 @@ def run( "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, @@ -279,7 +285,12 @@ def _run( 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"]) diff --git a/src/crispr_analysis_utils/guide_alignment/summary.py b/src/crispr_analysis_utils/guide_alignment/summary.py index fc962ba..b313dfd 100644 --- a/src/crispr_analysis_utils/guide_alignment/summary.py +++ b/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 @@ -463,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 @@ -474,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 ----- @@ -499,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: diff --git a/src/crispr_analysis_utils/skills/cau-guide-alignment/SKILL.md b/src/crispr_analysis_utils/skills/cau-guide-alignment/SKILL.md index 31dfd69..3479f95 100644 --- a/src/crispr_analysis_utils/skills/cau-guide-alignment/SKILL.md +++ b/src/crispr_analysis_utils/skills/cau-guide-alignment/SKILL.md @@ -58,6 +58,7 @@ the PAM, is the library designer's answer, not something to infer. | `--fix-orientation` | off | Align the reverse complement of flagged entries | | `--contigs` | all | Comma-separated contig names | | `--pam-only` | off | Write only the sites that have a PAM | +| `--bgzip` | off | Write `guide_mapping.bed.gz` and its `.tbi` | | `--id-col`, `--spacer-col` | `0`, `1` | Position or name | | `--sep`, `--header` | tab, off | `tab`, `comma`, or one character | @@ -112,6 +113,11 @@ Written out, that table is mostly rows nothing can be done with. - Leave it off for a 20-mer library: the non-PAM rows are what tell you a spacer is promiscuous, and there are few enough of them to read. +`--bgzip` goes with it. It writes `guide_mapping.bed.gz` and its tabix index +in place of the plain BED, so `tabix region` or a browser can seek into a +file too big to load. The rows are already in coordinate order, and they go +straight into the compressed file: the plain text never reaches the disk. + ## Orientation Both readings of every guide go into the same scan, so an entry is diff --git a/src/crispr_analysis_utils/skills/cau-guide-alignment/references/outputs.md b/src/crispr_analysis_utils/skills/cau-guide-alignment/references/outputs.md index d1ad2a5..fdd1acf 100644 --- a/src/crispr_analysis_utils/skills/cau-guide-alignment/references/outputs.md +++ b/src/crispr_analysis_utils/skills/cau-guide-alignment/references/outputs.md @@ -75,6 +75,10 @@ gets two rows. `start`/`end` are the protospacer, without the PAM. The name field is the `spacer`, which is the same on every row of one guide, so rows group by name. Matches with no PAM are not here; they are in the TSV. +With `--bgzip` this is `guide_mapping.bed.gz` with a `.tbi` beside it and no +plain BED, same rows in the same order. An empty index is still an index: a +library with nowhere to cut gets both files and no contigs. + ## `qc.tsv` Three columns, `section`, `name`, `value`. The `class` and `orientation` diff --git a/tests/test_cli_guide_alignment.py b/tests/test_cli_guide_alignment.py index 9355aba..5335bdc 100644 --- a/tests/test_cli_guide_alignment.py +++ b/tests/test_cli_guide_alignment.py @@ -60,6 +60,7 @@ def test_parse_run_defaults(): assert args.fix_orientation is False assert args.contigs is None assert args.pam_only is False + assert args.bgzip is False # The columns stay strings here: `run` turns a digit into a position. assert (args.id_col, args.spacer_col, args.sep, args.header) == ( "0", @@ -91,6 +92,7 @@ def test_parse_run_every_flag(): ",", "--header", "--pam-only", + "--bgzip", ) assert args.pam == "NRG" assert args.spacer_length == 19 @@ -99,6 +101,7 @@ def test_parse_run_every_flag(): assert args.fix_orientation is True assert args.contigs == "chr1, chr2" assert args.pam_only is True + assert args.bgzip is True assert (args.id_col, args.spacer_col, args.sep, args.header) == ( "id", "2", diff --git a/tests/test_guide_alignment_pipeline.py b/tests/test_guide_alignment_pipeline.py index f76f77c..70e5921 100644 --- a/tests/test_guide_alignment_pipeline.py +++ b/tests/test_guide_alignment_pipeline.py @@ -762,3 +762,50 @@ def test_qc_counts_the_matches_without_a_pam_the_way_the_table_marks_them(blunt) assert len(annotation(only)) - 2 == ( only.scan.counts()["n_sites"] - only.scan.counts()["n_no_pam"] ) + + +def test_bgzip_writes_an_indexed_bed_a_browser_can_seek_into(setup): + import pysam + + result = run(setup["guides"], setup["reference"], setup["outdir"], bgzip=True) + + tables = result.outdir / "alignments" + assert not (tables / "guide_mapping.bed").exists() + assert (tables / "guide_mapping.bed.gz").is_file() + assert (tables / "guide_mapping.bed.gz.tbi").is_file() + assert result.metadata["parameters"]["bgzip"] is True + + # The index is the point: it has to answer a region query. + chrom, start, _ = setup["sites"]["g_a"] + with pysam.TabixFile(str(tables / "guide_mapping.bed.gz")) as handle: + assert chrom in handle.contigs + (row,) = handle.fetch(chrom, start, start + 1) + assert row.split("\t")[:2] == [chrom, str(start)] + + +def test_the_indexed_bed_holds_exactly_what_the_plain_one_would(setup): + import gzip + + plain = run(setup["guides"], setup["reference"], setup["outdir"] / "plain") + zipped = run( + setup["guides"], setup["reference"], setup["outdir"] / "gz", bgzip=True + ) + + expected = (plain.outdir / "alignments" / "guide_mapping.bed").read_bytes() + with gzip.open(zipped.outdir / "alignments" / "guide_mapping.bed.gz") as handle: + assert handle.read() == expected + + +def test_a_library_with_nowhere_to_cut_still_gets_a_readable_index(tmp_path): + import pysam + + reference = build_blunt_reference(tmp_path / "ref" / "genome.fa") + guides = write_guides(tmp_path / "guides.tsv", [("g_c", SPACER_C)]) + + result = run(guides, reference, tmp_path / "out", bgzip=True) + + # No PAM anywhere, so the BED is empty. An empty index is still an index. + gz = result.outdir / "alignments" / "guide_mapping.bed.gz" + assert gz.with_suffix(".gz.tbi").is_file() + with pysam.TabixFile(str(gz)) as handle: + assert handle.contigs == []