Carry the measured gene, the element's location, and the channel

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`.
 

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.

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

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

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,
         )

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

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"),

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)
 

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.

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
 

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

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. |

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

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,15 +464,15 @@ 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):
     bare = pd.DataFrame({"guide_id": ["a"], "target_class": [NON_TARGETING]})
     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):