From 45c8d1eeefb838410e9c849de9afe6c7571c8211 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Thu, 15 Oct 2026 12:02:47 -0400 Subject: [PATCH 1/2] ADD the measured gene, the element's location, and what the channel is Three gaps, each found by reading the assembled object and asking where something had gone. `guide_element_mapping.tsv` named no gene, so the mapping could not stand on its own. It now takes the same --genes pairs table igvf-export does, through the same lookup, so the two cannot disagree: 3,380 of 3,890 promoter rows and 180 of 379 e2g_validated rows gain one. `genomic_element` cannot say whether an element sits in a gene, because promoter and enhancer win on precedence before the gene body is tested. The cross-tab shows why that matters: exon, intron and distal are determined, but 499 promoter and 110 enhancer guides sit outside any gene. `element_location` is computed from the overlap igvf-export already had, so no new input and no second source of truth. And the object said nothing about which channel it was. --channel NAME=VALUE joins every cell's obs and is kept in uns, so replicate, day, chip and batch travel with the counts. Stated rather than parsed: the names live in FASTQ paths, and reading them there would tie the tool to one lab's scheme. Co-Authored-By: Claude Opus 5 --- .../cli/element_mapping.py | 10 +++++++ .../cli/mudata_builder.py | 23 ++++++++++++++++ src/crispr_analysis_utils/element_mapping.py | 27 +++++++++++++++++++ src/crispr_analysis_utils/igvf_export.py | 14 ++++++++++ src/crispr_analysis_utils/mudata_builder.py | 16 ++++++++++- tests/test_element_mapping.py | 8 +++--- 6 files changed, 93 insertions(+), 5 deletions(-) diff --git a/src/crispr_analysis_utils/cli/element_mapping.py b/src/crispr_analysis_utils/cli/element_mapping.py index fe6162f..c528a26 100644 --- a/src/crispr_analysis_utils/cli/element_mapping.py +++ b/src/crispr_analysis_utils/cli/element_mapping.py @@ -43,6 +43,15 @@ def add_arguments(parser: argparse.ArgumentParser) -> None: "which fills the distance_to_tss column." ), ) + run_parser.add_argument( + "--genes", + default=None, + help=( + "A pairs table of element to gene, which fills " + "putative_target_genes. The same file `cau igvf-export run` " + "takes, so the two tables name the same genes." + ), + ) def run(args: argparse.Namespace) -> int: @@ -56,6 +65,7 @@ def run(args: argparse.Namespace) -> int: args.elements, args.outdir, annotation=args.annotation, + genes=args.genes, ) except (FileNotFoundError, ElementError, ValueError, OSError, KeyError) as error: raise CommandError(str(error)) from error diff --git a/src/crispr_analysis_utils/cli/mudata_builder.py b/src/crispr_analysis_utils/cli/mudata_builder.py index 9304d3e..fd931a8 100644 --- a/src/crispr_analysis_utils/cli/mudata_builder.py +++ b/src/crispr_analysis_utils/cli/mudata_builder.py @@ -67,6 +67,17 @@ def add_arguments(parser: argparse.ArgumentParser) -> None: "--element-mapping", help="guide_element_mapping.tsv. Carried whole in guide.uns.", ) + run_parser.add_argument( + "--channel", + action="append", + default=[], + metavar="NAME=VALUE", + help=( + "What this channel is, repeatable: --channel replicate=1 " + "--channel day=4 --channel chip=A. One value for the whole " + "channel, joined to every cell's obs and kept in uns." + ), + ) run_parser.add_argument( "--covariates", help=( @@ -76,6 +87,17 @@ def add_arguments(parser: argparse.ArgumentParser) -> None: ) +def _channel(pairs: list[str]) -> dict[str, str]: + """NAME=VALUE arguments as a mapping, refusing one without the equals.""" + out = {} + for pair in pairs: + name, sep, value = pair.partition("=") + if not sep or not name: + raise CommandError(f"--channel takes NAME=VALUE, got {pair!r}.") + out[name] = value + return out + + def run(args: argparse.Namespace) -> int: """Run the action.""" from ..mudata_builder import BuildError, run @@ -92,6 +114,7 @@ def run(args: argparse.Namespace) -> int: gtf_accession=args.gtf_accession, genome_accession=args.genome_accession, covariates=args.covariates, + metadata=_channel(args.channel), guide_annotation=args.guide_annotation, element_mapping=args.element_mapping, ) diff --git a/src/crispr_analysis_utils/element_mapping.py b/src/crispr_analysis_utils/element_mapping.py index 79a2d39..db808db 100644 --- a/src/crispr_analysis_utils/element_mapping.py +++ b/src/crispr_analysis_utils/element_mapping.py @@ -18,6 +18,7 @@ import logging import os from collections import Counter +from collections.abc import Mapping from dataclasses import dataclass, field from datetime import UTC, datetime from pathlib import Path @@ -65,6 +66,12 @@ " than trusting the class alone. From guide_gene_annotation.tsv, so" " empty without it." ), + "putative_target_genes": ( + "Semicolon-separated ENSEMBL ids of the genes measured against this" + " element, from the pairs table given to --genes. The same genes the" + " IGVF export carries, repeated here so the mapping stands alone." + " Empty for an element with no measured pair, and without --genes." + ), } """Every column of the mapping table, and what it holds.""" @@ -137,6 +144,7 @@ def map_to_elements( elements: pd.DataFrame, guides: pd.Series | None = None, distances: pd.DataFrame | None = None, + genes: Mapping[str, list[str]] | None = None, ) -> pd.DataFrame: """Map every site to the elements it overlaps, and classify it. @@ -224,6 +232,15 @@ def map_to_elements( for column in ("start", "end"): frame[column] = frame[column].astype("Float64").astype("Int64") frame = _add_distance(frame, distances) + # The same lookup the IGVF export uses, so the two tables cannot disagree + # about which gene an element was measured against. + from .igvf_export import _measured_genes + + frame["putative_target_genes"] = ( + [";".join(_measured_genes(genes, eid)) or pd.NA for eid in frame["element_id"]] + if genes + else pd.NA + ) return frame.reset_index(drop=True).reindex(columns=list(COLUMNS)) @@ -515,6 +532,7 @@ def run( outdir: str | os.PathLike[str], *, annotation: str | os.PathLike[str] | pd.DataFrame | None = None, + genes: str | os.PathLike[str] | None = None, ) -> RunResult: """Map a library to the elements a screen targets, and write the table. @@ -599,11 +617,20 @@ def run( else: logger.warning("No annotation given, so distance_to_tss will be empty.") + measured = None + if genes is not None: + from .igvf_export import _read_genes + + measured = _read_genes(genes) + state["inputs"]["genes"] = str(genes) + state["inputs"]["genes_md5"] = file_md5(genes) + mapping = map_to_elements( sites, regions, guides=library, distances=distances, + genes=measured, ) write_mapping_tsv(outdir / RESULTS_DIR / TABLE, mapping) write_columns_yaml(outdir / RESULTS_DIR / "guide_element_mapping_columns.yaml") diff --git a/src/crispr_analysis_utils/igvf_export.py b/src/crispr_analysis_utils/igvf_export.py index 932fcb9..9896924 100644 --- a/src/crispr_analysis_utils/igvf_export.py +++ b/src/crispr_analysis_utils/igvf_export.py @@ -56,6 +56,11 @@ "The protospacer adjacent motif as the genome holds it, 5' to 3' on" " the targeted strand, immediately 3' of the spacer." ), + "element_location": ( + "genic when the element overlaps a protein-coding gene, intergenic" + " when it does not. Independent of genomic_element, which tests" + " promoter and enhancer before it ever looks at the gene body." + ), "genomic_element": ( "What the targeted element is: promoter, enhancer, distal element," " exon or intron. Decided from the genome." @@ -82,6 +87,8 @@ DISTAL = "distal element" EXON = "exon" INTRON = "intron" +GENIC = "genic" +INTERGENIC = "intergenic" ELEMENTS = (PROMOTER, ENHANCER, DISTAL, EXON, INTRON) """The `genomic_element` values this writes. The format allows more.""" @@ -192,6 +199,12 @@ def classify_elements( for row in frame.itertuples() ] frame["nearest_ensg"] = frame["element_id"].map(nearest_gene) + # Whether the element itself sits in a gene. genomic_element cannot say: + # promoter and enhancer win on precedence before the gene body is tested, + # so this is computed from the same overlap and kept beside it. + frame["element_location"] = np.where( + frame["element_id"].isin(in_gene), GENIC, INTERGENIC + ) exon_ids = _exon_ids(ranges, store, pr) introns = _intron_spans(ranges, store, pr) @@ -466,6 +479,7 @@ def build_table( ) for column, source in ( + ("element_location", "element_location"), ("genomic_element", "genomic_element"), ("intended_target_name", "intended_target_name"), ("intended_target_chr", "intended_target_chr"), diff --git a/src/crispr_analysis_utils/mudata_builder.py b/src/crispr_analysis_utils/mudata_builder.py index 0396cdb..afa7d44 100644 --- a/src/crispr_analysis_utils/mudata_builder.py +++ b/src/crispr_analysis_utils/mudata_builder.py @@ -279,6 +279,7 @@ def build( reference: dict | None = None, covariates: pd.DataFrame | None = None, tables: dict[str, pd.DataFrame] | None = None, + metadata: dict[str, str] | None = None, ): """Assemble the two matrices into one MuData. @@ -297,6 +298,11 @@ def build( covariates One row per cell, indexed by barcode. Its columns join the root ``obs``. Nothing is computed; whatever is given is carried. + metadata + What this channel is, as name to value: replicate, day, chip, batch, + whatever the screen records. One value for the whole channel, so it + joins every cell's ``obs`` and is also kept in ``uns`` for anything + reading the object rather than the table. tables Whole tables to carry in ``mod['guide'].uns``, by name. The gene annotation is one row per guide per site per gene, so it cannot go in @@ -349,6 +355,9 @@ def build( data = md.MuData({"gene": rna, "guide": guide}) data.uns["reference"] = dict(reference or {}) data.uns["measurement_set"] = measurement_set + data.uns["channel"] = dict(metadata or {}) + for name, value in (metadata or {}).items(): + data.obs[name] = value if covariates is not None and not covariates.empty: joined = covariates.reindex( @@ -393,6 +402,7 @@ def run( covariates: str | os.PathLike[str] | None = None, guide_annotation: str | os.PathLike[str] | None = None, element_mapping: str | os.PathLike[str] | None = None, + metadata: dict[str, str] | None = None, ) -> RunResult: """Assemble one channel's matrices into a MuData, and write it. @@ -443,7 +453,10 @@ def run( "finished": None, "cau_version": __version__, "inputs": {}, - "parameters": {"measurement_set": measurement_set}, + "parameters": { + "measurement_set": measurement_set, + "channel": dict(metadata or {}), + }, "counts": {}, } handler = logging.FileHandler(outdir / "logs" / LOG_FILE, encoding="utf-8") @@ -510,6 +523,7 @@ def run( reference=reference, covariates=covariate_table, tables=tables, + metadata=metadata, ) data.write(outdir / RESULTS_DIR / TABLE) diff --git a/tests/test_element_mapping.py b/tests/test_element_mapping.py index f4264c0..7c9709e 100644 --- a/tests/test_element_mapping.py +++ b/tests/test_element_mapping.py @@ -447,8 +447,8 @@ def test_whole_numbers_are_written_without_a_trailing_zero(tmp_path): ) write_mapping_tsv(tmp_path / "m.tsv", frame) lines = (tmp_path / "m.tsv").read_text().splitlines() - assert lines[1] == "a\tchr1\t1000\t1020\t+\tE1\tenhancer\ttargeting\t47" - assert lines[2] == "b\t\t\t\t\t\t\tnon_targeting\t" + assert lines[1] == "a\tchr1\t1000\t1020\t+\tE1\tenhancer\ttargeting\t47\t" + assert lines[2] == "b\t\t\t\t\t\t\tnon_targeting\t\t" def test_a_written_mapping_holds_no_decimal_point(with_distances, tmp_path): @@ -464,7 +464,7 @@ def test_a_written_mapping_holds_no_decimal_point(with_distances, tmp_path): def test_a_non_targeting_row_is_written_with_empty_cells(mapping, tmp_path): write_mapping_tsv(tmp_path / "m.tsv", mapping) lines = (tmp_path / "m.tsv").read_text().splitlines() - assert lines[-1] == "never\t\t\t\t\t\t\tnon_targeting\t" + assert lines[-1] == "never\t\t\t\t\t\t\tnon_targeting\t\t" def test_a_column_map_to_elements_never_made_is_written_empty(tmp_path): @@ -472,7 +472,7 @@ def test_a_column_map_to_elements_never_made_is_written_empty(tmp_path): write_mapping_tsv(tmp_path / "bare.tsv", bare) lines = (tmp_path / "bare.tsv").read_text().splitlines() assert lines[0].split("\t") == list(COLUMNS) - assert lines[1] == "a\t\t\t\t\t\t\tnon_targeting\t" + assert lines[1] == "a\t\t\t\t\t\t\tnon_targeting\t\t" def test_columns_yaml_parses_and_documents_the_table(mapping, tmp_path): From ecced4abcedd094b08143b5fb04476044ec0db7d Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Thu, 15 Oct 2026 12:04:57 -0400 Subject: [PATCH 2/2] UPDATE the skills and docs for the three new columns and flags Each addition is something a user has to be told, not inferred: --genes on element-mapping, --channel on mudata-builder, and why element_location exists when genomic_element already names a class. The last carries the cross-tab, because the argument for the column is the 499 promoter and 110 enhancer guides that sit outside any gene while every exon and intron guide is inside one. Co-Authored-By: Claude Opus 5 --- docs/element-mapping.md | 2 ++ docs/igvf-export.md | 9 +++++++++ docs/mudata-builder.md | 5 +++++ .../skills/cau-element-mapping/SKILL.md | 4 ++++ .../cau-element-mapping/references/outputs.md | 1 + .../skills/cau-igvf-export/SKILL.md | 9 +++++++++ .../cau-igvf-export/references/outputs.md | 1 + .../skills/cau-mudata-builder/SKILL.md | 20 ++++++++++++++++++- 8 files changed, 50 insertions(+), 1 deletion(-) diff --git a/docs/element-mapping.md b/docs/element-mapping.md index f78d25b..2eb0c82 100644 --- a/docs/element-mapping.md +++ b/docs/element-mapping.md @@ -18,6 +18,8 @@ cau element-mapping run \ --outdir results/ ``` +`--genes` is optional and fills `putative_target_genes` from the same pairs +table `cau igvf-export` takes, so the two tables name the same genes. `--annotation` is optional and fills `distance_to_tss`; the classes do not depend on it. Nothing is written outside `--outdir`. diff --git a/docs/igvf-export.md b/docs/igvf-export.md index 1d6e7a9..7860fd8 100644 --- a/docs/igvf-export.md +++ b/docs/igvf-export.md @@ -89,3 +89,12 @@ passed to it as it is. `intended_target_name` is what joins the guide table to the pipeline's pairs-to-test file, so the two have to name elements the same way. + +## genomic_element does not say genic + +`genomic_element` is decided by precedence -- promoter, then exon, then +intron, then distal element -- so promoter and enhancer win before the gene +body is tested and the class cannot say whether the element sits in a gene. +`element_location` does, computed from the same overlap. On the hiPSC-to-EC +library 499 promoter and 110 enhancer guides are intergenic, while every exon +and intron guide is genic. diff --git a/docs/mudata-builder.md b/docs/mudata-builder.md index bc9c735..8c439bc 100644 --- a/docs/mudata-builder.md +++ b/docs/mudata-builder.md @@ -21,6 +21,11 @@ cau mudata-builder run \ --outdir results/ ``` +`--channel NAME=VALUE` is repeatable and says what the channel is: replicate, +day, chip, batch. It joins every cell's `obs` and is kept in `uns['channel']`. +`--guide-annotation` and `--element-mapping` carry those tables whole into +`mod['guide'].uns`, since neither is one row per guide. + Nothing is written outside `--outdir`. ## What it assembles 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 e723078..ac306f5 100644 --- a/src/crispr_analysis_utils/skills/cau-element-mapping/SKILL.md +++ b/src/crispr_analysis_utils/skills/cau-element-mapping/SKILL.md @@ -70,6 +70,10 @@ in full. 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. +- **`--genes` fills `putative_target_genes`**, from the same pairs table + `cau igvf-export run` takes and through the same lookup, so the two tables + name the same genes. Without it the column is empty; with it, only elements + the pairs table covers get one. - **`element_id` comes from the BED's fourth column** when it has one, so an id used elsewhere keeps its spelling. Otherwise it is `chr:start-end`, 1-based, as a locus is read aloud. 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 6c0d97b..0414da0 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 @@ -24,6 +24,7 @@ site and no element. | `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 target site to the nearest protein-coding TSS, carried from the annotation. | +| `putative_target_genes` | Semicolon-separated ENSEMBL ids measured against the element, from `--genes`. The same genes the IGVF export carries. | ### The classes are per site diff --git a/src/crispr_analysis_utils/skills/cau-igvf-export/SKILL.md b/src/crispr_analysis_utils/skills/cau-igvf-export/SKILL.md index 6072310..fe40130 100644 --- a/src/crispr_analysis_utils/skills/cau-igvf-export/SKILL.md +++ b/src/crispr_analysis_utils/skills/cau-igvf-export/SKILL.md @@ -75,6 +75,15 @@ a promoter and the genome agrees. | `intron` | the intron's coordinates | | `enhancer`, `distal element` | the element's coordinates | +## genomic_element does not say genic + +`genomic_element` is decided by precedence: promoter, then exon, then intron, +then distal element. Promoter and enhancer win before the gene body is ever +tested, so the class cannot say whether the element sits in a gene. +`element_location` does, from the same overlap. On one library 499 promoter +and 110 enhancer guides are intergenic, while every exon and intron guide is +genic. + ## Reading the result `--outdir` gets `igvf/`, `run_metadata/` and `logs/`. The results are one diff --git a/src/crispr_analysis_utils/skills/cau-igvf-export/references/outputs.md b/src/crispr_analysis_utils/skills/cau-igvf-export/references/outputs.md index dc5128a..1502158 100644 --- a/src/crispr_analysis_utils/skills/cau-igvf-export/references/outputs.md +++ b/src/crispr_analysis_utils/skills/cau-igvf-export/references/outputs.md @@ -25,6 +25,7 @@ guide and one spacer. | `guide_chr`, `guide_start`, `guide_end` | The protospacer, 0-based and half-open, **without** the PAM. | | `strand` | The strand targeted. | | `pam` | The PAM as the genome holds it, 5' to 3' on that strand. | +| `element_location` | `genic` or `intergenic`: whether the element overlaps a protein-coding gene. Independent of `genomic_element`, which tests promoter and enhancer first. | | `genomic_element` | `promoter`, `enhancer`, `distal element`, `exon` or `intron`. | | `intended_target_name` | Names the target; what it holds depends on the element. | | `intended_target_chr/start/end` | The element's coordinates. | diff --git a/src/crispr_analysis_utils/skills/cau-mudata-builder/SKILL.md b/src/crispr_analysis_utils/skills/cau-mudata-builder/SKILL.md index e4739c5..47b72ab 100644 --- a/src/crispr_analysis_utils/skills/cau-mudata-builder/SKILL.md +++ b/src/crispr_analysis_utils/skills/cau-mudata-builder/SKILL.md @@ -20,7 +20,10 @@ cau mudata-builder run \ --guide-metadata results/igvf/guide_metadata.tsv \ --store gencode.v43.annotation.parquet \ --measurement-set IGVFDS8842PXCC \ + --guide-annotation results/annotations/guide_gene_annotation.tsv \ + --element-mapping results/elements/guide_element_mapping.tsv \ --gtf-accession IGVFFI9561BASO \ + --channel replicate=1 --channel day=4 --channel chip=A \ --outdir results/ ``` @@ -34,6 +37,10 @@ cau mudata-builder run \ - **Which channel?** `--measurement-set` goes into the cell index as `{barcode}_{measurement_set}`, so two channels concatenate without their barcodes colliding. There is no sensible default. +- **What is this channel?** `--channel NAME=VALUE`, repeatable: replicate, + day, chip, batch, whatever the screen records. Ask rather than read it out + of a FASTQ path, which ties the answer to one lab's naming. One value for + the whole channel, so it joins every cell. - **The portal ids.** `--gtf-accession` and `--genome-accession` take IGVF accessions; the CRISPR pipeline's inputs carry them. Without one the object records the file's name instead, so the field is never empty and never a @@ -49,7 +56,9 @@ cau mudata-builder run \ | `mod['gene'].var` | `symbol`, `gene_type`, `gene_chr`, `gene_start`, `gene_end` | | `mod['guide'].X` | the guide counts, untouched | | `mod['guide'].var` | every column of `guide_metadata.tsv` | -| `obs` | `barcode` and `measurement_set`, plus any `--covariates` given | +| `obs` | `barcode`, `measurement_set`, every `--channel` value, and any `--covariates` given | +| `mod['guide'].uns` | `annotation` and `element_mapping`, whole, from `--guide-annotation` and `--element-mapping` | +| `uns['channel']` | the same `--channel` values, for anything reading the object rather than the table | | `obsm['gene']`, `obsm['guide']` | which modality holds each cell | | `uns['reference']` | the GTF and genome: name, md5 and accession | @@ -60,6 +69,15 @@ statistics, no QC, no filtering, no doublet calls, no guide assignment, no pairs to test. Those are later steps, and each is a judgement this one has no business making. +## The tables that cannot go in var + +`guide_gene_annotation.tsv` is one row per guide per site per gene, 45,501 +rows for 43,766 guides, and `guide_element_mapping.tsv` one per guide per +element. Neither is one row per guide, so neither fits `var`. Both are +carried whole in `mod['guide'].uns`: the signed distances, the promoter +transcripts and the other-gene-type overlaps are there to be joined, and +nothing about how to collapse them has been decided. + ## Four things worth knowing - **Cells are the union, not the intersection.** A barcode seen in one