From 7acaf3894ae9790563deadd8786316300b5fb090 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Thu, 15 Oct 2026 09:49:14 -0400 Subject: [PATCH] REWRITE the annotation into separate axes, and widen what it reports `location` answered two questions at once, so a guide could not be genic and exonic and in a promoter, which is the ordinary case for a first exon. It is now `genic`/`intergenic`, with `feature` holding `exon`/`intron` beside it. A gene has as many promoters as it has transcript starts, and the TSS set was one per gene taken from the reference transcript, so a guide at an alternative promoter read FALSE. Every transcript now contributes a site, and `promoter_enst`, `promoter_ensg` and `promoter_gene` name the ones a protospacer sits in. A `--promoters` BED supplies transcript ids in an eighth column; without them those columns are empty. The calls stay protein-coding. The store keeps one gene row for every other type so `other_gene_types` and `other_genes` can report a lncRNA overlap without it deciding where the guide sits; `rows()` defaults to the annotated types so nothing widens by accident. An old store lacks `gene_type` and is refused by name. `distance_to_gene_start` tells a canonical promoter from an alternative one: both are flagged, and only the distance from the gene's annotated 5' end separates 27 from 497. Distances now measure from the protospacer midpoint, since this is CRISPRi and dCas9 occupies the whole protospacer; `--distance-from target_site` restores the nuclease convention, and `cds_position` keeps it regardless. `cut_site` is renamed `target_site` throughout for the same reason. Co-Authored-By: Claude Opus 5 --- CLAUDE.md | 15 +- README.md | 2 +- docs/guide-alignment.md | 4 +- docs/guide-annotation.md | 28 ++- .../cli/guide_annotation.py | 21 +- src/crispr_analysis_utils/element_mapping.py | 2 +- .../guide_alignment/summary.py | 50 ++--- .../guide_annotation/annotate.py | 188 ++++++++++++++++-- .../guide_annotation/db.py | 77 ++++++- .../guide_annotation/pipeline.py | 13 +- .../guide_annotation/promoters.py | 58 +++--- .../guide_annotation/summary.py | 62 +++++- .../skills/cau-element-mapping/SKILL.md | 2 +- .../cau-element-mapping/references/outputs.md | 2 +- .../skills/cau-guide-alignment/SKILL.md | 4 +- .../references/migration.md | 6 +- .../cau-guide-alignment/references/outputs.md | 6 +- .../skills/cau-guide-annotation/SKILL.md | 37 +++- .../references/outputs.md | 8 +- tests/guide_annotation_fixtures.py | 9 +- tests/test_annotation.py | 147 +++++++++++--- tests/test_annotation_db.py | 65 +++++- tests/test_annotation_pipeline.py | 23 ++- tests/test_guide_alignment_pipeline.py | 16 +- tests/test_guide_alignment_summary.py | 50 +++-- 25 files changed, 698 insertions(+), 197 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 65d9c1a..b7e94d7 100644 --- a/CLAUDE.md +++ b/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); diff --git a/README.md b/README.md index 4870629..16ec7af 100644 --- a/README.md +++ b/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) | diff --git a/docs/guide-alignment.md b/docs/guide-alignment.md index 0fb290c..9f7be58 100644 --- a/docs/guide-alignment.md +++ b/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 diff --git a/docs/guide-annotation.md b/docs/guide-annotation.md index 2a059b6..1f943cb 100644 --- a/docs/guide-annotation.md +++ b/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. diff --git a/src/crispr_analysis_utils/cli/guide_annotation.py b/src/crispr_analysis_utils/cli/guide_annotation.py index 13b370e..3c542da 100644 --- a/src/crispr_analysis_utils/cli/guide_annotation.py +++ b/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}") diff --git a/src/crispr_analysis_utils/element_mapping.py b/src/crispr_analysis_utils/element_mapping.py index ea5ad97..79a2d39 100644 --- a/src/crispr_analysis_utils/element_mapping.py +++ b/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" diff --git a/src/crispr_analysis_utils/guide_alignment/summary.py b/src/crispr_analysis_utils/guide_alignment/summary.py index abe737f..a9c919e 100644 --- a/src/crispr_analysis_utils/guide_alignment/summary.py +++ b/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: diff --git a/src/crispr_analysis_utils/guide_annotation/annotate.py b/src/crispr_analysis_utils/guide_annotation/annotate.py index ab08914..a1af97e 100644 --- a/src/crispr_analysis_utils/guide_annotation/annotate.py +++ b/src/crispr_analysis_utils/guide_annotation/annotate.py @@ -2,7 +2,7 @@ Overlaps are decided on the protospacer, so a guide reaching into an exon counts as exonic. Distances and the position along a CDS are measured from -the cut site, which is where the break falls. +the target site, which is where the break falls. Every interval operation ignores strand: a guide on one strand targets a gene on either. @@ -16,10 +16,18 @@ from .db import Store from .promoters import PromoterError +GENIC = "genic" +INTERGENIC = "intergenic" EXON = "exon" INTRON = "intron" -INTERGENIC = "intergenic" -LOCATIONS = (EXON, INTRON, INTERGENIC) +DELIMITER = ";" +"""Separates the ids in a list-valued column. + +A semicolon rather than a comma, because a comma already separates ids inside +the IGVF format's own ``putative_target_genes``. +""" +LOCATIONS = (GENIC, INTERGENIC) +FEATURES = (EXON, INTRON) """Where a site falls, in the order qc.tsv counts them.""" CDS = "CDS" @@ -28,9 +36,20 @@ SENSE = "sense" ANTISENSE = "antisense" -SITE_COLUMNS = ("chr", "start", "end", "strand", "cut_site") +SITE_COLUMNS = ("chr", "start", "end", "strand", "target_site") """What `annotate_sites` needs of every site.""" +MIDPOINT = "midpoint" +TARGET_SITE = "target_site" +DISTANCE_FROM = (MIDPOINT, TARGET_SITE) +"""Where a distance is measured from. + +The midpoint of the protospacer by default. A dCas9 occupies the whole +protospacer rather than acting at one base, so its centre is the more honest +reference for how far a guide sits from a TSS. `target_site` measures from +the boundary 3 bases 5' of the PAM instead, which is where a nuclease cuts. +""" + _IGNORE = "ignore" """Strand behaviour for every join: a guide targets a gene on either strand.""" @@ -45,6 +64,7 @@ def annotate_sites( promoters: pd.DataFrame, *, promoter_window: int = 1000, + distance_from: str = MIDPOINT, ) -> pd.DataFrame: """Annotate every site against the genes of the store. @@ -81,6 +101,11 @@ def annotate_sites( raise ValueError( f"promoter_window must not be negative, got {promoter_window}." ) + if distance_from not in DISTANCE_FROM: + raise ValueError( + f"distance_from must be one of {', '.join(DISTANCE_FROM)}, got" + f" {distance_from!r}." + ) sites = sites.reset_index(drop=True) _check_contigs(sites, store) @@ -94,11 +119,30 @@ def annotate_sites( "_row": np.arange(len(sites)), } ) - cuts = pr.PyRanges( + # The CDS position asks where along the coding sequence Cas9 acts, so it + # stays on the target site whatever distances are measured from. + targets = pr.PyRanges( + { + "Chromosome": sites["chr"].to_numpy(), + "Start": sites["target_site"].to_numpy(), + "End": sites["target_site"].to_numpy() + 1, + "Strand": sites["strand"].to_numpy(), + "_row": np.arange(len(sites)), + } + ) + + # Where distances are measured from. The midpoint rounds down, so a + # protospacer of even length takes the base left of centre. + anchor_at = ( + (sites["start"].to_numpy() + sites["end"].to_numpy()) // 2 + if distance_from == MIDPOINT + else sites["target_site"].to_numpy() + ) + anchors = pr.PyRanges( { "Chromosome": sites["chr"].to_numpy(), - "Start": sites["cut_site"].to_numpy(), - "End": sites["cut_site"].to_numpy() + 1, + "Start": anchor_at, + "End": anchor_at + 1, "Strand": sites["strand"].to_numpy(), "_row": np.arange(len(sites)), } @@ -108,28 +152,51 @@ def annotate_sites( pairs = pairs.merge( _exon_calls(spacers, store, pr), on=["_row", "ensg"], how="left" ) - pairs = pairs.merge(_cds_position(cuts, store, pr), on=["_row", "ensg"], how="left") + pairs = pairs.merge( + _cds_position(targets, store, pr), on=["_row", "ensg"], how="left" + ) - pairs["location"] = np.where( + # Two axes, not one. Whether a gene is there at all, and where inside it: + # a guide can be genic and in a promoter and in an exon at once, which a + # single column could not say. + pairs["location"] = np.where(pairs["ensg"].isna(), INTERGENIC, GENIC) + pairs["feature"] = np.where( pairs["ensg"].isna(), - INTERGENIC, + "", np.where(pairs["enst"].notna(), EXON, INTRON), ) pairs["region"] = np.where( - pairs["location"] == EXON, + pairs["feature"] == EXON, np.where(pairs["cds_position"].notna(), CDS, UTR), "", ) pairs["exonic_any_transcript"] = pairs["exonic_any_transcript"].fillna(False) + # How far into the gene the target site falls, from the gene's own 5' + # end. Near zero is the annotated start; a large one is deeper in, which + # is what distinguishes a canonical promoter from an alternative one. + five_prime = np.where( + pairs["gene_strand"] == "+", pairs["gene_start"], pairs["gene_end"] - 1 + ) + toward = anchor_at[pairs["_row"].to_numpy()] - five_prime + pairs["distance_to_gene_start"] = np.where( + pairs["gene_strand"].isna(), + np.nan, + np.where(pairs["gene_strand"] == "+", toward, -toward), + ) + pairs["sense"] = np.where( pairs["gene_strand"].isna(), "", np.where(pairs["strand_site"] == pairs["gene_strand"], SENSE, ANTISENSE), ) - tss = _tss_columns(cuts, spacers, promoters, promoter_window, pr) - frame = sites.join(pairs.set_index("_row").drop(columns=["strand_site"])) + tss = _tss_columns(anchors, spacers, promoters, promoter_window, pr) + others = _other_gene_overlaps(spacers, store, pr, len(sites)) + frame = sites.join( + pairs.set_index("_row").drop(columns=["strand_site", "gene_start", "gene_end"]) + ) frame = frame.merge(tss, left_index=True, right_index=True, how="left") + frame = frame.merge(others, left_index=True, right_index=True, how="left") frame["gc_percent"] = _gc_percent(frame) frame["max_polyt_run"] = _max_polyt(frame) return frame.reset_index(drop=True) @@ -170,11 +237,69 @@ def _gene_pairs(spacers, store: Store, pr) -> pd.DataFrame: "ensg": _joined(frame, "ensg", "_gene"), "gene_symbol": _joined(frame, "gene_symbol", "_gene"), "gene_strand": _joined(frame, "Strand", "_gene"), + "gene_start": _joined(frame, "Start", "_gene"), + "gene_end": _joined(frame, "End", "_gene"), "strand_site": frame["Strand"], } ) +def _other_gene_overlaps(spacers, store: Store, pr, n_rows: int) -> pd.DataFrame: + """Per site: the genes of other types the protospacer overlaps. + + Site-level, not per gene: an overlap with a lncRNA says the guide sits in + ambiguous territory whichever protein-coding gene the row is about. + """ + empty = pd.DataFrame( + { + "other_gene_types": pd.Series([""] * n_rows, dtype="string"), + "other_genes": pd.Series([""] * n_rows, dtype="string"), + }, + index=range(n_rows), + ) + others = store.other_genes() + if others.empty: + return empty + + hits = spacers.join_overlaps( + pr.PyRanges( + others[ + ["Chromosome", "Start", "End", "Strand", "ensg", "gene_type"] + ].reset_index(drop=True) + ), + join_type="inner", + strand_behavior=_IGNORE, + suffix="_other", + ) + frame = pd.DataFrame(hits) + if frame.empty: + return empty + frame = pd.DataFrame( + { + "_row": frame["_row"], + "ensg": _joined(frame, "ensg", "_other"), + "gene_type": _joined(frame, "gene_type", "_other"), + } + ) + grouped = frame.groupby("_row").agg( + other_gene_types=("gene_type", _joined_uniques), + other_genes=("ensg", _joined_uniques), + ) + rows = range(n_rows) + return pd.DataFrame( + { + name: grouped[name].reindex(rows).fillna("").astype("string") + for name in ("other_gene_types", "other_genes") + }, + index=rows, + ) + + +def _joined_uniques(values: pd.Series) -> str: + """The distinct values, sorted, as one delimited string.""" + return DELIMITER.join(sorted({str(v) for v in values if pd.notna(v) and str(v)})) + + def _exon_calls(spacers, store: Store, pr) -> pd.DataFrame: """Per site and gene: the reference exon hit, and whether any transcript has one.""" reference = store.rows("exon", reference_only=True) @@ -214,14 +339,14 @@ def _exon_calls(spacers, store: Store, pr) -> pd.DataFrame: return best.merge(anywhere, on=["_row", "ensg"], how="outer") -def _cds_position(cuts, store: Store, pr) -> pd.DataFrame: - """How far along the reference CDS the cut falls, and the CDS length.""" +def _cds_position(targets, store: Store, pr) -> pd.DataFrame: + """How far along the reference CDS the target site falls, and its length.""" blocks = _cds_blocks(store) if blocks.empty: return pd.DataFrame(columns=["_row", "ensg", "cds_position", "cds_fraction"]) columns = ["Chromosome", "Start", "End", "Strand", "ensg", "_offset", "_total"] hits = pd.DataFrame( - cuts.join_overlaps( + targets.join_overlaps( pr.PyRanges(blocks[columns].reset_index(drop=True)), join_type="inner", strand_behavior=_IGNORE, @@ -265,7 +390,7 @@ def _cds_blocks(store: Store) -> pd.DataFrame: def _tss_columns( - cuts, spacers, promoters: pd.DataFrame, window: int, pr + anchors, spacers, promoters: pd.DataFrame, window: int, pr ) -> pd.DataFrame: """Promoter flag, and the signed distance to the nearest protein-coding TSS.""" if promoters.empty: @@ -276,12 +401,21 @@ def _tss_columns( ["Chromosome", "Start", "End", "Strand", "tss_gene", "tss_ensg"] ].reset_index(drop=True) ) + # The windows carry their own ids, so an overlap says which transcripts' + # promoters a guide sits in, not merely that it sits in one. windows = pr.PyRanges( { "Chromosome": promoters["Chromosome"].to_numpy(), "Start": np.maximum(promoters["Start"].to_numpy() - window, 0), "End": promoters["End"].to_numpy() + window, "Strand": promoters["Strand"].to_numpy(), + "w_gene": promoters["gene_symbol"].to_numpy(), + "w_ensg": promoters["ensg"].to_numpy(), + "w_enst": ( + promoters["enst"].to_numpy() + if "enst" in promoters + else np.full(len(promoters), "") + ), } ) inside = pd.DataFrame( @@ -289,12 +423,26 @@ def _tss_columns( windows, join_type="inner", strand_behavior=_IGNORE, suffix="_promoter" ) ) - promoter = pd.Series(False, index=range(len(cuts)), name="promoter") + rows = range(len(anchors)) + promoter = pd.Series(False, index=rows, name="promoter") + hits = { + name: pd.Series([""] * len(anchors), index=rows, dtype="string") + for name in ("promoter_enst", "promoter_ensg", "promoter_gene") + } if not inside.empty: promoter.iloc[np.unique(inside["_row"].to_numpy())] = True + grouped = inside.groupby("_row").agg( + promoter_enst=("w_enst", _joined_uniques), + promoter_ensg=("w_ensg", _joined_uniques), + promoter_gene=("w_gene", _joined_uniques), + ) + hits = { + name: grouped[name].reindex(rows).fillna("").astype("string") + for name in hits + } near = pd.DataFrame( - cuts.nearest_ranges( + anchors.nearest_ranges( tss, strand_behavior=_IGNORE, suffix="_tss", @@ -303,13 +451,13 @@ def _tss_columns( ) ) near = near.drop_duplicates("_row").set_index("_row") - rows = range(len(promoter)) # Signed in the gene's orientation: negative upstream of the TSS. toward = near["Start"] - _joined(near, "Start", "_tss") signed = np.where(_joined(near, "Strand", "_tss") == "+", toward, -toward) return pd.DataFrame( { "promoter": promoter.to_numpy(), + **hits, "distance_to_tss": pd.Series(signed, index=near.index).reindex(rows), "nearest_tss_gene": _joined(near, "tss_gene", "_tss").reindex(rows), "nearest_tss_ensg": _joined(near, "tss_ensg", "_tss").reindex(rows), diff --git a/src/crispr_analysis_utils/guide_annotation/db.py b/src/crispr_analysis_utils/guide_annotation/db.py index e910a66..44c32d6 100644 --- a/src/crispr_analysis_utils/guide_annotation/db.py +++ b/src/crispr_analysis_utils/guide_annotation/db.py @@ -19,7 +19,12 @@ """The feature rows the store keeps. Everything else is dropped.""" GENE_TYPES = ("protein_coding",) -"""The gene types kept by default.""" +"""The gene types a guide is annotated against. + +Every call of genic, exon, intron and promoter follows these. Genes of other +types are kept as one row each, so a protospacer overlapping one can be +reported without that overlap deciding where the guide sits. +""" COLUMNS = ( "Chromosome", @@ -33,6 +38,7 @@ "exon_id", "exon_number", "is_reference", + "gene_type", ) """The store's columns. The first four are what pyranges needs.""" @@ -58,13 +64,38 @@ class Store: features: pd.DataFrame provenance: dict - def rows(self, feature: str, *, reference_only: bool = False) -> pd.DataFrame: - """The rows of one feature, optionally only the reference transcript's.""" + @property + def gene_types(self) -> tuple[str, ...]: + """The gene types this store's annotation is called against.""" + return tuple(self.provenance.get("gene_types") or GENE_TYPES) + + def rows( + self, feature: str, *, reference_only: bool = False, primary: bool = True + ) -> pd.DataFrame: + """The rows of one feature. + + Parameters + ---------- + feature + One of `FEATURES`. + reference_only + Keep only the reference transcript's rows. + primary + Keep only genes of `gene_types`. On by default, so a store that + also carries other gene types cannot widen what counts as genic. + """ rows = self.features[self.features["Feature"] == feature] + if primary: + rows = rows[rows["gene_type"].isin(self.gene_types)] if reference_only: rows = rows[rows["is_reference"]] return rows + def other_genes(self) -> pd.DataFrame: + """Gene rows whose type is not one the annotation is called against.""" + genes = self.features[self.features["Feature"] == "gene"] + return genes[~genes["gene_type"].isin(self.gene_types)] + @property def contigs(self) -> set[str]: """Every contig the store holds.""" @@ -87,7 +118,10 @@ def build_store( out The parquet file to write. gene_types - The gene types to keep. + The gene types a guide is annotated against: genic, exon, intron and + promoter all follow these, and every feature row is kept for them. + Genes of every other type are kept as one row each, so an overlap can + be reported without it deciding where the guide sits. Returns ------- @@ -108,14 +142,20 @@ def build_store( # duplicate_attr is not optional: a GENCODE row carries several `tag` # keys, and without it the reader keeps one and drops most of MANE_Select. frame = pr.read_gtf(gtf, duplicate_attr=True) + # Genes of the annotated types keep every feature row; genes of any other + # type keep one row each, enough to report an overlap without that overlap + # deciding where a guide sits. + primary = frame["gene_type"].isin(gene_types) frame = frame.loc[ - frame["gene_type"].isin(gene_types) & frame["Feature"].isin(FEATURES) + (primary & frame["Feature"].isin(FEATURES)) + | (~primary & (frame["Feature"] == "gene")) ] - if frame.empty: + kept_primary = frame["gene_type"].isin(gene_types) + if not kept_primary.any(): raise StoreError( f"{gtf} holds no {' or '.join(gene_types)} gene. Is it the right GTF?" ) - if not (frame["Feature"] == "exon").any(): + if not ((frame["Feature"] == "exon") & kept_primary).any(): raise StoreError( f"{gtf} holds no exon rows, so it cannot say exon, intron or" " transcript. A genes-only GTF is not enough; use the full" @@ -141,11 +181,23 @@ def build_store( provenance = { "gtf": str(gtf), "gene_types": list(gene_types), + "other_gene_types": sorted( + frame["gene_type"][~frame["gene_type"].isin(gene_types)].unique().tolist() + ), "cau_version": __version__, "counts": { feature: int((frame["Feature"] == feature).sum()) for feature in FEATURES }, - "genes": int(frame["ensg"][frame["Feature"] == "gene"].nunique()), + "genes": int( + frame["ensg"][ + (frame["Feature"] == "gene") & frame["gene_type"].isin(gene_types) + ].nunique() + ), + "other_genes": int( + frame["ensg"][ + (frame["Feature"] == "gene") & ~frame["gene_type"].isin(gene_types) + ].nunique() + ), "reference_transcripts": int(len(reference)), "mane_select": int( frame["mane_select"][frame["Feature"] == "transcript"].sum() @@ -189,9 +241,16 @@ def load_store(path: str | os.PathLike[str]) -> Store: table = pq.read_table(path) missing = set(COLUMNS) - set(table.column_names) if missing: + # gene_type arrived with the other-gene-type reporting, so a store + # built before it is simply old rather than foreign. + rebuild = ( + "Rebuild it with `cau guide-annotation build-db`." + if missing == {"gene_type"} + else "Rebuild it." + ) raise StoreError( f"{path} is missing {', '.join(sorted(missing))}, so it was not" - " written by this tool. Rebuild it." + f" written by this version of the tool. {rebuild}" ) raw = (table.schema.metadata or {}).get(_METADATA_KEY) return Store( diff --git a/src/crispr_analysis_utils/guide_annotation/pipeline.py b/src/crispr_analysis_utils/guide_annotation/pipeline.py index de645ea..f91a8ef 100644 --- a/src/crispr_analysis_utils/guide_annotation/pipeline.py +++ b/src/crispr_analysis_utils/guide_annotation/pipeline.py @@ -20,7 +20,7 @@ from .._files import file_md5, write_json from .._version import __version__ -from .annotate import AnnotationError, annotate_sites +from .annotate import MIDPOINT, AnnotationError, annotate_sites from .db import Store, load_store from .promoters import promoters_from_store, read_promoters from .summary import TABLE, counts, write_annotation_tsv, write_columns_yaml @@ -69,6 +69,7 @@ def run( *, promoters: str | os.PathLike[str] | None = None, promoter_window: int = PROMOTER_WINDOW, + distance_from: str = MIDPOINT, ) -> RunResult: """Annotate every aligned site against the genes of the store. @@ -84,10 +85,12 @@ def run( Where the table goes. Created if it does not exist; nothing is ever written outside it. promoters - A TSS BED. Without one, each gene's reference transcript start is - used instead. + A TSS BED. Without one, a site per transcript is taken from the store. promoter_window How far either side of a TSS counts as its promoter. + distance_from + Where distances are measured from: the protospacer midpoint, or the + target site 3 bases 5' of the PAM. Returns ------- @@ -110,6 +113,7 @@ def run( parameters = { "promoters": None if promoters is None else str(promoters), "promoter_window": promoter_window, + "distance_from": distance_from, } state: dict = { "step": STEP, @@ -187,6 +191,7 @@ def _run( loaded, tss, promoter_window=parameters["promoter_window"], + distance_from=parameters["distance_from"], ) with _timed(state, "outputs"): @@ -222,7 +227,7 @@ def _read_sites(alignment: str | os.PathLike[str] | pd.DataFrame) -> pd.DataFram "The alignment table holds no mapped site, so there is nothing to" " annotate. Every guide matched nowhere, or had no PAM." ) - for column in ("start", "end", "cut_site"): + for column in ("start", "end", "target_site"): frame[column] = frame[column].astype("int64") return frame.reset_index(drop=True) diff --git a/src/crispr_analysis_utils/guide_annotation/promoters.py b/src/crispr_analysis_utils/guide_annotation/promoters.py index 8a8a90f..7f8f1cd 100644 --- a/src/crispr_analysis_utils/guide_annotation/promoters.py +++ b/src/crispr_analysis_utils/guide_annotation/promoters.py @@ -21,16 +21,19 @@ def read_promoters(path: str | os.PathLike[str]) -> pd.DataFrame: Parameters ---------- path - A BED of transcription start sites, one per gene. Six columns at - least; a seventh is read as the gene id, which is what lets a TSS be - tied back to the store. A first line beginning with ``#`` is a header - and is skipped. + A BED of transcription start sites, one row per site: a gene with + several promoters needs one row each. Six columns at least; a seventh + is read as the gene id, which is what lets a TSS be tied back to the + store, and an eighth as the transcript id, which is what fills + ``promoter_enst``. A first line beginning with ``#`` is a header and + is skipped. Returns ------- pandas.DataFrame - ``Chromosome``, ``Start``, ``End``, ``Strand``, ``gene_symbol`` and - ``ensg``, with ``ensg`` empty when the file has no seventh column. + ``Chromosome``, ``Start``, ``End``, ``Strand``, ``gene_symbol``, + ``ensg`` and ``enst``, the last two empty when the file has no + seventh or eighth column. Raises ------ @@ -61,41 +64,48 @@ def read_promoters(path: str | os.PathLike[str]) -> pd.DataFrame: if frame.empty: raise PromoterError(f"{path} holds no transcription start site.") - names = list(_BED6) + ["ensg"] + [f"extra_{i}" for i in range(frame.shape[1] - 7)] + names = ( + list(_BED6) + + ["ensg", "enst"] + + [f"extra_{i}" for i in range(frame.shape[1] - 8)] + ) frame.columns = names[: frame.shape[1]] - if "ensg" not in frame: - frame["ensg"] = "" - return frame.assign(ensg=_unversioned(frame["ensg"]))[ - ["Chromosome", "Start", "End", "Strand", "gene_symbol", "ensg"] - ] + for column in ("ensg", "enst"): + if column not in frame: + frame[column] = "" + return frame.assign( + ensg=_unversioned(frame["ensg"]), enst=_unversioned(frame["enst"]) + )[["Chromosome", "Start", "End", "Strand", "gene_symbol", "ensg", "enst"]] + +def promoters_from_store(store, *, reference_only: bool = False) -> pd.DataFrame: + """Derive the TSS of every transcript from the store. -def promoters_from_store(store, *, reference_only: bool = True) -> pd.DataFrame: - """Derive one TSS per gene from the store, for when no BED is given. + 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, + so every transcript contributes a site by default. Parameters ---------- store A loaded `Store`. reference_only - Take the start of each gene's reference transcript. False takes the - start of the gene itself. + Take only each gene's reference transcript, one TSS per gene. This is + the older behaviour and loses alternative promoters. Returns ------- pandas.DataFrame - The same columns `read_promoters` returns, each TSS one base wide. + The same columns `read_promoters` returns, each TSS one base wide, + deduplicated so two transcripts sharing a start give one row per + transcript but not two identical ones. """ - rows = ( - store.rows("transcript", reference_only=True) - if reference_only - else store.rows("gene") - ) + rows = store.rows("transcript", reference_only=reference_only) if rows.empty: rows = store.rows("gene") forward = rows["Strand"] == "+" start = rows["Start"].where(forward, rows["End"] - 1) - return pd.DataFrame( + frame = pd.DataFrame( { "Chromosome": rows["Chromosome"].to_numpy(), "Start": start.to_numpy(), @@ -103,8 +113,10 @@ def promoters_from_store(store, *, reference_only: bool = True) -> pd.DataFrame: "Strand": rows["Strand"].to_numpy(), "gene_symbol": rows["gene_symbol"].to_numpy(), "ensg": rows["ensg"].to_numpy(), + "enst": rows["enst"].to_numpy() if "enst" in rows else "", } ) + return frame.drop_duplicates().reset_index(drop=True) def _unversioned(ids: pd.Series) -> pd.Series: diff --git a/src/crispr_analysis_utils/guide_annotation/summary.py b/src/crispr_analysis_utils/guide_annotation/summary.py index fa939b8..ae7d703 100644 --- a/src/crispr_analysis_utils/guide_annotation/summary.py +++ b/src/crispr_analysis_utils/guide_annotation/summary.py @@ -8,12 +8,21 @@ import pandas as pd -from .annotate import ANTISENSE, CDS, INTERGENIC, LOCATIONS, SENSE, UTR +from .annotate import ( + ANTISENSE, + CDS, + FEATURES, + INTERGENIC, + LOCATIONS, + SENSE, + UTR, +) GROUPS = ( ("which site this row is about, and the key back to the alignment table", 5), - ("the gene it falls in; empty when it falls in none", 10), - ("the nearest protein-coding promoter", 4), + ("the gene it falls in; empty when it falls in none", 12), + ("the genes of other types it overlaps", 2), + ("the promoters it sits in, and the nearest one", 7), ("the protospacer itself", 2), ) """The column groups, in order, and how many columns each has.""" @@ -41,8 +50,13 @@ " longest CDS. Empty unless the protospacer overlaps one of its exons." ), "location": ( - "exon, intron or intergenic, from whether the protospacer overlaps an" - " exon of the reference transcript. Independent of promoter." + "genic when the protospacer overlaps a gene of an annotated type," + " intergenic when it overlaps none. Independent of feature and" + " promoter: a guide can be all three at once." + ), + "feature": ( + "exon or intron, from whether the protospacer overlaps an exon of the" + " reference transcript. Empty on an intergenic row." ), "exonic_any_transcript": ( "TRUE when any transcript of this gene has an exon here, so a row" @@ -51,26 +65,47 @@ "region": "CDS or UTR for an exonic row, which is what decides a frameshift.", "exon_number": "Which exon of the reference transcript, counting from its 5' end.", "cds_position": ( - "How many coding bases precede the cut in the reference transcript," + "How many coding bases precede the target site in the reference" + " transcript," " so a small number is an early truncation." ), "cds_fraction": "cds_position as a fraction of the whole CDS, 0 to 1.", + "distance_to_gene_start": ( + "Bases from the target site to this gene's annotated 5' end, signed in" + " the gene's orientation: negative upstream of it, positive inside." + " Near zero is the canonical promoter; a large positive distance in a" + " promoter window is an alternative one." + ), "sense": ( "sense when the protospacer reads on the gene's own strand," " antisense otherwise." ), - # the nearest promoter + # the genes of other types + "other_gene_types": ( + "Semicolon-separated types of the genes the protospacer overlaps that" + " are not of an annotated type, such as lncRNA. Empty when none. A" + " guide here sits where a read cannot be assigned unambiguously." + ), + "other_genes": "Semicolon-separated Ensembl ids of those genes.", + # the promoters "promoter": ( "TRUE when the protospacer falls within the promoter window of any" " protein-coding TSS, whichever gene it belongs to." ), "distance_to_tss": ( - "Bases from the cut site (end - 3 on +, start + 3 on -) to the" + "Bases from the target site (end - 3 on +, start + 3 on -) to the" " nearest protein-coding TSS, signed in that gene's orientation:" " negative upstream, positive downstream." ), "nearest_tss_gene": "The gene whose TSS that is; not always this row's gene.", "nearest_tss_ensg": "That gene's Ensembl id.", + "promoter_enst": ( + "Semicolon-separated transcripts whose promoter window the" + " protospacer overlaps. Empty when the TSS file carries no transcript" + " id, which a BED of six or seven columns does not." + ), + "promoter_ensg": "Semicolon-separated gene ids of those transcripts.", + "promoter_gene": "Semicolon-separated gene symbols of those transcripts.", # the protospacer "gc_percent": "The protospacer's GC content, as a percentage.", "max_polyt_run": ( @@ -103,7 +138,13 @@ def write_annotation_tsv(path: str | os.PathLike[str], annotated: pd.DataFrame) frame = annotated.reindex(columns=annotation_columns()) for column in ("promoter", "exonic_any_transcript"): frame[column] = frame[column].map({True: "TRUE", False: "FALSE"}) - for column in ("exon_number", "cds_position", "distance_to_tss", "max_polyt_run"): + for column in ( + "exon_number", + "cds_position", + "distance_to_tss", + "distance_to_gene_start", + "max_polyt_run", + ): frame[column] = _whole_numbers(frame[column]) frame.to_csv(path, sep="\t", index=False, na_rep="") @@ -159,6 +200,9 @@ def counts(annotated: pd.DataFrame) -> dict: else 0, "genes": int(annotated["ensg"].nunique()), "location": {name: int(location.get(name, 0)) for name in LOCATIONS}, + "feature": { + name: int(Counter(annotated["feature"]).get(name, 0)) for name in FEATURES + }, "region": {name: int(region.get(name, 0)) for name in (CDS, UTR)}, "promoter": int(annotated["promoter"].sum()), "sense": {name: int(sense.get(name, 0)) for name in (SENSE, ANTISENSE)}, diff --git a/src/crispr_analysis_utils/skills/cau-element-mapping/SKILL.md b/src/crispr_analysis_utils/skills/cau-element-mapping/SKILL.md index 2c8da0a..6a8d592 100644 --- a/src/crispr_analysis_utils/skills/cau-element-mapping/SKILL.md +++ b/src/crispr_analysis_utils/skills/cau-element-mapping/SKILL.md @@ -67,7 +67,7 @@ in full. single site can appear against two of them. - **Any overlap counts, down to one base.** For CRISPRi that is right: the dCas9-KRAB footprint covers the whole protospacer, so a guide hanging off - an element's edge still represses it. For a cutting screen the cut site + an element's edge still represses it. For a cutting screen the target site would be the better test, and roughly 4% of rows are partial overlaps, so say which kind of screen it is before reading the counts. - **`element_id` comes from the BED's fourth column** when it has one, so an diff --git a/src/crispr_analysis_utils/skills/cau-element-mapping/references/outputs.md b/src/crispr_analysis_utils/skills/cau-element-mapping/references/outputs.md index bb459f9..6c0d97b 100644 --- a/src/crispr_analysis_utils/skills/cau-element-mapping/references/outputs.md +++ b/src/crispr_analysis_utils/skills/cau-element-mapping/references/outputs.md @@ -23,7 +23,7 @@ site and no element. | `element_id` | The element this site overlaps. Empty unless targeting. | | `element_category` | The elements BED's fifth column, when it has one. | | `target_class` | `targeting`, `safe_targeting` or `non_targeting`. | -| `distance_to_tss` | Signed bases from the cut site to the nearest protein-coding TSS, carried from the annotation. | +| `distance_to_tss` | Signed bases from the target site to the nearest protein-coding TSS, carried from the annotation. | ### The classes are per site 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 66a1e1f..88310d4 100644 --- a/src/crispr_analysis_utils/skills/cau-guide-alignment/SKILL.md +++ b/src/crispr_analysis_utils/skills/cau-guide-alignment/SKILL.md @@ -5,7 +5,7 @@ description: >- (the `cau guide-alignment` command, or `crispr_analysis_utils.guide_alignment.run`). Use when a user wants each guide's genomic site: which guides are unique, multi-mapping or absent; - their protospacer coordinates and cut sites, for sceptre, element annotation + their protospacer coordinates and target sites, for sceptre, element annotation or IGVF-style guide metadata; whether a guide table is reverse-complemented; or whether non-targeting controls really match nothing. Also for reading guide_genome_annotation.tsv, guide_mapping.bed, the columns file or @@ -71,7 +71,7 @@ tables, so it always describes the folder it sits in. See Three things trip people up: - **`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_cut_sites` is + 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 specificity. - **`no_match` is the right answer for a non-targeting control.** They are diff --git a/src/crispr_analysis_utils/skills/cau-guide-alignment/references/migration.md b/src/crispr_analysis_utils/skills/cau-guide-alignment/references/migration.md index 088f2b9..747ff97 100644 --- a/src/crispr_analysis_utils/skills/cau-guide-alignment/references/migration.md +++ b/src/crispr_analysis_utils/skills/cau-guide-alignment/references/migration.md @@ -28,7 +28,7 @@ coordinates at contig ends — cannot arise. | --- | --- | | `guides.tsv` and `sites.tsv` | one `guide_genome_annotation.tsv` | | `guides.bed` | `guide_mapping.bed` | -| `cut_sites.bed` | gone; `cut_site` is a column | +| `target_sites.bed` | gone; `target_site` is a column | | `orientation.tsv` | gone; its verdict and site are columns | | `summary.tsv` | `qc.tsv`, counts only | | `run.json` | `run_metadata/guide_alignment_.json` | @@ -41,9 +41,9 @@ Everything now lives under `alignments/`, `run_metadata/` and `logs/`. | --- | --- | | `guide_id` (the lab's name) | `alias` | | — | `guide_id` is now the aligned sequence | -| `targeting` | gone; use `n_cut_sites` | +| `targeting` | gone; use `n_target_sites` | | `alignment_class` | `site_class`, and `off_target_only` no longer exists | -| `n_ngg_mm` | `n_cut_sites`, exact matches only | +| `n_ngg_mm` | `n_target_sites`, exact matches only | | `min_ngg_mismatches` | gone | | `leading_g` | gone; compare the lengths of `guide_id` and `spacer` | | `pam_class` | `pam_matches`, `TRUE`/`FALSE` | 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 866b42a..4d53e72 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 @@ -28,11 +28,11 @@ empty. The guide's own columns repeat on each of its rows. | `spacer` | The protospacer: `guide_id` without any base 5' of it or the PAM. | | `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_cut_sites` | How many of those have a PAM. | +| `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`. | | `chr`, `start`, `end`, `strand` | This row's match. 0-based, half-open, the protospacer alone. | -| `cut_site` | The boundary the cut falls on: between `cut_site - 1` and `cut_site`. | +| `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. | | `mapped_sequence` | What the genome holds here, to compare against `guide_id`. | @@ -57,7 +57,7 @@ A guide can occur many times and still be cuttable once. This is the most common surprise: ```text -n_matches 89 n_cut_sites 1 site_class unique +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 diff --git a/src/crispr_analysis_utils/skills/cau-guide-annotation/SKILL.md b/src/crispr_analysis_utils/skills/cau-guide-annotation/SKILL.md index 411710b..4eea25e 100644 --- a/src/crispr_analysis_utils/skills/cau-guide-annotation/SKILL.md +++ b/src/crispr_analysis_utils/skills/cau-guide-annotation/SKILL.md @@ -4,7 +4,8 @@ description: >- Annotate aligned CRISPR guides against genes with crispr-analysis-utils (the `cau guide-annotation` command, or `crispr_analysis_utils.guide_annotation`). Use after `cau guide-alignment` when a user wants to know what each guide - hits: exon, intron or intergenic, whether it falls in a promoter, which gene + hits: genic or intergenic, exon or intron, whether it falls in a promoter + and whose, which gene and transcript, how far it is from the nearest protein-coding TSS, whether it reads sense or antisense, where in the CDS it cuts, its GC content and its poly-T run. Also for building the GENCODE annotation store, for choosing a @@ -38,7 +39,8 @@ the store in about a second, and a 43,000-site library annotates in another. - **A full annotation GTF, not a genes-only one.** The store needs `transcript` and `exon` rows. A genes-only GTF gives gene and intergenic calls and nothing else, and `build-db` refuses it. -- **Which promoter file?** Ask for the lab's TSS BED. Without `--promoters`, +- **Which promoter file?** Ask for the lab's TSS BED. A file with one site + per gene loses alternative promoters. Without `--promoters`, each gene's reference transcript start is used instead, which is a defensible default but not the same file the rest of the lab uses. - **Which gene types?** `--gene-types` defaults to `protein_coding`, so a @@ -50,8 +52,9 @@ the store in about a second, and a 43,000-site library annotates in another. | Option | Default | Meaning | | --- | --- | --- | | `--gene-types` | `protein_coding` | Which genes go in the store | -| `--promoters` | the store's TSS | A TSS BED, one site per gene | +| `--promoters` | the store's TSS | A TSS BED, one row per site. A seventh column is the gene id, an eighth the transcript id, which is what fills `promoter_enst` | | `--promoter-window` | `1000` | How far either side of a TSS is its promoter | +| `--distance-from` | `midpoint` | Where distances are measured from: the protospacer midpoint, or `target_site` 3 bases 5' of the PAM | ## Reading the result @@ -61,24 +64,42 @@ one table, `guide_gene_annotation.tsv`, and the describes it, generated from the same definitions. See [`references/outputs.md`](references/outputs.md) for the columns in full. -Four things trip people up: +Seven things trip people up: - **One row per site per gene.** A site in two overlapping genes gets two rows, so the row count is not the site count. Only the sites with a PAM are annotated, so a multi-mapper appears once per cuttable site. -- **`location` and `promoter` are separate.** A cut 300 bp into exon 1 is - both exonic and in a promoter; one column would have to hide one of them. +- **`location`, `feature` and `promoter` are three axes, not one.** + `location` is `genic` or `intergenic`, `feature` is `exon` or `intron` + inside the reference transcript, and `promoter` is independent of both: a + cut 300 bp into exon 1 is genic, exonic and in a promoter at once. - **The table keys back to the alignment table** on `guide_id` plus the four site columns. The alias, accession and spacer are not repeated. A guide with two sites needs all five columns to join unambiguously. - **`intergenic` means "not in a protein-coding gene"**, not "in no gene at - all", whenever `--gene-types` is left at its default. + all", whenever `--gene-types` is left at its default. A protospacer inside + a lncRNA is `intergenic` with `other_gene_types` naming the type and + `other_genes` its id, so the overlap is visible without deciding where the + guide sits. +- **Distances are measured from the protospacer midpoint**, not from + `target_site`. A dCas9 occupies the whole protospacer rather than acting at + one base, so its centre is the honest reference; `--distance-from + target_site` restores the nuclease convention. `cds_position` always uses + `target_site`, because it asks where along the coding sequence Cas9 acts. +- **`distance_to_gene_start` tells a canonical promoter from an alternative + one.** Both are flagged `promoter` and both name their transcripts, so the + flag alone cannot separate them; the distance from the gene's annotated 5' + end can. +- **A gene has as many promoters as transcript starts.** `promoter_enst`, + `promoter_ensg` and `promoter_gene` name every transcript whose window the + protospacer overlaps, semicolon-separated. `nearest_tss_*` stays singular, + because the nearest is one. ## The reference transcript A gene has many transcripts, so one is chosen for `enst`, `exon_number` and the CDS position: **MANE Select, else Ensembl canonical, else the longest -CDS**. `location` is called against it, and `exonic_any_transcript` says +CDS**. `feature` is called against it, and `exonic_any_transcript` says whether any other transcript has an exon there, so a row called `intron` that is someone else's exon is still visible. diff --git a/src/crispr_analysis_utils/skills/cau-guide-annotation/references/outputs.md b/src/crispr_analysis_utils/skills/cau-guide-annotation/references/outputs.md index c8d94bf..d3215e4 100644 --- a/src/crispr_analysis_utils/skills/cau-guide-annotation/references/outputs.md +++ b/src/crispr_analysis_utils/skills/cau-guide-annotation/references/outputs.md @@ -26,14 +26,18 @@ empty. Only sites with a PAM are annotated. | `chr`, `start`, `end`, `strand` | The protospacer, 0-based and half-open. | | `ensg`, `gene_symbol` | The gene this row is about, without a version. | | `enst` | Its reference transcript, when the protospacer hits one of that transcript's exons. | -| `location` | `exon`, `intron` or `intergenic`. | +| `location` | `genic` or `intergenic`, from whether a gene of an annotated type overlaps. | +| `feature` | `exon` or `intron` inside the reference transcript; empty when intergenic. | | `exonic_any_transcript` | TRUE when any transcript of the gene has an exon here. | | `region` | `CDS` or `UTR` for an exonic row. | | `exon_number` | Which exon of the reference transcript, from its 5' end. | | `cds_position`, `cds_fraction` | How far along the CDS the cut falls, in bases and as a fraction. | +| `distance_to_gene_start` | Signed bases from the anchor to this gene's annotated 5' end, in the gene's orientation. Near zero is the canonical promoter. | | `sense` | Whether the protospacer reads on the gene's own strand. | +| `other_gene_types`, `other_genes` | Semicolon-separated types and ids of overlapping genes of other types, such as lncRNA. Empty when none. | | `promoter` | TRUE within the promoter window of any protein-coding TSS. | -| `distance_to_tss` | Signed bases from the cut site to the nearest TSS. | +| `promoter_enst`, `promoter_ensg`, `promoter_gene` | Semicolon-separated transcripts whose promoter window the protospacer overlaps, and their genes. Empty when the TSS file carries no transcript id. | +| `distance_to_tss` | Signed bases from the anchor to the nearest TSS. The anchor is the protospacer midpoint unless `--distance-from target_site`. | | `nearest_tss_gene`, `nearest_tss_ensg` | Whose TSS that is. | | `gc_percent`, `max_polyt_run` | The protospacer's GC, and its longest T run. | diff --git a/tests/guide_annotation_fixtures.py b/tests/guide_annotation_fixtures.py index 4372d87..4e41628 100644 --- a/tests/guide_annotation_fixtures.py +++ b/tests/guide_annotation_fixtures.py @@ -180,7 +180,12 @@ "ENSG00000000003": 2999, "ENSG00000000005": 5000, } -"""The reference transcript's start, which is what `promoters_from_store` takes.""" +"""Each gene's reference transcript start, one per gene. + +`promoters_from_store(reference_only=True)` gives exactly these. The default +gives a site per transcript, so a gene with an alternative start contributes +more than one. +""" def _exon_ids(genes): @@ -314,7 +319,7 @@ def site(guide_id, chrom, start, strand, *, spacer=SPACER, length=20, pam_matche "start": start, "end": end, "strand": strand, - "cut_site": end - 3 if strand == "+" else start + 3, + "target_site": end - 3 if strand == "+" else start + 3, "spacer": spacer, "pam_matches": pam_matches, } diff --git a/tests/test_annotation.py b/tests/test_annotation.py index 2899e42..fca53e2 100644 --- a/tests/test_annotation.py +++ b/tests/test_annotation.py @@ -42,6 +42,12 @@ site("downstream_of_minus_tss", "chr1", 4390, "+"), # Upstream of AAA's TSS at 1000, and so in no gene but in a promoter. site("upstream_of_plus_tss", "chr1", 960, "+"), + # Inside DDD, which is lncRNA: no store annotates against it, so the row + # is intergenic and the overlap is reported on its own columns. + site("in_a_lncrna", "chr1", 7500, "+"), + # EEE's second transcript starts at 5500; its first starts at 5000, so + # this sits in the alternative promoter and no other. + site("alternative_promoter", "chr2", 5480, "+"), ) @@ -83,7 +89,7 @@ def test_the_input_columns_come_first_and_untouched(annotated): assert_frame_equal( annotated[carried].drop_duplicates(ignore_index=True), sites_frame(*SITES) ) - assert annotated["cut_site"].tolist()[:2] == [2047, 2033] + assert annotated["target_site"].tolist()[:2] == [2047, 2033] def test_one_row_per_site_except_where_two_genes_overlap(annotated): @@ -95,20 +101,21 @@ def test_one_row_per_site_except_where_two_genes_overlap(annotated): @pytest.mark.parametrize( - ("guide_id", "location"), + ("guide_id", "location", "feature"), [ - ("exon_plus", "exon"), - ("exon_minus", "exon"), - ("other_transcript_exon", "intron"), - ("deep_intron", "intron"), - ("minus_cds_first_block", "exon"), - ("utr_and_promoter", "exon"), - ("intergenic", "intergenic"), - ("upstream_of_plus_tss", "intergenic"), + ("exon_plus", "genic", "exon"), + ("exon_minus", "genic", "exon"), + ("other_transcript_exon", "genic", "intron"), + ("deep_intron", "genic", "intron"), + ("minus_cds_first_block", "genic", "exon"), + ("utr_and_promoter", "genic", "exon"), + ("intergenic", "intergenic", ""), + ("upstream_of_plus_tss", "intergenic", ""), ], ) -def test_where_each_site_falls(annotated, guide_id, location): - assert one(annotated, guide_id)["location"] == location +def test_where_each_site_falls(annotated, guide_id, location, feature): + row = one(annotated, guide_id) + assert (row["location"], row["feature"]) == (location, feature) def test_an_exonic_row_names_the_reference_transcript_and_its_exon(annotated): @@ -137,7 +144,8 @@ def test_a_site_in_two_genes_gets_a_row_for_each(annotated): rows = rows_for(annotated, "two_genes") assert rows["gene_symbol"].tolist() == ["AAA", "BBB"] # The same bases are AAA's intron and BBB's third exon. - assert rows["location"].tolist() == ["intron", "exon"] + assert rows["feature"].tolist() == ["intron", "exon"] + assert rows["location"].tolist() == ["genic", "genic"] assert rows["enst"].tolist()[1] == "ENST00000000021" assert rows["exon_number"].tolist()[1] == "3" @@ -145,10 +153,10 @@ def test_a_site_in_two_genes_gets_a_row_for_each(annotated): # --- another transcript's exon -------------------------------------------- -def test_a_non_reference_exon_is_flagged_without_changing_the_location(annotated): +def test_a_non_reference_exon_is_flagged_without_changing_the_feature(annotated): row = one(annotated, "other_transcript_exon") # A2 has an exon at 1500-1600; the reference transcript A1 has an intron. - assert row["location"] == "intron" + assert row["feature"] == "intron" assert row["exonic_any_transcript"] is True assert pd.isna(row["enst"]) assert row["region"] == "" @@ -169,14 +177,14 @@ def test_a_guide_on_the_other_strand_is_still_in_the_gene(annotated): # AAA reads on +, this guide on -. Every join ignores strand, so the site # is annotated; a same-strand join would call it intergenic. row = one(annotated, "exon_minus") - assert (row["gene_symbol"], row["location"]) == ("AAA", "exon") + assert (row["gene_symbol"], row["feature"]) == ("AAA", "exon") assert row["enst"] == "ENST00000000011" def test_a_minus_strand_gene_annotates_a_plus_strand_guide(annotated): # The other direction: BBB reads on -, this guide on +. row = rows_for(annotated, "two_genes").iloc[1] - assert (row["gene_symbol"], row["location"]) == ("BBB", "exon") + assert (row["gene_symbol"], row["feature"]) == ("BBB", "exon") @pytest.mark.parametrize( @@ -233,7 +241,7 @@ def test_the_lower_block_of_a_minus_strand_gene_carries_the_offset(annotated): def test_an_exon_whose_cut_misses_the_cds_is_utr(annotated): # The cut at 1047 is in AAA's exon 1 but 5' of the start codon at 1100. row = one(annotated, "utr_and_promoter") - assert (row["location"], row["region"], row["exon_number"]) == ("exon", "UTR", "1") + assert (row["feature"], row["region"], row["exon_number"]) == ("exon", "UTR", "1") assert pd.isna(row["cds_position"]) @@ -253,15 +261,16 @@ def test_region_is_empty_off_an_exon(annotated, guide_id): ("guide_id", "distance", "gene"), [ # AAA's TSS is at 1000 and reads on +, so a later cut is downstream. - ("utr_and_promoter", 47, "AAA"), - ("exon_plus", 1047, "AAA"), + # Measured from the protospacer midpoint, which is the default. + ("utr_and_promoter", 40, "AAA"), + ("exon_plus", 1040, "AAA"), # Before the TSS of a + gene: upstream, so negative. - ("upstream_of_plus_tss", -23, "AAA"), - # BBB's TSS is at 4499 and reads on -, so an earlier cut is downstream. - ("downstream_of_minus_tss", 92, "BBB"), - ("minus_cds_first_block", 956, "BBB"), + ("upstream_of_plus_tss", -30, "AAA"), + # BBB's TSS is at 4499 and reads on -, so an earlier site is downstream. + ("downstream_of_minus_tss", 99, "BBB"), + ("minus_cds_first_block", 949, "BBB"), # After the TSS of a - gene: upstream, so negative. - ("intergenic", -15518, "BBB"), + ("intergenic", -15511, "BBB"), ], ) def test_the_distance_to_the_tss_is_signed_in_the_genes_orientation( @@ -299,7 +308,7 @@ def test_the_promoter_flag_is_the_window_around_any_tss(annotated, guide_id, pro def test_a_site_can_be_an_exon_and_a_promoter_at_once(annotated): row = one(annotated, "utr_and_promoter") - assert row["location"] == "exon" + assert row["feature"] == "exon" assert bool(row["promoter"]) is True @@ -374,10 +383,10 @@ def test_the_contig_message_shows_both_namings(store, tss): def test_one_shared_contig_is_enough(store, tss): sites = sites_frame(site("here", "chr1", 2030, "+"), site("away", "chrZ", 10, "+")) frame = annotate_sites(sites, store, tss, promoter_window=WINDOW) - assert frame["location"].tolist() == ["exon", "intergenic"] + assert frame["location"].tolist() == ["genic", "intergenic"] -@pytest.mark.parametrize("missing", ["chr", "start", "end", "strand", "cut_site"]) +@pytest.mark.parametrize("missing", ["chr", "start", "end", "strand", "target_site"]) def test_a_missing_site_column_is_named(store, tss, missing): sites = sites_frame(site("g", "chr1", 2030, "+")).drop(columns=[missing]) with pytest.raises(AnnotationError, match=f"sites are missing .*{missing}"): @@ -392,3 +401,85 @@ def test_a_negative_promoter_window_is_refused(store, tss): def test_a_promoter_table_with_no_tss_is_refused(store, tss): with pytest.raises(PromoterError, match="holds no transcription start site"): annotate_sites(sites_frame(*SITES), store, tss.iloc[:0]) + + +def test_a_lncrna_overlap_is_reported_without_making_the_row_genic(annotated): + row = one(annotated, "in_a_lncrna") + assert (row["location"], row["feature"]) == ("intergenic", "") + assert pd.isna(row["ensg"]) + assert row["other_gene_types"] == "lncRNA" + assert row["other_genes"] == "ENSG00000000004" + + +def test_a_row_over_no_other_gene_type_leaves_those_columns_empty(annotated): + row = one(annotated, "exon_plus") + assert (row["other_gene_types"], row["other_genes"]) == ("", "") + + +def test_an_alternative_promoter_is_found_and_names_its_transcript(annotated): + row = one(annotated, "alternative_promoter") + assert row["promoter"] + # The one-TSS-per-gene file would have had only ENST00000000051 at 5000. + assert row["promoter_enst"] == "ENST00000000052" + assert row["promoter_ensg"] == "ENSG00000000005" + assert row["promoter_gene"] == "EEE" + + +def test_a_site_in_no_promoter_leaves_the_promoter_ids_empty(annotated): + row = one(annotated, "deep_intron") + assert not row["promoter"] + assert (row["promoter_enst"], row["promoter_ensg"], row["promoter_gene"]) == ( + "", + "", + "", + ) + + +def test_the_distance_to_the_gene_start_is_signed_in_the_genes_orientation(annotated): + # AAA runs + from 1000; this cut is in its exon 2, well inside. + plus = one(annotated, "exon_plus") + assert plus["distance_to_gene_start"] > 0 + + # upstream_of_plus_tss sits before AAA entirely, so it is in no gene and + # the column is empty rather than negative. + outside = one(annotated, "upstream_of_plus_tss") + assert pd.isna(outside["distance_to_gene_start"]) + + +def test_a_minus_strand_gene_measures_from_its_own_5_prime_end(annotated, store): + row = one(annotated, "minus_cds_first_block") + gene = store.rows("gene") + bbb = gene[gene["gene_symbol"] == row["gene_symbol"]].iloc[0] + assert bbb["Strand"] == "-" + # Measured from End - 1, the gene's 5'-most base, counting inwards. + midpoint = (row["start"] + row["end"]) // 2 + assert row["distance_to_gene_start"] == bbb["End"] - 1 - midpoint + + +def test_the_gene_start_distance_separates_canonical_from_alternative(annotated): + # EEE starts at 5000; the alternative promoter sits 480 bases in, which is + # what tells the two apart when both are flagged promoter. + alt = one(annotated, "alternative_promoter") + assert alt["promoter"] + # EEE starts at 5000 and the protospacer midpoint is 5490. + assert alt["distance_to_gene_start"] == 490 + + +@pytest.mark.parametrize( + ("distance_from", "expected"), + [("midpoint", 40), ("target_site", 47)], +) +def test_distance_from_chooses_the_reference_point(store, tss, distance_from, expected): + frame = annotate_sites( + sites_frame(*SITES), + store, + tss, + promoter_window=WINDOW, + distance_from=distance_from, + ) + assert one(frame, "utr_and_promoter")["distance_to_tss"] == expected + + +def test_an_unknown_distance_reference_is_refused(store, tss): + with pytest.raises(ValueError, match="distance_from must be one of"): + annotate_sites(sites_frame(*SITES), store, tss, distance_from="pam") diff --git a/tests/test_annotation_db.py b/tests/test_annotation_db.py index 4512dfc..c0ad80b 100644 --- a/tests/test_annotation_db.py +++ b/tests/test_annotation_db.py @@ -52,21 +52,38 @@ def test_the_store_holds_the_five_features_and_their_columns(store): assert set(store.features["Feature"].unique()) == set(FEATURES) -def test_only_the_asked_for_gene_types_survive(store): - # DDD is lncRNA, so the default protein_coding store knows nothing of it. - assert sorted(store.features["gene_symbol"].unique()) == [ +def test_only_the_asked_for_gene_types_are_annotated_against(store): + # DDD is lncRNA. It is kept, so an overlap with it can be reported, but + # only as a gene row: nothing of it can decide where a guide sits. + assert sorted(store.rows("gene")["gene_symbol"].unique()) == [ "AAA", "BBB", "CCC", "EEE", ] + assert sorted(store.other_genes()["gene_symbol"].unique()) == ["DDD"] assert store.provenance["genes"] == 4 + assert store.provenance["other_genes"] == 1 assert store.contigs == {"chr1", "chr2"} -def test_another_gene_type_keeps_what_the_default_drops(gtf, tmp_path): +def test_a_gene_of_another_type_keeps_no_feature_but_its_own_row(store): + ddd = store.features[store.features["gene_symbol"] == "DDD"] + assert ddd["Feature"].tolist() == ["gene"] + # So it can never reach the exon, intron or promoter calls. + for feature in ("transcript", "exon", "CDS", "UTR"): + assert "DDD" not in set(store.rows(feature, primary=False)["gene_symbol"]) + + +def test_another_gene_type_swaps_which_genes_are_annotated_against(gtf, tmp_path): store = build_store(gtf, tmp_path / "lnc.parquet", gene_types=("lncRNA",)) - assert sorted(store.features["gene_symbol"].unique()) == ["DDD"] + assert sorted(store.rows("gene")["gene_symbol"].unique()) == ["DDD"] + assert sorted(store.other_genes()["gene_symbol"].unique()) == [ + "AAA", + "BBB", + "CCC", + "EEE", + ] assert store.provenance["gene_types"] == ["lncRNA"] @@ -92,7 +109,9 @@ def test_the_counts_in_the_provenance_are_the_rows_kept(store): counts = store.provenance["counts"] for feature, count in counts.items(): assert count == int((store.features["Feature"] == feature).sum()) - assert counts == {"gene": 4, "transcript": 9, "exon": 21, "CDS": 13, "UTR": 23} + # Five genes, because the lncRNA keeps its one row; every other feature + # count is of the annotated types alone. + assert counts == {"gene": 5, "transcript": 9, "exon": 21, "CDS": 13, "UTR": 23} def test_exons_are_numbered_from_the_transcripts_own_5_prime_end(store): @@ -272,11 +291,13 @@ def test_a_six_column_bed_is_read_with_no_gene_id(tmp_path): "Strand", "gene_symbol", "ensg", + "enst", ] assert frame["Start"].tolist() == [1000, 4499] assert frame["gene_symbol"].tolist() == ["AAA", "BBB"] assert frame["Strand"].tolist() == ["+", "-"] assert frame["ensg"].tolist() == ["", ""] + assert frame["enst"].tolist() == ["", ""] def test_a_leading_hash_line_is_a_header_and_not_a_site(tmp_path): @@ -286,13 +307,25 @@ def test_a_leading_hash_line_is_a_header_and_not_a_site(tmp_path): assert_frame_equal(with_header, without) -def test_a_seventh_column_is_the_gene_id_and_an_eighth_is_ignored(tmp_path): +def test_a_seventh_column_is_the_gene_id_and_an_eighth_the_transcript(tmp_path): bed = write_bed( tmp_path / "tss8.bed", - [("chr1", "1000", "1001", "AAA", "0", "+", "ENSG00000000001.5", "note")], + [ + ( + "chr1", + "1000", + "1001", + "AAA", + "0", + "+", + "ENSG00000000001.5", + "ENST00000000011.2", + ) + ], ) frame = read_promoters(bed) assert frame["ensg"].tolist() == ["ENSG00000000001"] + assert frame["enst"].tolist() == ["ENST00000000011"] assert list(frame.columns) == [ "Chromosome", "Start", @@ -300,6 +333,7 @@ def test_a_seventh_column_is_the_gene_id_and_an_eighth_is_ignored(tmp_path): "Strand", "gene_symbol", "ensg", + "enst", ] @@ -323,8 +357,17 @@ def test_too_few_columns_names_the_six_a_tss_bed_needs(tmp_path, width): # --- promoters from the store --------------------------------------------- -def test_the_store_gives_one_tss_per_gene_at_the_reference_start(store): +def test_the_store_gives_a_tss_for_every_transcript(store): frame = promoters_from_store(store) + # EEE has a second transcript starting at 5500, and a guide at that + # promoter is invisible to a one-site-per-gene file. + assert sorted(frame[frame["ensg"] == "ENSG00000000005"]["Start"]) == [5000, 5500] + assert len(frame) == len(store.rows("transcript")) + assert (frame["End"] - frame["Start"]).unique().tolist() == [1] + + +def test_reference_only_gives_back_one_tss_per_gene(store): + frame = promoters_from_store(store, reference_only=True) assert dict(zip(frame["ensg"], frame["Start"], strict=True)) == TSS assert (frame["End"] - frame["Start"]).tolist() == [1, 1, 1, 1] @@ -347,5 +390,7 @@ def test_promoters_from_the_store_carry_the_gene_columns(store): "Strand", "gene_symbol", "ensg", + "enst", ] - assert frame["gene_symbol"].tolist() == ["AAA", "BBB", "CCC", "EEE"] + assert sorted(set(frame["gene_symbol"])) == ["AAA", "BBB", "CCC", "EEE"] + assert frame["enst"].str.startswith("ENST").all() diff --git a/tests/test_annotation_pipeline.py b/tests/test_annotation_pipeline.py index b72a2ae..8c41661 100644 --- a/tests/test_annotation_pipeline.py +++ b/tests/test_annotation_pipeline.py @@ -202,7 +202,8 @@ def test_the_record_says_the_run_finished_and_what_it_found(store, tmp_path): assert state["complete"] is True assert state["error"] is None assert state["counts"]["rows"] == len(result.annotated) - assert state["counts"]["location"] == {"exon": 2, "intron": 1, "intergenic": 1} + assert state["counts"]["location"] == {"genic": 3, "intergenic": 1} + assert state["counts"]["feature"] == {"exon": 2, "intron": 1} assert state["counts"]["sites"] == 3 assert state["counts"]["genes"] == 2 @@ -214,7 +215,11 @@ def test_the_record_carries_the_inputs_and_the_stores_provenance(store, tmp_path assert state["inputs"]["alignment"] == str(alignment) assert len(state["inputs"]["alignment_md5"]) == 32 assert state["store"]["gene_types"] == ["protein_coding"] - assert state["parameters"] == {"promoters": None, "promoter_window": 1000} + assert state["parameters"] == { + "promoters": None, + "promoter_window": 1000, + "distance_from": "midpoint", + } assert set(state["seconds"]) >= {"store", "sites", "promoters", "annotate", "total"} @@ -229,8 +234,10 @@ def test_a_store_handed_over_in_memory_is_recorded_as_such(store, tmp_path): def test_without_a_promoter_file_the_store_supplies_the_tss(store, tmp_path): run(write_alignment(tmp_path / "aln.tsv"), store, tmp_path / "out") + # One per transcript now, not one per gene. Transcripts sharing a start + # stay as separate sites, which is what lets promoter_enst name them all. assert ( - "4 TSS from the store" + "9 TSS from the store" in (tmp_path / "out" / "logs" / "guide_annotation.log").read_text() ) @@ -246,7 +253,11 @@ def test_a_promoter_file_is_read_and_recorded(store, tmp_path): promoter_window=50, ) state = record_of(tmp_path / "out") - assert state["parameters"] == {"promoters": str(bed), "promoter_window": 50} + assert state["parameters"] == { + "promoters": str(bed), + "promoter_window": 50, + "distance_from": "midpoint", + } assert len(state["inputs"]["promoters_md5"]) == 32 # Every distance is now measured from the one TSS the file gives. assert set(result.annotated["nearest_tss_gene"]) == {"XXX"} @@ -306,6 +317,7 @@ def test_a_column_annotate_sites_never_made_is_written_empty(tmp_path): "end": [30], "strand": ["+"], "location": ["intergenic"], + "feature": [""], "region": [""], "sense": [""], "ensg": [None], @@ -424,7 +436,8 @@ def test_build_db_then_run_from_the_command_line(tmp_path, capsys): assert status == 0 assert capsys.readouterr().out.splitlines() == [ - "4 rows for 3 sites in 2 genes: 2 exon, 1 intron, 1 intergenic, 0 in a promoter", + "4 rows for 3 sites in 2 genes: 3 genic (2 exon, 1 intron)," + " 1 intergenic, 0 in a promoter", "Outputs in out", ] assert matches(written(tmp_path / "out"), EXPECTED_FILES) diff --git a/tests/test_guide_alignment_pipeline.py b/tests/test_guide_alignment_pipeline.py index 5ff97a6..1b4491b 100644 --- a/tests/test_guide_alignment_pipeline.py +++ b/tests/test_guide_alignment_pipeline.py @@ -48,7 +48,7 @@ "start", "end", "strand", - "cut_site", + "target_site", "pam_matches", "mapped_sequence", ) @@ -156,11 +156,11 @@ def test_run_writes_every_output(setup): assert g_a["end"] == str(a_start + 20) assert g_a["strand"] == "+" assert g_a["pam_matches"] == "TRUE" - assert g_a["cut_site"] == str(a_start + 17) + assert g_a["target_site"] == str(a_start + 17) assert g_a["guide_id"] == SPACER_A + "NGG" # The genomic PAM is the tail of mapped_sequence; the pattern ends guide_id. assert g_a["mapped_sequence"] == SPACER_A + "TGG" - assert (g_a["n_matches"], g_a["n_cut_sites"]) == ("1", "1") + assert (g_a["n_matches"], g_a["n_target_sites"]) == ("1", "1") assert result.summaries[0].site_pam == "TGG" assert result.summaries[0].has_perfect_site is True assert result.summaries[0].extra_5p is False @@ -168,10 +168,10 @@ def test_run_writes_every_output(setup): (g_c,) = guides["g_c"] assert g_c["site_class"] == "unique" assert (g_c["start"], g_c["strand"]) == (str(c_start), "-") - assert g_c["cut_site"] == str(c_start + 3) + assert g_c["target_site"] == str(c_start + 3) # Two sites, so two rows, with the guide's own columns repeated on each. assert {row["site_class"] for row in guides["g_g"]} == {"multi"} - assert {row["n_cut_sites"] for row in guides["g_g"]} == {"2"} + assert {row["n_target_sites"] for row in guides["g_g"]} == {"2"} assert len(guides["g_g"]) == 2 (g_nt,) = guides["g_nt"] assert g_nt["site_class"] == "no_site" @@ -246,7 +246,7 @@ def test_the_scan_is_exact_so_a_mismatched_site_is_never_reported(setup): rows = table(setup["outdir"] / "alignments" / "guide_genome_annotation.tsv") assert str(nag_start) not in {row["start"] for row in rows} guides = by_alias(rows) - assert guides["g_c"][0]["n_cut_sites"] == "1" + assert guides["g_c"][0]["n_target_sites"] == "1" assert guides["g_c"][0]["n_matches"] == "1" assert [s.min_mismatches for s in result.summaries if s.alias == "g_c"] == [0] assert { @@ -331,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["n_matches"], row["n_cut_sites"]) == ("0", "0") + 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" @@ -346,7 +346,7 @@ def test_a_match_without_a_pam_is_kept_out_of_the_bed_but_written_to_the_tsv(tmp assert (site.pam, site.pam_class) == ("TTT", "none") (row,) = table(result.outdir / "alignments" / "guide_genome_annotation.tsv") assert (row["chr"], row["start"], row["pam_matches"]) == ("c1", "4", "FALSE") - assert (row["n_matches"], row["n_cut_sites"]) == ("1", "0") + assert (row["n_matches"], row["n_target_sites"]) == ("1", "0") assert row["mapped_sequence"] == SPACER_A + "TTT" # Nothing can be cut there, so the BED is empty. assert (result.outdir / "alignments" / "guide_mapping.bed").read_text() == "" diff --git a/tests/test_guide_alignment_summary.py b/tests/test_guide_alignment_summary.py index 10b28e0..9992f8f 100644 --- a/tests/test_guide_alignment_summary.py +++ b/tests/test_guide_alignment_summary.py @@ -61,14 +61,14 @@ "spacer", "site_class", "n_matches", - "n_cut_sites", + "n_target_sites", "orientation", "orientation_evidence", "chr", "start", "end", "strand", - "cut_site", + "target_site", "pam_matches", "mapped_sequence", ] @@ -196,7 +196,7 @@ def test_classes_counts_and_coordinates(guides, sites): "dup": "unique", } uniq = by_alias["uniq"] - # The PAM-less match counts as a match, but neither as a cut site nor + # The PAM-less match counts as a match, but neither as a target site nor # towards multi. assert (uniq.n_matches, uniq.n_sites) == (2, (1, 0, 0, 0)) assert (uniq.guide_chr, uniq.guide_start, uniq.guide_end, uniq.strand) == ( @@ -205,15 +205,19 @@ def test_classes_counts_and_coordinates(guides, sites): 120, "+", ) - assert (uniq.site_pam, uniq.cut_site, uniq.has_perfect_site) == ("TGG", 117, True) + assert (uniq.site_pam, uniq.target_site, uniq.has_perfect_site) == ( + "TGG", + 117, + True, + ) assert (uniq.mapped_sequence, uniq.input_sequence) == (SPACER_A + "TGG", SPACER_A) - assert by_alias["uniq_minus"].cut_site == 203 + assert by_alias["uniq_minus"].target_site == 203 assert by_alias["uniq_minus"].site_pam == "AGG" multi = by_alias["multi"] assert ( multi.has_perfect_site, multi.guide_chr, - multi.cut_site, + multi.target_site, multi.site_pam, multi.mapped_sequence, ) == (True, None, None, None, None) @@ -281,7 +285,7 @@ def test_a_duplicate_gets_the_first_guide_s_sites_and_verdict(guides, sites): assert guides[6].duplicate_of == "uniq" assert guides[0].duplicate_of is None assert dup.orientation == uniq.orientation == "forward" - for field in ("site_class", "guide_start", "cut_site", "n_sites", "site_pam"): + for field in ("site_class", "guide_start", "target_site", "n_sites", "site_pam"): assert getattr(dup, field) == getattr(uniq, field) @@ -363,7 +367,7 @@ def test_guides_tsv(tmp_path, guides, sites): ("700", "-"), ] assert {row["site_class"] for row in multi} == {"multi"} - assert {row["n_cut_sites"] for row in multi} == {"2"} + assert {row["n_target_sites"] for row in multi} == {"2"} def test_a_protospacer_matching_twice_with_one_pam_is_unique(tmp_path, guides, sites): @@ -375,7 +379,7 @@ def test_a_protospacer_matching_twice_with_one_pam_is_unique(tmp_path, guides, s assert (uniq.n_matches, uniq.n_sites[0], uniq.site_class) == (2, 1, "unique") rows = rows_by_alias(table)["uniq"] assert [ - (row["n_matches"], row["n_cut_sites"], row["site_class"]) for row in rows + (row["n_matches"], row["n_target_sites"], row["site_class"]) for row in rows ] == [("2", "1", "unique")] * 2 assert [row["pam_matches"] for row in rows] == ["TRUE", "FALSE"] @@ -388,7 +392,7 @@ def test_two_pam_bearing_matches_are_multi_with_a_bed_row_each(tmp_path, guides, multi = summaries[2] assert (multi.site_class, multi.n_sites[0], multi.n_matches) == ("multi", 2, 2) rows = rows_by_alias(table)["multi"] - assert [row["n_cut_sites"] for row in rows] == ["2", "2"] + assert [row["n_target_sites"] for row in rows] == ["2", "2"] named = [row for row in bed if row[3] == SPACER_G] assert named == [ ["chr1", "600", "620", SPACER_G, "0", "+"], @@ -405,7 +409,7 @@ def test_a_guide_matching_nowhere_gets_one_row_with_empty_site_columns( (row,) = rows_by_alias(table)["nothing"] assert [row[name] for name in SITE_COLUMN_NAMES] == [""] * 7 - assert (row["n_matches"], row["n_cut_sites"]) == ("0", "0") + 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 assert not [line for line in bed if line[3] == SPACER_NOTHING] @@ -476,7 +480,7 @@ def test_guides_tsv_rejects_summaries_that_are_not_one_per_guide( write_guides_tsv(tmp_path / "guides.tsv", guides, summaries[:-1], sites) -def test_the_cut_site_column_is_three_bases_from_the_pam(tmp_path, guides, sites): +def test_the_target_site_column_is_three_bases_from_the_pam(tmp_path, guides, sites): # The boundary is end - 3 on + and start + 3 on -, counted from the PAM # either way. summaries = summarized(guides, sites) @@ -485,12 +489,15 @@ def test_the_cut_site_column_is_three_bases_from_the_pam(tmp_path, guides, sites rows = rows_by_alias(table) assert by_alias["uniq"].strand == "+" - assert by_alias["uniq"].guide_end - by_alias["uniq"].cut_site == 3 + assert by_alias["uniq"].guide_end - by_alias["uniq"].target_site == 3 assert by_alias["uniq_minus"].strand == "-" - assert by_alias["uniq_minus"].cut_site - by_alias["uniq_minus"].guide_start == 3 - assert (by_alias["uniq"].cut_site, by_alias["uniq_minus"].cut_site) == (117, 203) - assert rows["uniq"][0]["cut_site"] == "117" - assert rows["uniq_minus"][0]["cut_site"] == "203" + assert by_alias["uniq_minus"].target_site - by_alias["uniq_minus"].guide_start == 3 + assert (by_alias["uniq"].target_site, by_alias["uniq_minus"].target_site) == ( + 117, + 203, + ) + assert rows["uniq"][0]["target_site"] == "117" + assert rows["uniq_minus"][0]["target_site"] == "203" def test_qc_tsv_counts_sum_to_the_entries(tmp_path, guides, sites): @@ -545,7 +552,7 @@ def test_qc_tsv_rejects_a_clashing_section(tmp_path, guides, sites): ) -def test_coordinates_and_cut_sites_from_a_hand_built_reference(tmp_path): +def test_coordinates_and_target_sites_from_a_hand_built_reference(tmp_path): guides, contig, fasta = planted(tmp_path) assert [guide.extra_5p for guide in guides] == [True] * 4 @@ -556,7 +563,7 @@ def test_coordinates_and_cut_sites_from_a_hand_built_reference(tmp_path): # - : PAM [150, 153), protospacer [153, 173), b = 156. # + with the added G on a C: protospacer [201, 221), b = 218. # + with the added G on a G: protospacer [235, 255), b = 252. - assert [(s.guide_start, s.guide_end, s.cut_site) for s in summaries] == [ + assert [(s.guide_start, s.guide_end, s.target_site) for s in summaries] == [ (100, 120, 117), (153, 173, 156), (201, 221, 218), @@ -578,7 +585,10 @@ def test_coordinates_and_cut_sites_from_a_hand_built_reference(tmp_path): table, bed = write_tables(tmp_path, guides, summaries, scan.sites) rows = rows_by_alias(table) - assert [rows[alias][0]["cut_site"] for alias in ("plus", "minus")] == ["117", "156"] + assert [rows[alias][0]["target_site"] for alias in ("plus", "minus")] == [ + "117", + "156", + ] assert sorted(line[1] for line in bed) == ["100", "153", "201", "235"]