diff --git a/CLAUDE.md b/CLAUDE.md index db9f600..6e3d5bb 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -322,8 +322,13 @@ and narrowing it would buy nothing. added G of a `GG` 21-mer and then called it a duplicate; - `mapped_sequence` is the guide's own bases that pair with the genome: the spacer, keeping a base 5' of it only where the genome holds the - same base. The coordinates stay on the declared protospacer, so a G that - pairs by chance never shifts a start; + same base, and `paired_5p` says whether it did (empty when the guide + has no such base, or the contig ended before it). **The table's + coordinates stay on the declared protospacer** -- `spacer`, + `target_site` and `guide_annotation`'s join are built on them -- + **while the BED's cover what pairs**, so a G sitting on a genomic G + makes that row one base longer. On the minus strand 5' is the higher + coordinate, so that moves `end`, not `start`; - coordinates are 0-based half-open protospacers, without the PAM or the added base, and `target_site` is a column, the boundary `b` (`end - 3` on `+`, `start + 3` on `-`). The IGVF-named columns keep these diff --git a/src/crispr_analysis_utils/guide_alignment/sites.py b/src/crispr_analysis_utils/guide_alignment/sites.py index 2985b1d..a4a03a2 100644 --- a/src/crispr_analysis_utils/guide_alignment/sites.py +++ b/src/crispr_analysis_utils/guide_alignment/sites.py @@ -60,3 +60,17 @@ class Site: mapped_sequence: str genomic_protospacer: str genomic_5p_base: str | None + + @property + def paired_5p(self) -> bool | None: + """Whether the guide's base 5' of the spacer pairs with the genome. + + None when there is nothing to ask: a guide with no 5' base, or a + match so close to a contig end that the base under it was never + read. A 5' base that pairs is recognised like any other, which is + why the BED's interval covers it and the declared protospacer does + not move. + """ + if self.genomic_5p_base is None: + return None + return len(self.mapped_sequence) > len(self.genomic_protospacer) + len(self.pam) diff --git a/src/crispr_analysis_utils/guide_alignment/summary.py b/src/crispr_analysis_utils/guide_alignment/summary.py index b313dfd..2d39ef7 100644 --- a/src/crispr_analysis_utils/guide_alignment/summary.py +++ b/src/crispr_analysis_utils/guide_alignment/summary.py @@ -294,7 +294,7 @@ def fan_out_sites( GROUPS = ( ("from the guide table, before anything is aligned", 6), ("about the whole guide, repeated on each of its rows", 5), - ("about this row's match; empty when the guide matches nowhere", 7), + ("about this row's match; empty when the guide matches nowhere", 8), ) """The column groups, in order, and how many columns each has.""" @@ -356,6 +356,12 @@ def fan_out_sites( "pam_matches": ( "TRUE when the genomic PAM fits the pattern, so Cas9 can act here." ), + "paired_5p": ( + "TRUE when the guide's base 5' of the protospacer is the base the" + " genome holds there, so it pairs and the site is one longer than" + " the spacer. Empty when the guide has no such base, or the match" + " sits too close to a contig end for it to have been read." + ), "mapped_sequence": ( "What the genome holds here, to compare against guide_id: the" " protospacer and the real PAM, keeping the 5' base only where the" @@ -366,7 +372,10 @@ def fan_out_sites( BED_COLUMNS = { "chr": "Contig of the match.", - "start": "First base of the protospacer, 0-based.", + "start": ( + "First base of what pairs, 0-based: the protospacer, and the base 5'" + " of it where that one pairs too." + ), "end": "One past its last base.", "name": "The protospacer, as spacer in the metadata table.", "score": "Always 0; BED6 wants the field.", @@ -444,7 +453,7 @@ def write_guides_tsv( ] found = sites.get(source, ()) if not found: - _write_row(handle, [*head, *[None] * 7]) + _write_row(handle, [*head, *[None] * 8]) for site in found: _write_row( handle, @@ -456,6 +465,7 @@ def write_guides_tsv( site.strand, _target_site(site), site.pam_class != NO_PAM_CLASS, + site.paired_5p, site.mapped_sequence, ], ) @@ -491,8 +501,11 @@ def write_guides_bed( ----- BED6 without a header, fields as in `BED_COLUMNS`, one row per site whose genomic PAM matches, so a spacer with two PAM-bearing sites gets a row for - each. The name field is the protospacer, which is the same on every row - of a guide that matches more than once. Sorted by contig, then position. + each. The interval covers what pairs, so it is one longer wherever the + guide's 5' base pairs too, and the table's `paired_5p` says where that + is. The name field stays the protospacer either way, which is what + groups a guide's rows when one of its sites pairs there and another does + not. Sorted by contig, then position. A match with no PAM is not somewhere Cas9 can act, so it is left out; guide_genome_annotation.tsv has them all. @@ -506,7 +519,7 @@ def write_guides_bed( continue seen.add(guide.guide_id) rows.extend( - (site.chr, site.start, site.end, guide.spacer, site.strand) + (site.chr, *_pairing_extent(site), guide.spacer, site.strand) for site in sites.get(source, ()) if site.pam_class != NO_PAM_CLASS ) @@ -532,6 +545,21 @@ def write_guides_bed( return path +def _pairing_extent(site: Site) -> tuple[int, int]: + """The site's bases that pair, as a 0-based half-open interval. + + The declared protospacer, widened by one at its 5' end where the guide's + own 5' base pairs: the U6 transcription start is not part of the design, + but a genome holding that base is recognised there like any other. On + the minus strand 5' is the higher coordinate, so it is `end` that moves. + """ + if not site.paired_5p: + return site.start, site.end + return ( + (site.start - 1, site.end) if site.strand == "+" else (site.start, site.end + 1) + ) + + def write_columns_yaml(path: str | os.PathLike[str]) -> None: """Write guide_genome_annotation_columns.yaml: what every column holds. 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 2dd5d45..39fbc44 100644 --- a/src/crispr_analysis_utils/skills/cau-guide-alignment/SKILL.md +++ b/src/crispr_analysis_utils/skills/cau-guide-alignment/SKILL.md @@ -121,6 +121,31 @@ 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. +## The base 5' of the protospacer + +A U6 construct starts its transcript with a G, so a library is often ordered +as G plus the protospacer. That G is not part of the design, and the tool +never infers it from a first base: `n_extra` is the length over +`--spacer-length`, so give the real spacer length and leave +`--add-leading-g` for a table of bare protospacers. + +- **A genome holding that base pairs with it.** The G is there for the + promoter, but where the genome has a G it is recognised like any other + base, and the site is one longer than the spacer. +- **`paired_5p` says where that happened**: `TRUE`, `FALSE`, or empty when + the guide has no such base or the match sat too close to a contig end to + read it. Three states, so a bare-protospacer library says nothing rather + than no. +- **The table's `start`/`end` never move.** They are the protospacer that + was declared, which is what `spacer`, `target_site` and every downstream + join are built on. +- **The BED's interval does.** It covers what pairs, so a row is 10 bases + where `paired_5p` is `TRUE` and 9 where it is not, and its `start` differs + from the table's by exactly one base there. On the minus strand 5' is the + higher coordinate, so it is `end` that moves. +- The cut site does not move either way: Cas9 cuts 3 bases 5' of the PAM, + at the other end of the protospacer. + ## 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 fdd1acf..9883249 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 @@ -36,6 +36,7 @@ empty. The guide's own columns repeat on each of its rows. | `chr`, `start`, `end`, `strand` | This row's match. 0-based, half-open, the protospacer alone. | | `target_site` | The boundary the cut falls on: between `target_site - 1` and `target_site`. | | `pam_matches` | `TRUE` when the genomic PAM fits the pattern. | +| `paired_5p` | `TRUE` when the guide's 5' base is the base the genome holds there, so the site is one longer than the spacer. Empty when the guide has no such base, or the contig ended before it. | | `mapped_sequence` | What the genome holds here, to compare against `guide_id`. | ### The two sequence columns @@ -75,6 +76,15 @@ 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. +**The interval covers what pairs, not what was declared.** Where the guide's +base 5' of the protospacer is the base the genome holds there, it pairs, and +the row is one longer: 10 bases for a G plus a 9-mer, 9 where the genome has +something else. So `start` here differs from the TSV's by exactly one base +wherever `paired_5p` is `TRUE`, and matches it everywhere else. On the minus +strand 5' is the higher coordinate, so `end` moves instead. The TSV's +coordinates are the declared protospacer and never move, because `spacer`, +`target_site` and every downstream join are built on them. + 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. diff --git a/tests/test_guide_alignment_pipeline.py b/tests/test_guide_alignment_pipeline.py index 70e5921..455ed19 100644 --- a/tests/test_guide_alignment_pipeline.py +++ b/tests/test_guide_alignment_pipeline.py @@ -50,6 +50,7 @@ "strand", "target_site", "pam_matches", + "paired_5p", "mapped_sequence", ) @@ -175,7 +176,7 @@ def test_run_writes_every_output(setup): assert len(guides["g_g"]) == 2 (g_nt,) = guides["g_nt"] assert g_nt["site_class"] == "no_site" - assert [g_nt[name] for name in SITE_COLUMNS] == [""] * 7 + assert [g_nt[name] for name in SITE_COLUMNS] == [""] * 8 (g_dup,) = guides["g_dup"] assert g_dup["guide_id"] == g_a["guide_id"] assert g_dup["start"] == str(a_start) @@ -330,7 +331,7 @@ def test_a_pam_running_off_a_contig_end_is_counted_not_reported(tmp_path): ) # Nothing was kept, so the guide gets its one row with no site in it. (row,) = table(result.outdir / "alignments" / "guide_genome_annotation.tsv") - assert [row[name] for name in SITE_COLUMNS] == [""] * 7 + assert [row[name] for name in SITE_COLUMNS] == [""] * 8 assert (row["n_matches"], row["n_target_sites"]) == ("0", "0") assert result.summaries[0].site_class == "no_site" assert result.summaries[0].orientation == "no_match" @@ -809,3 +810,86 @@ def test_a_library_with_nowhere_to_cut_still_gets_a_readable_index(tmp_path): assert gz.with_suffix(".gz.tbi").is_file() with pysam.TabixFile(str(gz)) as handle: assert handle.contigs == [] + + +SPACER_9 = SPACER_A[:9] + + +def write_ninemer(path, contigs): + reference = path / "ref" / "genome.fa" + reference.parent.mkdir(parents=True, exist_ok=True) + write_fasta(reference, contigs) + return reference + + +def ninemer_run(tmp_path, contigs, **kwargs): + """A G + 9-mer library, as a U6 construct orders it, against `contigs`.""" + guides = write_guides(tmp_path / "guides.tsv", [("g9", "G" + SPACER_9)]) + return run( + guides, + write_ninemer(tmp_path, contigs), + tmp_path / "out", + spacer_length=9, + orientation_check="off", + **kwargs, + ) + + +def test_a_leading_g_the_genome_also_holds_is_mapped_with_the_protospacer(tmp_path): + # The G is there for the U6 transcription start, but a genome holding one + # pairs with it like any other base, so the site is ten long, not nine. + result = ninemer_run(tmp_path, {"chr1": "TTTG" + SPACER_9 + "TGG" + "TTT"}) + + (row,) = annotation(result) + assert row["paired_5p"] == "TRUE" + # The table reports the protospacer that was declared, which does not move. + assert (row["start"], row["end"]) == ("4", "13") + bed = (result.outdir / "alignments" / "guide_mapping.bed").read_text().split("\t") + assert bed[:3] == ["chr1", "3", "13"] + + +def test_a_leading_g_the_genome_does_not_hold_leaves_the_site_nine_long(tmp_path): + result = ninemer_run(tmp_path, {"chr1": "TTTA" + SPACER_9 + "TGG" + "TTT"}) + + (row,) = annotation(result) + assert row["paired_5p"] == "FALSE" + bed = (result.outdir / "alignments" / "guide_mapping.bed").read_text().split("\t") + assert bed[:3] == ["chr1", "4", "13"] + + +def test_on_the_minus_strand_it_is_the_end_that_moves(tmp_path): + # 5' of the protospacer is the higher coordinate there. + from crispr_analysis_utils.guide_alignment.iupac import reverse_complement + + planted = reverse_complement("G" + SPACER_9 + "TGG") + result = ninemer_run(tmp_path, {"chr1": "TTT" + planted + "TTT"}) + + (row,) = annotation(result) + assert (row["strand"], row["paired_5p"]) == ("-", "TRUE") + assert (row["start"], row["end"]) == ("6", "15") + bed = (result.outdir / "alignments" / "guide_mapping.bed").read_text().split("\t") + assert bed[:3] == ["chr1", "6", "16"] + + +def test_a_base_the_contig_ended_before_is_neither_paired_nor_unpaired(tmp_path): + # The protospacer starts at base 0, so there is nothing under the G to + # compare: that is not the same as a base that failed to pair. + result = ninemer_run(tmp_path, {"chr1": SPACER_9 + "TGG" + "TTT"}) + + (row,) = annotation(result) + assert row["paired_5p"] == "" + bed = (result.outdir / "alignments" / "guide_mapping.bed").read_text().split("\t") + assert bed[:3] == ["chr1", "0", "9"] + + +def test_a_library_with_no_fifth_prime_base_says_nothing_either_way(setup): + # A bare-protospacer library has no such base, so the column is empty and + # no interval moves. + result = go(setup) + + rows = annotation(result) + assert {row["paired_5p"] for row in rows} == {""} + for row in rows: + if not row["start"]: + continue + assert int(row["end"]) - int(row["start"]) == 20 diff --git a/tests/test_guide_alignment_summary.py b/tests/test_guide_alignment_summary.py index de143a3..ea84cfd 100644 --- a/tests/test_guide_alignment_summary.py +++ b/tests/test_guide_alignment_summary.py @@ -72,6 +72,7 @@ "strand", "target_site", "pam_matches", + "paired_5p", "mapped_sequence", ] SITE_COLUMN_NAMES = COLUMN_NAMES[11:] @@ -350,12 +351,13 @@ def test_guides_tsv(tmp_path, guides, sites): # The fixture carries SPACER_A twice, as uniq and dup. + [SPACER_A, "2"] + ["unique", "2", "1", "forward", ""] - + ["chr1", "100", "120", "+", "117", "TRUE", SPACER_A + "TGG"] + # paired_5p is empty: these hand-built sites carry no 5' base. + + ["chr1", "100", "120", "+", "117", "TRUE", "", SPACER_A + "TGG"] ) # The guide's own columns repeat on its second row; only the site changes. assert lines[2].split("\t")[:11] == lines[1].split("\t")[:11] assert lines[2].split("\t")[11:] == ( - ["chr1", "500", "520", "+", "517", "FALSE", SPACER_A + "TTT"] + ["chr1", "500", "520", "+", "517", "FALSE", "", SPACER_A + "TTT"] ) # One row per entry per site, so the 11 sites of the library are 11 rows. assert len(lines) == 1 + 11 @@ -412,7 +414,7 @@ def test_a_guide_matching_nowhere_gets_one_row_with_empty_site_columns( table, bed = write_tables(tmp_path, guides, summaries, sites) (row,) = rows_by_alias(table)["nothing"] - assert [row[name] for name in SITE_COLUMN_NAMES] == [""] * 7 + assert [row[name] for name in SITE_COLUMN_NAMES] == [""] * 8 assert (row["n_matches"], row["n_target_sites"]) == ("0", "0") assert (row["site_class"], row["orientation"]) == ("no_site", "no_match") assert row["spacer"] == SPACER_NOTHING @@ -593,7 +595,16 @@ def test_coordinates_and_target_sites_from_a_hand_built_reference(tmp_path): "117", "156", ] - assert sorted(line[1] for line in bed) == ["100", "153", "201", "235"] + # 'kept_g' starts at 235 and its added G sits on the genomic G at 234, + # so the BED covers 21 bases there and 20 everywhere else: the declared + # protospacer in the table does not move, but what pairs is what maps. + assert sorted(line[1] for line in bed) == ["100", "153", "201", "234"] + widths = {line[3]: int(line[2]) - int(line[1]) for line in bed} + assert widths[SPACER_T] == 21 + assert {widths[s] for s in (SPACER_G, SPACER_C, SPACER_A)} == {20} + assert {row["paired_5p"] for row in rows["kept_g"]} == {"TRUE"} + assert {row["paired_5p"] for row in rows["plus"]} == {"FALSE"} + assert [rows["kept_g"][0]["start"], rows["kept_g"][0]["end"]] == ["235", "255"] def test_guide_id_ends_in_the_pam_pattern_and_mapped_sequence_in_the_real_pam(tmp_path):