Split the annotation into separate axes, and widen what it reports

CLAUDE.md

@@ -291,8 +291,17 @@ and narrowing it would buy nothing.
       position: MANE Select, else Ensembl canonical, else the longest CDS,
       ordered down to the id so a rebuild picks the same one;
     - overlaps are decided on the protospacer, distances and the CDS position
-      from the cut site. `location` and `promoter` are separate columns
-      because a cut in exon 1 is usually both;
+      from the target site. `location` (genic or intergenic), `feature` (exon or
+      intron) and `promoter` are three separate columns because a site in
+      exon 1 is usually all three;
+    - the store keeps every feature row of the annotated gene types and one
+      gene row for every other type, so `other_gene_types` can report a
+      lncRNA overlap without that overlap deciding where the guide sits.
+      `rows()` defaults to the annotated types for exactly that reason;
+    - a gene has as many promoters as transcript starts, so the TSS set is
+      one site per transcript and `promoter_enst` names every one a
+      protospacer sits in. A `--promoters` BED supplies transcript ids in an
+      eighth column, and without them those columns are empty;
     - the table keys back to `guide_genome_annotation.tsv` on `guide_id`
       plus the four site columns, and repeats nothing else. A guide with two
       sites needs all five to join unambiguously;
@@ -316,7 +325,7 @@ and narrowing it would buy nothing.
       same base. The coordinates stay on the declared protospacer, so a G that
       pairs by chance never shifts a start;
     - coordinates are 0-based half-open protospacers, without the PAM or the
-      added base, and `cut_site` is a column, the boundary `b` (`end - 3` on
+      added base, and `target_site` is a column, the boundary `b` (`end - 3` on
       `+`, `start + 3` on `-`). The IGVF-named columns keep these
       definitions, which differ from the lab's DC_TAP tables (the guide page
       has the comparison);

README.md

@@ -50,7 +50,7 @@ cpm = cau.counts_per_million(counts)
 | --- | --- |
 | `cau.accession` | Accession ids for SpCas9 guides, encoded from the guide's own sequence |
 | `cau.counts_per_million` | Counts-per-million normalization of count matrices |
-| `cau.guide_alignment` | Where a guide library maps: each guide's exact sites, class, coordinates and cut site (also the `cau guide-alignment` command) |
+| `cau.guide_alignment` | Where a guide library maps: each guide's exact sites, class, coordinates and target site (also the `cau guide-alignment` command) |
 | `cau.guide_annotation` | What those sites hit: exon, intron, promoter or intergenic, the gene and transcript, and the distance to the nearest TSS (also the `cau guide-annotation` command) |
 | `cau.element_mapping` | Which element of a screen each guide tests, and which guides are controls (also the `cau element-mapping` command) |
 | `cau.igvf_export` | The guide design in the IGVF per-guide metadata format, ready to submit or to feed the CRISPR pipeline (also the `cau igvf-export` command) |

docs/guide-alignment.md

@@ -74,8 +74,8 @@ random scaffolds are read too, and a guide will match its own alt copy.
 0-based and half-open, on the reference's forward strand, and the protospacer
 alone — the PAM and any base 5' of it are outside `[start, end)`.
 
-`cut_site` is the boundary the blunt SpCas9 cut falls on, three bases from the
-PAM: the cut lies between `cut_site - 1` and `cut_site`. On `+` that is
+`target_site` is the boundary the blunt SpCas9 cut falls on, three bases from the
+PAM: the cut lies between `target_site - 1` and `target_site`. On `+` that is
 `end - 3`, on `-` it is `start + 3`.
 
 ## Orientation

docs/guide-annotation.md

@@ -67,23 +67,37 @@ the position along the CDS: **MANE Select, else Ensembl canonical, else the
 longest CDS**, and the order is deterministic down to the transcript id, so a
 rebuild chooses the same one.
 
-`location` is called against that transcript. `exonic_any_transcript` records
+`feature` is called against that transcript. `exonic_any_transcript` records
 whether any transcript of the gene has an exon there, so a row called `intron`
 that is another transcript's exon is still visible.
 
-## Location and promoter are separate columns
+## Three axes, not one
 
-A cut 300 bp into exon 1 is both exonic and inside a promoter. A single
-column would have to pick one and hide the other, so `location` is
-`exon`/`intron`/`intergenic` and `promoter` is its own TRUE/FALSE.
+A target site 300 bp into exon 1 is in a gene, in an exon and inside a
+promoter. One
+column would have to pick one and hide the others, so there are three:
+`location` is `genic`/`intergenic`, `feature` is `exon`/`intron` within the
+reference transcript, and `promoter` is its own TRUE/FALSE.
 
 With `--gene-types` at its default, `intergenic` means "not in a
-protein-coding gene", not "in no gene at all".
+protein-coding gene", not "in no gene at all". A protospacer inside a lncRNA
+is `intergenic`, and `other_gene_types` and `other_genes` name what it
+overlaps, so the ambiguity is on the row without it deciding the call.
+
+## Every transcript has a promoter
+
+A gene has as many promoters as it has distinct transcript starts, and a
+guide at an alternative promoter is invisible to a one-TSS-per-gene file.
+Without `--promoters` the store supplies a site for every transcript, and
+`promoter_enst`, `promoter_ensg` and `promoter_gene` name the ones whose
+window the protospacer overlaps. `nearest_tss_gene` and `nearest_tss_ensg`
+stay singular, because the nearest is one.
 
 ## What is measured from where
 
 Overlaps are decided on the **protospacer**, so a guide reaching into an exon
-counts as exonic. Distances and the CDS position are measured from the **cut
+counts as exonic. Distances and the CDS position are measured from the
+**target
 site**. Every interval test ignores strand: a guide on one strand targets a
 gene on either.
 

src/crispr_analysis_utils/cli/guide_annotation.py

@@ -8,6 +8,7 @@
 
 import argparse
 
+from ..guide_annotation.annotate import DISTANCE_FROM, MIDPOINT
 from . import CommandError
 
 
@@ -75,6 +76,17 @@ def add_arguments(parser: argparse.ArgumentParser) -> None:
         default=1000,
         help="How far either side of a TSS counts as its promoter (default: 1000).",
     )
+    parser.add_argument(
+        "--distance-from",
+        choices=DISTANCE_FROM,
+        default=MIDPOINT,
+        help=(
+            "Where distance_to_tss and distance_to_gene_start are measured"
+            " from. The protospacer midpoint by default, since dCas9 occupies"
+            " the whole protospacer; target_site is the boundary 3 bases 5'"
+            " of the PAM, where a nuclease cuts."
+        ),
+    )
 
 
 def run(args: argparse.Namespace) -> int:
@@ -111,7 +123,7 @@ def _annotate(args: argparse.Namespace) -> int:
     """Annotate the aligned sites."""
     from ..guide_annotation import AnnotationError, PromoterError, StoreError
     from ..guide_annotation import run as annotate
-    from ..guide_annotation.annotate import LOCATIONS
+    from ..guide_annotation.annotate import FEATURES, GENIC, INTERGENIC
 
     try:
         result = annotate(
@@ -120,6 +132,7 @@ def _annotate(args: argparse.Namespace) -> int:
             args.outdir,
             promoters=args.promoters,
             promoter_window=args.promoter_window,
+            distance_from=args.distance_from,
         )
     except (
         FileNotFoundError,
@@ -132,9 +145,13 @@ def _annotate(args: argparse.Namespace) -> int:
         raise CommandError(str(error)) from error
 
     found = result.metadata["counts"]
+    # Both axes, because one no longer implies the other: genic says a gene is
+    # there, exon and intron say where inside it.
+    inside = ", ".join(f"{found['feature'][name]} {name}" for name in FEATURES)
     print(
         f"{found['rows']} rows for {found['sites']} sites in {found['genes']} genes: "
-        + ", ".join(f"{found['location'][name]} {name}" for name in LOCATIONS)
+        + f"{found['location'][GENIC]} genic ({inside})"
+        + f", {found['location'][INTERGENIC]} intergenic"
         + f", {found['promoter']} in a promoter"
     )
     print(f"Outputs in {result.outdir}")

src/crispr_analysis_utils/element_mapping.py

@@ -59,7 +59,7 @@
         " nowhere in the genome."
     ),
     "distance_to_tss": (
-        "Bases from the cut site to the nearest protein-coding TSS, signed in"
+        "Bases from the target site to the nearest protein-coding TSS, signed in"
         " that gene's orientation. How far a safe-targeting site has to be to"
         " be safe depends on the screen, so threshold this yourself rather"
         " than trusting the class alone. From guide_gene_annotation.tsv, so"

src/crispr_analysis_utils/guide_alignment/summary.py

@@ -6,9 +6,10 @@
 that duplicates an earlier one's sequence gets that entry's sites.
 
 Coordinates are 0-based and half-open on the reference's forward strand, and
-the protospacer excludes the PAM and any base 5' of it. The cut site is a
-boundary ``b``: the cut lies between ``b - 1`` and ``b``, 3 bases 5' of the PAM
-(``end - 3`` on ``+``, ``start + 3`` on ``-``). In the tables, a value that
+the protospacer excludes the PAM and any base 5' of it. The target site is a
+boundary ``b`` 3 bases 5' of the PAM (``end - 3`` on ``+``, ``start + 3`` on
+``-``). A nuclease cuts between ``b - 1`` and ``b``; a dCas9 only sits there,
+so the name says where the site is, not what happens to it. In the tables, a value that
 does not apply is an empty field, and booleans are ``TRUE`` or ``FALSE``.
 """
 
@@ -38,8 +39,8 @@
 CLASSES = (UNIQUE, MULTI, NO_SITE)
 """Every site class, in report order."""
 
-CUT_OFFSET = 3
-"""Bases between the cut and the PAM."""
+PAM_OFFSET = 3
+"""Bases between the target site and the PAM."""
 
 
 _QC_SECTIONS = ("guides", "class", "orientation")
@@ -73,8 +74,8 @@ class GuideSummary:
         The unique site's protospacer end; None unless unique.
     strand
         The unique site's strand; None unless unique.
-    cut_site
-        The unique site's cut boundary (see the module description); None
+    target_site
+        The unique site's target boundary (see the module description); None
         unless unique.
     site_class
         ``unique``, ``multi`` or ``no_site``.
@@ -92,7 +93,7 @@ class GuideSummary:
     n_matches
         Every exact match of the spacer, whether or not a PAM sits beside it.
         More than one means the spacer occurs more than once in the genome,
-        which `site_class` does not say: that counts cuttable sites only.
+        which `site_class` does not say: that counts PAM-bearing sites only.
     n_sites
         Sites with a matching PAM by spacer mismatches: ``n_sites[m]`` sites
         have ``m``, for ``m`` from 0 to `max_mismatches`.
@@ -110,7 +111,7 @@ class GuideSummary:
     guide_start: int | None
     guide_end: int | None
     strand: str | None
-    cut_site: int | None
+    target_site: int | None
     site_class: str
     orientation: str
     orientation_evidence: str | None
@@ -213,7 +214,7 @@ class that is neither one of the patterns nor ``none``, a guide is a
                 guide_start=unique.start if unique else None,
                 guide_end=unique.end if unique else None,
                 strand=unique.strand if unique else None,
-                cut_site=_cut_site(unique) if unique else None,
+                target_site=_target_site(unique) if unique else None,
                 site_class=site_class,
                 orientation=verdict,
                 orientation_evidence=(
@@ -287,12 +288,14 @@ def fan_out_sites(
         "The protospacer: guide_id without the PAM pattern or any base 5' of it."
     ),
     # about the whole guide
-    "site_class": "unique, multi or no_site: how many places this guide can cut.",
+    "site_class": (
+        "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."
     ),
-    "n_cut_sites": "How many of those have a PAM, so Cas9 could cut there.",
+    "n_target_sites": "How many of those have a PAM, so Cas9 can act there.",
     "orientation": (
         "Whether the sequence reads 5' to 3': forward, reverse_complement,"
         " no_pam, no_match, or not_checked."
@@ -307,11 +310,12 @@ def fan_out_sites(
     "start": "First base of the protospacer, 0-based.",
     "end": "One past its last base, so the protospacer is [start, end).",
     "strand": "Strand the protospacer reads 5' to 3' on.",
-    "cut_site": (
-        "Where Cas9 cuts: between cut_site - 1 and cut_site, 3 bases from the PAM."
+    "target_site": (
+        "Where Cas9 acts: the boundary 3 bases from the PAM. A nuclease cuts"
+        " between target_site - 1 and target_site; dCas9 binds without cutting."
     ),
     "pam_matches": (
-        "TRUE when the genomic PAM fits the pattern, so this match is cuttable."
+        "TRUE when the genomic PAM fits the pattern, so Cas9 can act here."
     ),
     "mapped_sequence": (
         "What the genome holds here, to compare against guide_id: the"
@@ -378,7 +382,7 @@ 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 be cuttable at one.
+    a guide can match many places and still bear a PAM at one.
     """
     if len(summaries) != len(guides):
         raise ValueError(f"Got {len(summaries)} summaries for {len(guides)} guides.")
@@ -407,7 +411,7 @@ def write_guides_tsv(
                         site.start,
                         site.end,
                         site.strand,
-                        _cut_site(site),
+                        _target_site(site),
                         site.pam_class != NO_PAM_CLASS,
                         site.mapped_sequence,
                     ],
@@ -419,7 +423,7 @@ def write_guides_bed(
     guides: Sequence[Guide],
     sites: Mapping[int, Sequence[Site]],
 ) -> None:
-    """Write guide_mapping.bed: every site a guide could be cut at.
+    """Write guide_mapping.bed: every site a guide could act at.
 
     Parameters
     ----------
@@ -433,11 +437,11 @@ def write_guides_bed(
     Notes
     -----
     BED6 without a header, fields as in `BED_COLUMNS`, one row per site whose
-    genomic PAM matches, so a spacer that cuts in two places gets a row for
+    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.
     A match with no PAM
-    is not somewhere anything can be cut, so it is left out;
+    is not somewhere Cas9 can act, so it is left out;
     guide_genome_annotation.tsv has them all.
     """
     # One row per distinct guide, not per entry: two entries of one sequence
@@ -564,9 +568,9 @@ def _source_rows(guides: Sequence[Guide]) -> list[int]:
     return rows
 
 
-def _cut_site(site: Site) -> int:
-    """The cut boundary of a site, `CUT_OFFSET` bases 5' of its PAM."""
-    return site.end - CUT_OFFSET if site.strand == "+" else site.start + CUT_OFFSET
+def _target_site(site: Site) -> int:
+    """The target boundary of a site, `PAM_OFFSET` bases 5' of its PAM."""
+    return site.end - PAM_OFFSET if site.strand == "+" else site.start + PAM_OFFSET
 
 
 def _open_table(path: str | os.PathLike[str], columns: Sequence[str] | None) -> TextIO: