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"]