Track merged guides, fix the measured-gene lookup, and carry the annotation tables

src/crispr_analysis_utils/cli/mudata_builder.py

@@ -56,6 +56,17 @@ def add_arguments(parser: argparse.ArgumentParser) -> None:
         ),
     )
     run_parser.add_argument("--genome-accession", help="The genome's IGVF portal id.")
+    run_parser.add_argument(
+        "--guide-annotation",
+        help=(
+            "guide_gene_annotation.tsv. Carried whole in guide.uns: it is "
+            "one row per guide per site per gene, so it cannot go in var."
+        ),
+    )
+    run_parser.add_argument(
+        "--element-mapping",
+        help="guide_element_mapping.tsv. Carried whole in guide.uns.",
+    )
     run_parser.add_argument(
         "--covariates",
         help=(
@@ -81,6 +92,8 @@ def run(args: argparse.Namespace) -> int:
             gtf_accession=args.gtf_accession,
             genome_accession=args.genome_accession,
             covariates=args.covariates,
+            guide_annotation=args.guide_annotation,
+            element_mapping=args.element_mapping,
         )
     except (BuildError, FileNotFoundError, OSError) as error:
         raise CommandError(str(error)) from error
@@ -97,5 +110,7 @@ def run(args: argparse.Namespace) -> int:
             f"Unannotated: {found['genes_unannotated']} genes, "
             f"{found['guides_unannotated']} guides"
         )
+    for name, rows in found.get("uns_tables", {}).items():
+        print(f"guide.uns[{name!r}]: {rows} rows")
     print(f"Outputs in {Path(result.outdir)}")
     return 0

src/crispr_analysis_utils/guide_alignment/summary.py

@@ -58,6 +58,13 @@ class GuideSummary:
         The name the guide table gives this entry.
     spacer
         The spacer, 5' to 3'.
+    oligo
+        The sequence the library ordered, without the PAM: the spacer plus
+        any base 5' of it.
+    n_oligos
+        How many library entries carry this spacer. More than one means a
+        feature reference built on the spacer merges them, so the counts are
+        a sum and are not comparable with a guide ordered once.
     mapped_sequence
         The guide's bases that pair with the genome at the unique site, as in
         `Site.mapped_sequence`; None unless the class is ``unique``.
@@ -104,6 +111,8 @@ class GuideSummary:
     guide_id: str
     alias: str
     spacer: str
+    oligo: str
+    n_oligos: int
     mapped_sequence: str | None
     site_pam: str | None
     has_perfect_site: bool
@@ -164,6 +173,10 @@ class that is neither one of the patterns nor ``none``, a guide is a
     patterns = validate_pams([pam])
     _check_max_mismatches(max_mismatches)
     primary = patterns[0]
+    # Two entries of one spacer are one feature in any reference built on the
+    # spacer, so their reads land together. Count them here so the table can
+    # say which guides carry a sum.
+    per_spacer = Counter(guide.spacer for guide in guides)
     summaries = []
     for guide, source in zip(guides, _source_rows(guides), strict=True):
         n_matches = 0
@@ -205,6 +218,8 @@ class that is neither one of the patterns nor ``none``, a guide is a
         summaries.append(
             GuideSummary(
                 guide_id=guide.guide_id,
+                oligo=guide.aligned_sequence,
+                n_oligos=per_spacer[guide.spacer],
                 alias=guide.alias,
                 spacer=guide.spacer,
                 mapped_sequence=unique.mapped_sequence if unique else None,
@@ -267,7 +282,7 @@ def fan_out_sites(
 
 
 GROUPS = (
-    ("from the guide table, before anything is aligned", 4),
+    ("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),
 )
@@ -287,6 +302,19 @@ def fan_out_sites(
     "spacer": (
         "The protospacer: guide_id without the PAM pattern or any base 5' of it."
     ),
+    "oligo": (
+        "The sequence the library ordered, without the PAM. It differs from"
+        " the spacer when the entry carries a base 5' of the protospacer,"
+        " such as a G added for U6 transcription."
+    ),
+    "n_oligos": (
+        "How many library entries share this spacer. Above one, a feature"
+        " reference built on the spacer cannot tell them apart, so the"
+        " counts are a sum: expect roughly that multiple of a guide ordered"
+        " once. Here the second entry is the same spacer with a 5' G, which"
+        " no read can be assigned to on its own, because the leader ends in"
+        " a poly-G run."
+    ),
     # about the whole guide
     "site_class": (
         "unique, multi or no_site: how many PAM-bearing places this guide matches."
@@ -393,6 +421,8 @@ def write_guides_tsv(
                 summary.alias,
                 guide_accession(summary.spacer),
                 summary.spacer,
+                summary.oligo,
+                summary.n_oligos,
                 summary.site_class,
                 summary.n_matches,
                 summary.n_sites[0],

src/crispr_analysis_utils/igvf_export.py

@@ -14,6 +14,7 @@
 
 import logging
 import os
+import re
 from collections.abc import Mapping
 from dataclasses import dataclass, field
 from datetime import UTC, datetime
@@ -30,6 +31,15 @@
 COLUMNS = {
     "guide_id": "Unique name of the guide.",
     "spacer": "The targeting sequence of the gRNA, read 5' to 3', in DNA bases.",
+    "oligo": (
+        "The sequence the library ordered, without the PAM: the spacer plus"
+        " any base 5' of it, such as a G added for U6 transcription."
+    ),
+    "n_oligos": (
+        "How many library entries share this spacer. Above one, they are one"
+        " feature in the counts and the reads cannot be told apart, so the"
+        " counts are a sum."
+    ),
     "targeting": (
         "TRUE when the guide targets a locus relevant to the experiment."
         " FALSE for a safe-targeting or non-targeting control."
@@ -212,9 +222,33 @@ def _element_kind(row, in_gene: set, in_exon: set, near: Mapping, window: int) -
     return DISTAL
 
 
+_COORDS = re.compile(r"^(?P<chrom>[^:]+):(?P<start>\d+)-(?P<end>\d+)$")
+
+
+def _measured_genes(known: Mapping[str, list[str]], element_id: str) -> list[str]:
+    """The measured genes of an element, matching on coordinates when needed.
+
+    ``chr:start-end`` is what a genome browser takes, and a browser is
+    1-based, so an element named that way starts one base after the BED row
+    it came from. A pairs table that concatenated the BED's own columns is
+    0-based instead, and the two spellings of one element differ by a single
+    base. All 26 validated enhancers of the hiPSC-to-EC library missed this
+    way. Only that one direction is tried: the reverse would not arise from
+    this convention, and would risk matching the element next door.
+    """
+    hit = known.get(element_id, [])
+    if hit:
+        return hit
+    found = _COORDS.match(str(element_id))
+    if not found:
+        return []
+    chrom, start, end = found["chrom"], int(found["start"]), int(found["end"])
+    return known.get(f"{chrom}:{start - 1}-{end}", [])
+
+
 def _genes_for(row, known: Mapping[str, list[str]]) -> list[str]:
     """Every gene associated with an element, measured ones first."""
-    measured = known.get(row.element_id, [])
+    measured = _measured_genes(known, row.element_id)
     if measured:
         return list(dict.fromkeys(measured))
     if row.genomic_element == PROMOTER and isinstance(row.nearest_ensg, str):
@@ -391,6 +425,12 @@ def build_table(
     out = pd.DataFrame(index=guides.index)
     out["guide_id"] = guides["alias"]
     out["spacer"] = guides["spacer"]
+    # Carried so the submitted table says which guides share a feature: their
+    # counts are a sum and are not comparable with a guide ordered once.
+    for column in ("oligo", "n_oligos"):
+        # An alignment table written before these existed has neither, and
+        # the column is then written empty rather than the table refused.
+        out[column] = guides[column] if column in guides.columns else pd.NA
 
     element = chosen.reindex(out.index)
     kind = element.map(info["genomic_element"])

src/crispr_analysis_utils/mudata_builder.py

@@ -278,6 +278,7 @@ def build(
     measurement_set: str,
     reference: dict | None = None,
     covariates: pd.DataFrame | None = None,
+    tables: dict[str, pd.DataFrame] | None = None,
 ):
     """Assemble the two matrices into one MuData.
 
@@ -296,6 +297,11 @@ def build(
     covariates
         One row per cell, indexed by barcode. Its columns join the root
         ``obs``. Nothing is computed; whatever is given is carried.
+    tables
+        Whole tables to carry in ``mod['guide'].uns``, by name. The gene
+        annotation is one row per guide per site per gene, so it cannot go in
+        ``var`` without collapsing it; keeping it whole loses nothing and
+        decides nothing.
 
     Returns
     -------
@@ -337,6 +343,9 @@ def build(
 
     # MuData takes the union of the two indexes and records membership in
     # obsm, so a barcode seen in one modality only is kept and visible.
+    for name, table in (tables or {}).items():
+        guide.uns[name] = table
+
     data = md.MuData({"gene": rna, "guide": guide})
     data.uns["reference"] = dict(reference or {})
     data.uns["measurement_set"] = measurement_set
@@ -382,6 +391,8 @@ def run(
     gtf_accession: str | None = None,
     genome_accession: str | None = None,
     covariates: str | os.PathLike[str] | None = None,
+    guide_annotation: str | os.PathLike[str] | None = None,
+    element_mapping: str | os.PathLike[str] | None = None,
 ) -> RunResult:
     """Assemble one channel's matrices into a MuData, and write it.
 
@@ -405,6 +416,10 @@ def run(
     covariates
         A table of per-cell covariates, indexed by barcode in its first
         column. Carried into the root ``obs`` unchanged.
+    guide_annotation, element_mapping
+        ``guide_gene_annotation.tsv`` and ``guide_element_mapping.tsv``.
+        Carried whole in ``mod['guide'].uns``, since neither is one row per
+        guide. ``var`` keeps the IGVF columns alone.
 
     Returns
     -------
@@ -467,6 +482,17 @@ def run(
             state["inputs"]["covariates"] = str(covariates)
             state["inputs"]["covariates_md5"] = file_md5(covariates)
 
+        tables = {}
+        for name, path in (
+            ("annotation", guide_annotation),
+            ("element_mapping", element_mapping),
+        ):
+            if path is None:
+                continue
+            tables[name] = pd.read_csv(path, sep="\t", low_memory=False)
+            state["inputs"][name] = str(path)
+            state["inputs"][f"{name}_md5"] = file_md5(path)
+
         reference = reference_block(
             loaded,
             gtf_accession=gtf_accession,
@@ -483,6 +509,7 @@ def run(
             measurement_set=measurement_set,
             reference=reference,
             covariates=covariate_table,
+            tables=tables,
         )
         data.write(outdir / RESULTS_DIR / TABLE)
 
@@ -499,6 +526,7 @@ def run(
             "guides": int(data["guide"].n_vars),
             "genes_unannotated": int(data["gene"].var["symbol"].isna().sum()),
             "guides_unannotated": int(data["guide"].var["spacer"].isna().sum()),
+            "uns_tables": {k: int(len(v)) for k, v in tables.items()},
         }
         state["complete"] = True
         state["finished"] = _now()

src/crispr_analysis_utils/skills/cau-element-mapping/SKILL.md

@@ -73,6 +73,12 @@ in full.
 - **`element_id` comes from the BED's fourth column** when it has one, so an
   id used elsewhere keeps its spelling. Otherwise it is `chr:start-end`,
   1-based, as a locus is read aloud.
+- **That id is 1-based over a 0-based interval, and the two cannot be
+  concatenated.** `chr:start-end` is what a genome browser takes, so its
+  start is the BED's `Start` plus one. A table that pastes a BED's own
+  columns back together is 0-based and will miss every element by a single
+  base: all 26 validated enhancers of one library matched nothing that way,
+  silently, until someone asked where the target genes had gone.
 - **It keys back** to `guide_genome_annotation.tsv` on `guide_id` and the
   four site columns.
 

src/crispr_analysis_utils/skills/cau-guide-alignment/SKILL.md

@@ -70,6 +70,13 @@ tables, so it always describes the folder it sits in. See
 
 Three things trip people up:
 
+- **`n_oligos` above one means those counts will be a sum.** Two library
+  entries of one spacer, typically a 20-mer and the same spacer with a 5' G,
+  are one feature in any reference built on the spacer, and no read can be
+  assigned to one or the other when the leader ends in a poly-G run: the
+  extra base is indistinguishable from the leader's own. `oligo` records what
+  was ordered, `spacer` is what the reads can resolve, and a guide flagged
+  this way is not comparable with one ordered once.
 - **`n_matches` is not `site_class`.** A protospacer can occur 89 times and
   still be `unique`, because only one of those has a PAM. `n_target_sites` is
   the one that matters for cutting; `n_matches` is the one that matters for

src/crispr_analysis_utils/skills/cau-guide-alignment/references/outputs.md

@@ -26,6 +26,8 @@ empty. The guide's own columns repeat on each of its rows.
 | `guide_id` | The sequence that was aligned, ending in the PAM pattern. Identifies the guide. |
 | `alias` | The name the input table gives it. |
 | `spacer` | The protospacer: `guide_id` without any base 5' of it or the PAM. |
+| `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_target_sites` | How many of those have a PAM. |

src/crispr_analysis_utils/skills/cau-igvf-export/references/outputs.md

@@ -18,6 +18,8 @@ guide and one spacer.
 | --- | --- |
 | `guide_id` | The guide's name, from the input table's own name for it. |
 | `spacer` | The targeting sequence, 5' to 3', in DNA bases. |
+| `oligo` | The sequence ordered, without the PAM. |
+| `n_oligos` | Entries sharing this spacer; above one the counts are a sum. |
 | `targeting` | TRUE only when the guide targets a relevant locus. |
 | `type` | `targeting`, `safe-targeting`, `non-targeting` or `positive control`. |
 | `guide_chr`, `guide_start`, `guide_end` | The protospacer, 0-based and half-open, **without** the PAM. |

tests/test_guide_alignment_summary.py

@@ -59,6 +59,8 @@
     "alias",
     "accession",
     "spacer",
+    "oligo",
+    "n_oligos",
     "site_class",
     "n_matches",
     "n_target_sites",
@@ -72,7 +74,7 @@
     "pam_matches",
     "mapped_sequence",
 ]
-SITE_COLUMN_NAMES = COLUMN_NAMES[9:]
+SITE_COLUMN_NAMES = COLUMN_NAMES[11:]
 
 
 def site(spacer, start, strand="+", mismatches=0, pam_class="NGG", pam="TGG"):
@@ -345,12 +347,14 @@ def test_guides_tsv(tmp_path, guides, sites):
     assert lines[0].split("\t") == COLUMN_NAMES
     assert lines[1].split("\t") == (
         [SPACER_A + "NGG", "uniq", guide_accession(SPACER_A), SPACER_A]
+        # 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"]
     )
     # The guide's own columns repeat on its second row; only the site changes.
-    assert lines[2].split("\t")[:9] == lines[1].split("\t")[:9]
-    assert lines[2].split("\t")[9:] == (
+    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"]
     )
     # One row per entry per site, so the 11 sites of the library are 11 rows.
@@ -623,3 +627,32 @@ def test_guide_id_ends_in_the_pam_pattern_and_mapped_sequence_in_the_real_pam(tm
         SPACER_A + "GGG",
         "G" + SPACER_T + "TGG",
     ]
+
+
+def test_n_oligos_marks_the_entries_whose_counts_will_be_a_sum():
+    """Two entries of one spacer share a feature, so their reads land together.
+
+    The pair is a 20-mer and the same spacer with a 5' G. No read can be
+    assigned to one or the other: the leader ends in a poly-G run, so the
+    extra base is indistinguishable from the leader's own. The counts are a
+    sum, and n_oligos is what says so.
+    """
+    plain, with_g = SPACER_A, "G" + SPACER_A
+    guides = read_guides(pd.DataFrame([("plain", plain), ("with_g", with_g)]))
+    summaries = summarize(
+        guides,
+        {0: [site(plain, 100)], 1: [site(plain, 100)]},
+        orientation={0: "forward", 1: "forward"},
+    )
+    assert [s.spacer for s in summaries] == [plain, plain]
+    assert [s.oligo for s in summaries] == [plain, with_g]
+    assert [s.n_oligos for s in summaries] == [2, 2]
+    # One accession, because one feature is what the reads can distinguish.
+    assert len({guide_accession(s.spacer) for s in summaries}) == 1
+
+
+def test_a_guide_ordered_once_reports_one_oligo(guides, sites):
+    summaries = summarized(guides, sites)
+    singles = [s for s in summaries if s.alias not in {"uniq", "dup"}]
+    assert {s.n_oligos for s in singles} == {1}
+    assert all(s.oligo == s.spacer for s in singles)