From 5f58e9a341d8848b56b96d818c871a57597dfbfe Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Thu, 15 Oct 2026 11:40:32 -0400 Subject: [PATCH 1/3] ADD oligo and n_oligos, so a merged guide says its counts are a sum Two library entries of one spacer are one feature in any reference built on the spacer, and their reads land together. This library has 12 such pairs, each a 20-mer and the same spacer with a 5' G added for U6 transcription. They cannot be separated, and not for want of a better reference: the leader ends in a poly-G run, so the extra base is indistinguishable from the leader's own, and splitcode extracts a fixed 20 nt anchored on the scaffold, which both entries yield identically. Cell Ranger keeps all 43,778 features and assigns the 21-mer zero reads in every one of the 12 pairs; its counts for those guides are the 20-mer's alone. So the spacer stays what the reads can resolve, and the table says which guides carry a sum: `oligo` records what was ordered, `n_oligos` how many entries share the spacer. Above one, a guide is not comparable with one ordered once. Both columns travel into guide_metadata.tsv. Measured on the hiPSC-to-EC library: 24 entries flagged over 12 spacers, and 33,164 of 43,778 entries carry an extra 5' base, so oligo differs from spacer for three quarters of the library. Co-Authored-By: Claude Opus 5 --- .../guide_alignment/summary.py | 32 ++++++++++++++- src/crispr_analysis_utils/igvf_export.py | 15 +++++++ .../skills/cau-guide-alignment/SKILL.md | 7 ++++ .../cau-guide-alignment/references/outputs.md | 2 + .../cau-igvf-export/references/outputs.md | 2 + tests/test_guide_alignment_summary.py | 39 +++++++++++++++++-- tests/test_igvf_export.py | 12 ++++-- 7 files changed, 102 insertions(+), 7 deletions(-) diff --git a/src/crispr_analysis_utils/guide_alignment/summary.py b/src/crispr_analysis_utils/guide_alignment/summary.py index a9c919e..ecca9ab 100644 --- a/src/crispr_analysis_utils/guide_alignment/summary.py +++ b/src/crispr_analysis_utils/guide_alignment/summary.py @@ -58,6 +58,13 @@ class GuideSummary: The name the guide table gives this entry. spacer The spacer, 5' to 3'. + oligo + The sequence the library ordered, without the PAM: the spacer plus + any base 5' of it. + n_oligos + How many library entries carry this spacer. More than one means a + feature reference built on the spacer merges them, so the counts are + a sum and are not comparable with a guide ordered once. mapped_sequence The guide's bases that pair with the genome at the unique site, as in `Site.mapped_sequence`; None unless the class is ``unique``. @@ -104,6 +111,8 @@ class GuideSummary: guide_id: str alias: str spacer: str + oligo: str + n_oligos: int mapped_sequence: str | None site_pam: str | None has_perfect_site: bool @@ -164,6 +173,10 @@ class that is neither one of the patterns nor ``none``, a guide is a patterns = validate_pams([pam]) _check_max_mismatches(max_mismatches) primary = patterns[0] + # Two entries of one spacer are one feature in any reference built on the + # spacer, so their reads land together. Count them here so the table can + # say which guides carry a sum. + per_spacer = Counter(guide.spacer for guide in guides) summaries = [] for guide, source in zip(guides, _source_rows(guides), strict=True): n_matches = 0 @@ -205,6 +218,8 @@ class that is neither one of the patterns nor ``none``, a guide is a summaries.append( GuideSummary( guide_id=guide.guide_id, + oligo=guide.aligned_sequence, + n_oligos=per_spacer[guide.spacer], alias=guide.alias, spacer=guide.spacer, mapped_sequence=unique.mapped_sequence if unique else None, @@ -267,7 +282,7 @@ def fan_out_sites( GROUPS = ( - ("from the guide table, before anything is aligned", 4), + ("from the guide table, before anything is aligned", 6), ("about the whole guide, repeated on each of its rows", 5), ("about this row's match; empty when the guide matches nowhere", 7), ) @@ -287,6 +302,19 @@ def fan_out_sites( "spacer": ( "The protospacer: guide_id without the PAM pattern or any base 5' of it." ), + "oligo": ( + "The sequence the library ordered, without the PAM. It differs from" + " the spacer when the entry carries a base 5' of the protospacer," + " such as a G added for U6 transcription." + ), + "n_oligos": ( + "How many library entries share this spacer. Above one, a feature" + " reference built on the spacer cannot tell them apart, so the" + " counts are a sum: expect roughly that multiple of a guide ordered" + " once. Here the second entry is the same spacer with a 5' G, which" + " no read can be assigned to on its own, because the leader ends in" + " a poly-G run." + ), # about the whole guide "site_class": ( "unique, multi or no_site: how many PAM-bearing places this guide matches." @@ -393,6 +421,8 @@ def write_guides_tsv( summary.alias, guide_accession(summary.spacer), summary.spacer, + summary.oligo, + summary.n_oligos, summary.site_class, summary.n_matches, summary.n_sites[0], diff --git a/src/crispr_analysis_utils/igvf_export.py b/src/crispr_analysis_utils/igvf_export.py index 32b27c8..c141b0a 100644 --- a/src/crispr_analysis_utils/igvf_export.py +++ b/src/crispr_analysis_utils/igvf_export.py @@ -30,6 +30,15 @@ COLUMNS = { "guide_id": "Unique name of the guide.", "spacer": "The targeting sequence of the gRNA, read 5' to 3', in DNA bases.", + "oligo": ( + "The sequence the library ordered, without the PAM: the spacer plus" + " any base 5' of it, such as a G added for U6 transcription." + ), + "n_oligos": ( + "How many library entries share this spacer. Above one, they are one" + " feature in the counts and the reads cannot be told apart, so the" + " counts are a sum." + ), "targeting": ( "TRUE when the guide targets a locus relevant to the experiment." " FALSE for a safe-targeting or non-targeting control." @@ -391,6 +400,12 @@ def build_table( out = pd.DataFrame(index=guides.index) out["guide_id"] = guides["alias"] out["spacer"] = guides["spacer"] + # Carried so the submitted table says which guides share a feature: their + # counts are a sum and are not comparable with a guide ordered once. + for column in ("oligo", "n_oligos"): + # An alignment table written before these existed has neither, and + # the column is then written empty rather than the table refused. + out[column] = guides[column] if column in guides.columns else pd.NA element = chosen.reindex(out.index) kind = element.map(info["genomic_element"]) 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 88310d4..483ea98 100644 --- a/src/crispr_analysis_utils/skills/cau-guide-alignment/SKILL.md +++ b/src/crispr_analysis_utils/skills/cau-guide-alignment/SKILL.md @@ -70,6 +70,13 @@ tables, so it always describes the folder it sits in. See Three things trip people up: +- **`n_oligos` above one means those counts will be a sum.** Two library + entries of one spacer, typically a 20-mer and the same spacer with a 5' G, + are one feature in any reference built on the spacer, and no read can be + assigned to one or the other when the leader ends in a poly-G run: the + extra base is indistinguishable from the leader's own. `oligo` records what + was ordered, `spacer` is what the reads can resolve, and a guide flagged + this way is not comparable with one ordered once. - **`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_target_sites` is the one that matters for cutting; `n_matches` is the one that matters for 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 4d53e72..265a8f2 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 @@ -26,6 +26,8 @@ empty. The guide's own columns repeat on each of its rows. | `guide_id` | The sequence that was aligned, ending in the PAM pattern. Identifies the guide. | | `alias` | The name the input table gives it. | | `spacer` | The protospacer: `guide_id` without any base 5' of it or the PAM. | +| `oligo` | The sequence the library ordered, without the PAM: the spacer plus any base 5' of it. | +| `n_oligos` | How many entries share this spacer. Above one the counts are a sum, because one feature is all the reads can resolve. | | `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_target_sites` | How many of those have a PAM. | 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 596ae0a..dc5128a 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 @@ -18,6 +18,8 @@ guide and one spacer. | --- | --- | | `guide_id` | The guide's name, from the input table's own name for it. | | `spacer` | The targeting sequence, 5' to 3', in DNA bases. | +| `oligo` | The sequence ordered, without the PAM. | +| `n_oligos` | Entries sharing this spacer; above one the counts are a sum. | | `targeting` | TRUE only when the guide targets a relevant locus. | | `type` | `targeting`, `safe-targeting`, `non-targeting` or `positive control`. | | `guide_chr`, `guide_start`, `guide_end` | The protospacer, 0-based and half-open, **without** the PAM. | diff --git a/tests/test_guide_alignment_summary.py b/tests/test_guide_alignment_summary.py index 9992f8f..de143a3 100644 --- a/tests/test_guide_alignment_summary.py +++ b/tests/test_guide_alignment_summary.py @@ -59,6 +59,8 @@ "alias", "accession", "spacer", + "oligo", + "n_oligos", "site_class", "n_matches", "n_target_sites", @@ -72,7 +74,7 @@ "pam_matches", "mapped_sequence", ] -SITE_COLUMN_NAMES = COLUMN_NAMES[9:] +SITE_COLUMN_NAMES = COLUMN_NAMES[11:] def site(spacer, start, strand="+", mismatches=0, pam_class="NGG", pam="TGG"): @@ -345,12 +347,14 @@ def test_guides_tsv(tmp_path, guides, sites): assert lines[0].split("\t") == COLUMN_NAMES assert lines[1].split("\t") == ( [SPACER_A + "NGG", "uniq", guide_accession(SPACER_A), SPACER_A] + # The fixture carries SPACER_A twice, as uniq and dup. + + [SPACER_A, "2"] + ["unique", "2", "1", "forward", ""] + ["chr1", "100", "120", "+", "117", "TRUE", SPACER_A + "TGG"] ) # The guide's own columns repeat on its second row; only the site changes. - assert lines[2].split("\t")[:9] == lines[1].split("\t")[:9] - assert lines[2].split("\t")[9:] == ( + assert lines[2].split("\t")[:11] == lines[1].split("\t")[:11] + assert lines[2].split("\t")[11:] == ( ["chr1", "500", "520", "+", "517", "FALSE", SPACER_A + "TTT"] ) # One row per entry per site, so the 11 sites of the library are 11 rows. @@ -623,3 +627,32 @@ def test_guide_id_ends_in_the_pam_pattern_and_mapped_sequence_in_the_real_pam(tm SPACER_A + "GGG", "G" + SPACER_T + "TGG", ] + + +def test_n_oligos_marks_the_entries_whose_counts_will_be_a_sum(): + """Two entries of one spacer share a feature, so their reads land together. + + The pair is a 20-mer and the same spacer with a 5' G. No read can be + assigned to one or the other: the leader ends in a poly-G run, so the + extra base is indistinguishable from the leader's own. The counts are a + sum, and n_oligos is what says so. + """ + plain, with_g = SPACER_A, "G" + SPACER_A + guides = read_guides(pd.DataFrame([("plain", plain), ("with_g", with_g)])) + summaries = summarize( + guides, + {0: [site(plain, 100)], 1: [site(plain, 100)]}, + orientation={0: "forward", 1: "forward"}, + ) + assert [s.spacer for s in summaries] == [plain, plain] + assert [s.oligo for s in summaries] == [plain, with_g] + assert [s.n_oligos for s in summaries] == [2, 2] + # One accession, because one feature is what the reads can distinguish. + assert len({guide_accession(s.spacer) for s in summaries}) == 1 + + +def test_a_guide_ordered_once_reports_one_oligo(guides, sites): + summaries = summarized(guides, sites) + singles = [s for s in summaries if s.alias not in {"uniq", "dup"}] + assert {s.n_oligos for s in singles} == {1} + assert all(s.oligo == s.spacer for s in singles) diff --git a/tests/test_igvf_export.py b/tests/test_igvf_export.py index 00bfe1f..1826e04 100644 --- a/tests/test_igvf_export.py +++ b/tests/test_igvf_export.py @@ -442,8 +442,12 @@ def test_a_written_coordinate_has_no_decimal_point(table, tmp_path): def test_a_field_that_does_not_apply_is_an_empty_cell(table, tmp_path): write_table(tmp_path / "igvf.tsv", table) lines = (tmp_path / "igvf.tsv").read_text().splitlines() - empty = "\t" * (len(COLUMNS) - 4) - assert lines[-1] == f"non_targeting\t{spacer('non_targeting')}\tFALSE\t{NON}{empty}" + # guide_id, spacer, oligo, n_oligos, targeting, type: the fixture's + # alignment predates oligo and n_oligos, so both are empty. + empty = "\t" * (len(COLUMNS) - 6) + assert lines[-1] == ( + f"non_targeting\t{spacer('non_targeting')}\t\t\tFALSE\t{NON}{empty}" + ) # The reporter and the imperfect guide are unused, so every row ends empty. assert all(line.endswith("\t\t") for line in lines[1:]) @@ -453,7 +457,9 @@ def test_a_column_the_table_lacks_is_written_empty(tmp_path): write_table(tmp_path / "bare.tsv", bare) lines = (tmp_path / "bare.tsv").read_text().splitlines() assert lines[0].split("\t") == list(COLUMNS) - assert lines[1].split("\t") == ["a", "", "", NON] + [""] * (len(COLUMNS) - 4) + assert lines[1].split("\t") == ["a", "", "", "", "", NON] + [""] * ( + len(COLUMNS) - 6 + ) def test_columns_yaml_parses_and_documents_the_table(table, tmp_path): From 770ab8fee34d17f73e0f91ffbd13b5aca0c2beb5 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Thu, 15 Oct 2026 11:44:16 -0400 Subject: [PATCH 2/3] FIX the measured-gene lookup missing every element by one base chr:start-end is what a genome browser takes, and a browser is 1-based, so an element named that way starts one base after the BED row it came from. A pairs table that concatenated the BED's own columns is 0-based, and the exact string match between them found nothing: all 26 validated enhancers of the hiPSC-to-EC library carried no putative_target_genes, silently, until someone asked where the target genes had gone. The lookup now falls back to that one offset. Only that direction: the reverse does not arise from the convention and would risk matching the element next door. On the library, 180 of 340 enhancer guides gain a gene. The other 160 are 14 elements the BED marks e2g_validated but which have no row in the pairs table, so their genes are missing from the source rather than from the join. Co-Authored-By: Claude Opus 5 --- src/crispr_analysis_utils/igvf_export.py | 27 ++++++++++++++++++- .../skills/cau-element-mapping/SKILL.md | 6 +++++ tests/test_igvf_export.py | 17 ++++++++++++ 3 files changed, 49 insertions(+), 1 deletion(-) diff --git a/src/crispr_analysis_utils/igvf_export.py b/src/crispr_analysis_utils/igvf_export.py index c141b0a..932fcb9 100644 --- a/src/crispr_analysis_utils/igvf_export.py +++ b/src/crispr_analysis_utils/igvf_export.py @@ -14,6 +14,7 @@ import logging import os +import re from collections.abc import Mapping from dataclasses import dataclass, field from datetime import UTC, datetime @@ -221,9 +222,33 @@ def _element_kind(row, in_gene: set, in_exon: set, near: Mapping, window: int) - return DISTAL +_COORDS = re.compile(r"^(?P[^:]+):(?P\d+)-(?P\d+)$") + + +def _measured_genes(known: Mapping[str, list[str]], element_id: str) -> list[str]: + """The measured genes of an element, matching on coordinates when needed. + + ``chr:start-end`` is what a genome browser takes, and a browser is + 1-based, so an element named that way starts one base after the BED row + it came from. A pairs table that concatenated the BED's own columns is + 0-based instead, and the two spellings of one element differ by a single + base. All 26 validated enhancers of the hiPSC-to-EC library missed this + way. Only that one direction is tried: the reverse would not arise from + this convention, and would risk matching the element next door. + """ + hit = known.get(element_id, []) + if hit: + return hit + found = _COORDS.match(str(element_id)) + if not found: + return [] + chrom, start, end = found["chrom"], int(found["start"]), int(found["end"]) + return known.get(f"{chrom}:{start - 1}-{end}", []) + + def _genes_for(row, known: Mapping[str, list[str]]) -> list[str]: """Every gene associated with an element, measured ones first.""" - measured = known.get(row.element_id, []) + measured = _measured_genes(known, row.element_id) if measured: return list(dict.fromkeys(measured)) if row.genomic_element == PROMOTER and isinstance(row.nearest_ensg, str): 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 6a8d592..e723078 100644 --- a/src/crispr_analysis_utils/skills/cau-element-mapping/SKILL.md +++ b/src/crispr_analysis_utils/skills/cau-element-mapping/SKILL.md @@ -73,6 +73,12 @@ in full. - **`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. +- **That id is 1-based over a 0-based interval, and the two cannot be + concatenated.** `chr:start-end` is what a genome browser takes, so its + start is the BED's `Start` plus one. A table that pastes a BED's own + columns back together is 0-based and will miss every element by a single + base: all 26 validated enhancers of one library matched nothing that way, + silently, until someone asked where the target genes had gone. - **It keys back** to `guide_genome_annotation.tsv` on `guide_id` and the four site columns. diff --git a/tests/test_igvf_export.py b/tests/test_igvf_export.py index 1826e04..a916589 100644 --- a/tests/test_igvf_export.py +++ b/tests/test_igvf_export.py @@ -518,3 +518,20 @@ def test_a_targeting_guide_without_a_pam_is_counted(alignment, mapping, info): table = build_table(without, mapping, info) assert pd.isna(guide(table, "distal")["pam"]) assert counts(table, info)["missing_pam"] == 1 + + +def test_a_measured_pair_matches_across_a_one_base_start_difference(): + """chr:start-end is browser notation, so 1-based; a pairs table built + from the BED's own columns is 0-based.""" + from crispr_analysis_utils.igvf_export import _measured_genes + + known = {"chr11:35618822-35619124": ["ENSG00000179431"]} + # The element as a browser would name it, one base along. + assert _measured_genes(known, "chr11:35618823-35619124") == ["ENSG00000179431"] + # The reverse does not arise from this convention and is not tried, so a + # neighbouring element one base below cannot be matched by accident. + assert _measured_genes(known, "chr11:35618821-35619124") == [] + # An exact match still wins, and an unrelated element still finds nothing. + assert _measured_genes(known, "chr11:35618822-35619124") == ["ENSG00000179431"] + assert _measured_genes(known, "chr11:99999999-99999999") == [] + assert _measured_genes(known, "not_coordinates") == [] From 971722e0653a626366e5077b8c0873af3bf1c59d Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Thu, 15 Oct 2026 11:44:25 -0400 Subject: [PATCH 3/3] ADD the annotation tables to guide.uns, where they can be joined guide.var holds the IGVF columns, and the signed distances, promoter transcripts and other-gene-type overlaps were nowhere in the object. They cannot go in var either: guide_gene_annotation.tsv is one row per guide per site per gene, 45,501 for 43,766 guides, and collapsing it would force a rule on multi-site guides that nothing downstream asked for. So both tables are carried whole in mod['guide'].uns, by --guide-annotation and --element-mapping. Nothing is decided; the rows are there to be joined. Co-Authored-By: Claude Opus 5 --- .../cli/mudata_builder.py | 15 ++++++++++ src/crispr_analysis_utils/mudata_builder.py | 28 +++++++++++++++++++ 2 files changed, 43 insertions(+) diff --git a/src/crispr_analysis_utils/cli/mudata_builder.py b/src/crispr_analysis_utils/cli/mudata_builder.py index 6e105da..9304d3e 100644 --- a/src/crispr_analysis_utils/cli/mudata_builder.py +++ b/src/crispr_analysis_utils/cli/mudata_builder.py @@ -56,6 +56,17 @@ def add_arguments(parser: argparse.ArgumentParser) -> None: ), ) run_parser.add_argument("--genome-accession", help="The genome's IGVF portal id.") + run_parser.add_argument( + "--guide-annotation", + help=( + "guide_gene_annotation.tsv. Carried whole in guide.uns: it is " + "one row per guide per site per gene, so it cannot go in var." + ), + ) + run_parser.add_argument( + "--element-mapping", + help="guide_element_mapping.tsv. Carried whole in guide.uns.", + ) run_parser.add_argument( "--covariates", help=( @@ -81,6 +92,8 @@ def run(args: argparse.Namespace) -> int: gtf_accession=args.gtf_accession, genome_accession=args.genome_accession, covariates=args.covariates, + guide_annotation=args.guide_annotation, + element_mapping=args.element_mapping, ) except (BuildError, FileNotFoundError, OSError) as error: raise CommandError(str(error)) from error @@ -97,5 +110,7 @@ def run(args: argparse.Namespace) -> int: f"Unannotated: {found['genes_unannotated']} genes, " f"{found['guides_unannotated']} guides" ) + for name, rows in found.get("uns_tables", {}).items(): + print(f"guide.uns[{name!r}]: {rows} rows") print(f"Outputs in {Path(result.outdir)}") return 0 diff --git a/src/crispr_analysis_utils/mudata_builder.py b/src/crispr_analysis_utils/mudata_builder.py index 01063cf..0396cdb 100644 --- a/src/crispr_analysis_utils/mudata_builder.py +++ b/src/crispr_analysis_utils/mudata_builder.py @@ -278,6 +278,7 @@ def build( measurement_set: str, reference: dict | None = None, covariates: pd.DataFrame | None = None, + tables: dict[str, pd.DataFrame] | None = None, ): """Assemble the two matrices into one MuData. @@ -296,6 +297,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. + 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 + ``var`` without collapsing it; keeping it whole loses nothing and + decides nothing. Returns ------- @@ -337,6 +343,9 @@ def build( # MuData takes the union of the two indexes and records membership in # obsm, so a barcode seen in one modality only is kept and visible. + for name, table in (tables or {}).items(): + guide.uns[name] = table + data = md.MuData({"gene": rna, "guide": guide}) data.uns["reference"] = dict(reference or {}) data.uns["measurement_set"] = measurement_set @@ -382,6 +391,8 @@ def run( gtf_accession: str | None = None, genome_accession: str | None = None, covariates: str | os.PathLike[str] | None = None, + guide_annotation: str | os.PathLike[str] | None = None, + element_mapping: str | os.PathLike[str] | None = None, ) -> RunResult: """Assemble one channel's matrices into a MuData, and write it. @@ -405,6 +416,10 @@ def run( covariates A table of per-cell covariates, indexed by barcode in its first column. Carried into the root ``obs`` unchanged. + guide_annotation, element_mapping + ``guide_gene_annotation.tsv`` and ``guide_element_mapping.tsv``. + Carried whole in ``mod['guide'].uns``, since neither is one row per + guide. ``var`` keeps the IGVF columns alone. Returns ------- @@ -467,6 +482,17 @@ def run( state["inputs"]["covariates"] = str(covariates) state["inputs"]["covariates_md5"] = file_md5(covariates) + tables = {} + for name, path in ( + ("annotation", guide_annotation), + ("element_mapping", element_mapping), + ): + if path is None: + continue + tables[name] = pd.read_csv(path, sep="\t", low_memory=False) + state["inputs"][name] = str(path) + state["inputs"][f"{name}_md5"] = file_md5(path) + reference = reference_block( loaded, gtf_accession=gtf_accession, @@ -483,6 +509,7 @@ def run( measurement_set=measurement_set, reference=reference, covariates=covariate_table, + tables=tables, ) data.write(outdir / RESULTS_DIR / TABLE) @@ -499,6 +526,7 @@ def run( "guides": int(data["guide"].n_vars), "genes_unannotated": int(data["gene"].var["symbol"].isna().sum()), "guides_unannotated": int(data["guide"].var["spacer"].isna().sum()), + "uns_tables": {k: int(len(v)) for k, v in tables.items()}, } state["complete"] = True state["finished"] = _now()