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..d198f3d 100644 --- a/src/crispr_analysis_utils/cli/guide_alignment.py +++ b/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, diff --git a/src/crispr_analysis_utils/guide_alignment/pipeline.py b/src/crispr_analysis_utils/guide_alignment/pipeline.py index ab40c7f..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/``. """ @@ -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() diff --git a/src/crispr_analysis_utils/guide_alignment/scan.py b/src/crispr_analysis_utils/guide_alignment/scan.py index a51b09d..b17fdb3 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,9 +172,13 @@ 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 @@ -166,17 +186,29 @@ 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 + 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 diff --git a/src/crispr_analysis_utils/guide_alignment/summary.py b/src/crispr_analysis_utils/guide_alignment/summary.py index ecca9ab..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 @@ -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: 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..3479f95 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,8 @@ 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 | +| `--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 | @@ -87,6 +89,35 @@ 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. + +`--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 265a8f2..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 @@ -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` @@ -73,11 +75,17 @@ 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` 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..5335bdc 100644 --- a/tests/test_cli_guide_alignment.py +++ b/tests/test_cli_guide_alignment.py @@ -59,6 +59,8 @@ 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 + 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", @@ -89,6 +91,8 @@ def test_parse_run_every_flag(): "--sep", ",", "--header", + "--pam-only", + "--bgzip", ) assert args.pam == "NRG" assert args.spacer_length == 19 @@ -96,6 +100,8 @@ 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.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 1b4491b..70e5921 100644 --- a/tests/test_guide_alignment_pipeline.py +++ b/tests/test_guide_alignment_pipeline.py @@ -617,3 +617,195 @@ 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 + + +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") + + +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"] + ) + + +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 == []