From 7d9cbe1b69c0eeaf421796f1ba25c3292405d9ce Mon Sep 17 00:00:00 2001 From: hq <29302823+lhqing@users.noreply.github.com> Date: Sat, 22 Aug 2026 12:58:11 -0400 Subject: [PATCH 1/4] the genome pin moves forward to the commit that added the gene list MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The lock held liulab-genome at ce23e02, which predates `Genome.gene_list` — the accessor that hands back the gene ids under a biotype. The pin was fine for as long as nothing here asked an annotation which of its genes were which; the per-cell rRNA metric that follows is the first caller that does, and against the old pin it cannot import. Nothing else in the lock moves: the only lines that differ are the six per-platform pins of that one git dependency and the version string derived from it. Co-Authored-By: Claude Opus 5 (1M context) --- pixi.lock | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/pixi.lock b/pixi.lock index 1ba07d5..2a596d5 100644 --- a/pixi.lock +++ b/pixi.lock @@ -338,7 +338,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/zipp-4.1.0-pyhcf101f3_0.conda - pypi: ./ - pypi: git+https://github.com/liuhlab/liulab-data.git#97d0ad9ee3c4645ef758336d71bb10b03eda4ca3 - - pypi: git+https://github.com/liuhlab/liulab-genome.git#ce23e02fdf3100ab05a64b64e639085469fd8a79 + - pypi: git+https://github.com/liuhlab/liulab-genome.git#30e302516c653189108e3baffb8b83ddea919082 - pypi: https://files.pythonhosted.org/packages/59/a1/18fad0ca587eef5ba25a207db5f9730058365825bfdd1fc6650ab2b25609/py2bit-1.0.1-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl - pypi: https://files.pythonhosted.org/packages/c4/ea/066ce356c5df3c2d42b72801f768d41bf691f58f6d2a90f6334fbed39785/pysam-0.24.0-cp313-cp313-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl osx-64: @@ -628,7 +628,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-64/zstd-1.5.7-h3eecb57_6.conda - pypi: ./ - pypi: git+https://github.com/liuhlab/liulab-data.git#97d0ad9ee3c4645ef758336d71bb10b03eda4ca3 - - pypi: git+https://github.com/liuhlab/liulab-genome.git#ce23e02fdf3100ab05a64b64e639085469fd8a79 + - pypi: git+https://github.com/liuhlab/liulab-genome.git#30e302516c653189108e3baffb8b83ddea919082 - pypi: https://files.pythonhosted.org/packages/46/3b/9dedd2e35cebcd4eb3539514b5ac674fd45654404839e6f9eff6e25d67c7/py2bit-1.0.1.tar.gz - pypi: https://files.pythonhosted.org/packages/86/79/2f5151ac001d8c74fb047036bfea9e4e897939e6587d3c4d512e46c450b1/pysam-0.24.0-cp313-cp313-macosx_10_13_x86_64.whl osx-arm64: @@ -917,7 +917,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-arm64/zstd-1.5.7-hbf9d68e_6.conda - pypi: ./ - pypi: git+https://github.com/liuhlab/liulab-data.git#97d0ad9ee3c4645ef758336d71bb10b03eda4ca3 - - pypi: git+https://github.com/liuhlab/liulab-genome.git#ce23e02fdf3100ab05a64b64e639085469fd8a79 + - pypi: git+https://github.com/liuhlab/liulab-genome.git#30e302516c653189108e3baffb8b83ddea919082 - pypi: https://files.pythonhosted.org/packages/46/3b/9dedd2e35cebcd4eb3539514b5ac674fd45654404839e6f9eff6e25d67c7/py2bit-1.0.1.tar.gz - pypi: https://files.pythonhosted.org/packages/d1/2c/fd59b47677a1df3efa64172dcd9b99fa7db437de8c663f08120ebd4db835/pysam-0.24.0-cp313-cp313-macosx_11_0_arm64.whl docs: @@ -1214,7 +1214,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/zipp-4.1.0-pyhcf101f3_0.conda - pypi: ./ - pypi: git+https://github.com/liuhlab/liulab-data.git#97d0ad9ee3c4645ef758336d71bb10b03eda4ca3 - - pypi: git+https://github.com/liuhlab/liulab-genome.git#ce23e02fdf3100ab05a64b64e639085469fd8a79 + - pypi: git+https://github.com/liuhlab/liulab-genome.git#30e302516c653189108e3baffb8b83ddea919082 - pypi: https://files.pythonhosted.org/packages/59/a1/18fad0ca587eef5ba25a207db5f9730058365825bfdd1fc6650ab2b25609/py2bit-1.0.1-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl - pypi: https://files.pythonhosted.org/packages/c4/ea/066ce356c5df3c2d42b72801f768d41bf691f58f6d2a90f6334fbed39785/pysam-0.24.0-cp313-cp313-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl osx-64: @@ -1477,7 +1477,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-64/zstd-1.5.7-h3eecb57_6.conda - pypi: ./ - pypi: git+https://github.com/liuhlab/liulab-data.git#97d0ad9ee3c4645ef758336d71bb10b03eda4ca3 - - pypi: git+https://github.com/liuhlab/liulab-genome.git#ce23e02fdf3100ab05a64b64e639085469fd8a79 + - pypi: git+https://github.com/liuhlab/liulab-genome.git#30e302516c653189108e3baffb8b83ddea919082 - pypi: https://files.pythonhosted.org/packages/46/3b/9dedd2e35cebcd4eb3539514b5ac674fd45654404839e6f9eff6e25d67c7/py2bit-1.0.1.tar.gz - pypi: https://files.pythonhosted.org/packages/86/79/2f5151ac001d8c74fb047036bfea9e4e897939e6587d3c4d512e46c450b1/pysam-0.24.0-cp313-cp313-macosx_10_13_x86_64.whl osx-arm64: @@ -1739,7 +1739,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-arm64/zstd-1.5.7-hbf9d68e_6.conda - pypi: ./ - pypi: git+https://github.com/liuhlab/liulab-data.git#97d0ad9ee3c4645ef758336d71bb10b03eda4ca3 - - pypi: git+https://github.com/liuhlab/liulab-genome.git#ce23e02fdf3100ab05a64b64e639085469fd8a79 + - pypi: git+https://github.com/liuhlab/liulab-genome.git#30e302516c653189108e3baffb8b83ddea919082 - pypi: https://files.pythonhosted.org/packages/46/3b/9dedd2e35cebcd4eb3539514b5ac674fd45654404839e6f9eff6e25d67c7/py2bit-1.0.1.tar.gz - pypi: https://files.pythonhosted.org/packages/d1/2c/fd59b47677a1df3efa64172dcd9b99fa7db437de8c663f08120ebd4db835/pysam-0.24.0-cp313-cp313-macosx_11_0_arm64.whl test: @@ -2029,7 +2029,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/zipp-4.1.0-pyhcf101f3_0.conda - pypi: ./ - pypi: git+https://github.com/liuhlab/liulab-data.git#97d0ad9ee3c4645ef758336d71bb10b03eda4ca3 - - pypi: git+https://github.com/liuhlab/liulab-genome.git#ce23e02fdf3100ab05a64b64e639085469fd8a79 + - pypi: git+https://github.com/liuhlab/liulab-genome.git#30e302516c653189108e3baffb8b83ddea919082 - pypi: https://files.pythonhosted.org/packages/59/a1/18fad0ca587eef5ba25a207db5f9730058365825bfdd1fc6650ab2b25609/py2bit-1.0.1-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl - pypi: https://files.pythonhosted.org/packages/c4/ea/066ce356c5df3c2d42b72801f768d41bf691f58f6d2a90f6334fbed39785/pysam-0.24.0-cp313-cp313-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl osx-64: @@ -2285,7 +2285,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-64/zstd-1.5.7-h3eecb57_6.conda - pypi: ./ - pypi: git+https://github.com/liuhlab/liulab-data.git#97d0ad9ee3c4645ef758336d71bb10b03eda4ca3 - - pypi: git+https://github.com/liuhlab/liulab-genome.git#ce23e02fdf3100ab05a64b64e639085469fd8a79 + - pypi: git+https://github.com/liuhlab/liulab-genome.git#30e302516c653189108e3baffb8b83ddea919082 - pypi: https://files.pythonhosted.org/packages/46/3b/9dedd2e35cebcd4eb3539514b5ac674fd45654404839e6f9eff6e25d67c7/py2bit-1.0.1.tar.gz - pypi: https://files.pythonhosted.org/packages/86/79/2f5151ac001d8c74fb047036bfea9e4e897939e6587d3c4d512e46c450b1/pysam-0.24.0-cp313-cp313-macosx_10_13_x86_64.whl osx-arm64: @@ -2540,7 +2540,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-arm64/zstd-1.5.7-hbf9d68e_6.conda - pypi: ./ - pypi: git+https://github.com/liuhlab/liulab-data.git#97d0ad9ee3c4645ef758336d71bb10b03eda4ca3 - - pypi: git+https://github.com/liuhlab/liulab-genome.git#ce23e02fdf3100ab05a64b64e639085469fd8a79 + - pypi: git+https://github.com/liuhlab/liulab-genome.git#30e302516c653189108e3baffb8b83ddea919082 - pypi: https://files.pythonhosted.org/packages/46/3b/9dedd2e35cebcd4eb3539514b5ac674fd45654404839e6f9eff6e25d67c7/py2bit-1.0.1.tar.gz - pypi: https://files.pythonhosted.org/packages/d1/2c/fd59b47677a1df3efa64172dcd9b99fa7db437de8c663f08120ebd4db835/pysam-0.24.0-cp313-cp313-macosx_11_0_arm64.whl test-star: @@ -12368,9 +12368,9 @@ packages: - snakemake>=8 - typer>=0.12 requires_python: '>=3.12' -- pypi: git+https://github.com/liuhlab/liulab-genome.git#ce23e02fdf3100ab05a64b64e639085469fd8a79 +- pypi: git+https://github.com/liuhlab/liulab-genome.git#30e302516c653189108e3baffb8b83ddea919082 name: liulab-genome - version: 2026.8.1.dev1+gce23e02fd + version: 2026.8.1.dev2+g30e302516 requires_dist: - numpy>=2 - pandas>=2.1,<3 From d24362d8300704e1e956f03c069e22fb1ed354c6 Mon Sep 17 00:00:00 2001 From: hq <29302823+lhqing@users.noreply.github.com> Date: Sat, 22 Aug 2026 13:04:51 -0400 Subject: [PATCH 2/4] the counter says how much of each cell's yield was ribosomal MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A ribosomal fragment maps uniquely, reaches exactly one gene, and is credited to it like any other — nothing in this counter filters one out and nothing downstream of it does either. So a plate object can report an ordinary-looking molecule total over an ordinary-looking number of genes while almost none of it is the transcriptome anybody asked for, and until now there was nothing on the object that would say so. Two `obs` columns say it: `n_umis_rrna`, the cell's deduplicated molecules falling in ribosomal genes, and `rrna_fraction`, that count over `n_umis`. Both are counted over `umi_combined` and never over `X`. That is the same population `n_umis`, `genes_detected` and the saturation are all read off, so the share and the total it is a share OF are one number rather than two derivations nothing forces to agree; the primary matrix is exonic, and a numerator taken off it would report a different share of a different population. A cell with no molecules at all gets `nan` and never `0.0` — the denominator is the arithmetic with no answer, exactly as it is for saturation, and the reader already reads that back as absence. The ids arrive as a parameter rather than being looked up here. This module takes a resolved annotation database and never an assembly id, precisely so the `liulab-genome` import can live in the CLI verb and this stay strictly typed and testable against a synthetic annotation; the ribosomal gene set is the same kind of fact and enters the same way. An id that is not on the gene axis is dropped and the rest are counted, because the ribosomal set is curated per assembly while the axis is whatever GTF was registered, and the two overlap partially all the time. What an EMPTY overlap does is the decision this whole measurement rests on: neither column is written at all. A plate whose every cell reads `0.0% rRNA` is indistinguishable from a plate that really is ribosome-free, and it is the one answer nobody may act on — so where no gene was named, or none of the named ones is on this annotation, the object carries no ribosomal column and the report omits the metric, which is what an unmeasured thing looks like everywhere else on this path. It is absence rather than a refusal because the plate itself counted fine: what is missing is a QC number about it, and refusing a whole plate over a gene list would cost every cell's matrices for a column. The page gets the share and not the count. A raw ribosomal molecule total is not comparable across cells whose depths differ by three orders of magnitude, which is the trade the four fates already make against `n_fragments`; the count stays on the object for whoever wants to recompute. It is ungraded, like every other column this reader builds — how much rRNA a library should carry is a property of the prep and the organism, nobody has measured a bar for it, and a bar invented at review is worse than a number a reader compares across the plate themselves. Nothing else moves: no matrix, no `uns` note, no `var` column, and no gene the counter used to keep is dropped now. `WORKFLOW_VERSION` moves to 2026.8.24 because the ARTIFACT changed and not because a command line did. The rule's params, its declared outputs and every matrix it writes are byte-identical; what moves is the shape of the object that comes out, and an h5ad written before this beside one written after would otherwise both claim the same workflow produced them. Closes #471 Co-Authored-By: Claude Opus 5 (1M context) --- src/seqforge/workflows/__init__.py | 14 ++- src/seqforge/workflows/umite/count.py | 121 +++++++++++++++++++++++++- tests/test_workflows.py | 105 ++++++++++++++++++++++ 3 files changed, 235 insertions(+), 5 deletions(-) diff --git a/src/seqforge/workflows/__init__.py b/src/seqforge/workflows/__init__.py index 03b3bb8..176b84c 100644 --- a/src/seqforge/workflows/__init__.py +++ b/src/seqforge/workflows/__init__.py @@ -35,6 +35,18 @@ from ..kb.schema import Spec #: CalVer YYYY.M.PATCH; bump when any shipped module's rules/params change. +#: 2026.8.24 — `rule umi_count` says how much of each cell's yield was ribosomal (#471). Two more +#: `obs` columns on the plate object — `n_umis_rrna` and `rrna_fraction`, both counted over +#: `umi_combined`, the same population `n_umis` and the saturation are read off — and the page carries +#: the share beside them, ungraded. The gene set is liulab-genome's curated `rRNA` category for the +#: annotation the plate was actually counted against, so nothing here owns a list of ribosomal genes; +#: an annotation that cannot answer, or one whose axis shares no id with the answer, yields NEITHER +#: column rather than a plate of `0.0%`, which is the exact claim the measurement exists to prevent. +#: **The bump is owed by the ARTIFACT and not by a moved command line**, the same as 2026.8.18: the +#: rule's params, its declared outputs and every matrix it writes are byte-identical — these genes +#: were already counted as expression and still are, and nothing was filtered — and what moves is the +#: shape of the object that comes out, which an h5ad written before this could not otherwise be told +#: apart from one written after. #: 2026.8.23 — `rule genome_index` is DELETED from all five modules and the lookup it performed is a #: `params:` callable in each (#478, from #475). It owned `results/index/`, a path every #: concurrent instance over one results directory shares, and snakemake removes an output before @@ -611,7 +623,7 @@ #: dereferenced and never declared. The contract was wrong, not the module. #: 2026.7.1 — star.smk hardcodes --outSAMtype (it is a module detail, and starsolo.smk always #: hardcoded it); required_config gains primary_feature and drops bulk.outSAMtype. -WORKFLOW_VERSION = "2026.8.23" +WORKFLOW_VERSION = "2026.8.24" _MODULE_DIR = Path(__file__).parent diff --git a/src/seqforge/workflows/umite/count.py b/src/seqforge/workflows/umite/count.py index b133249..dc52bdd 100644 --- a/src/seqforge/workflows/umite/count.py +++ b/src/seqforge/workflows/umite/count.py @@ -229,6 +229,28 @@ class UmiCountError(RuntimeError): N_UMIS = "n_umis" GENES_DETECTED = "genes_detected" +#: How much of that yield was ribosomal: the cell's deduplicated molecules falling in rRNA genes, and +#: that count over :data:`N_UMIS`. Counted over ``umi_combined`` for the reason the two columns above +#: are — the total this is a share OF is that layer's, so numerator and denominator are one population +#: rather than two derivations that agree by coincidence. +#: +#: **These genes are counted as EXPRESSION in every matrix here, which is the whole reason the number +#: is worth writing down.** A ribosomal fragment maps uniquely, reaches exactly one gene and is +#: credited to it like any other; nothing in this counter filters one out, and nothing downstream of +#: it does either. So a library that is mostly rRNA reports an ordinary-looking molecule total over an +#: ordinary-looking number of genes while almost none of it is the transcriptome anybody asked for, +#: and the only way to see that from the object is a column that says how much of it is. This is a QC +#: number a consumer reads BEFORE deciding whether to drop those genes — the decision stays theirs, +#: and the matrices are unchanged by the measurement. +#: +#: A cell with no molecules at all gets ``nan`` and never a zero, the rule :data:`SATURATION` already +#: follows: the denominator is the arithmetic with no answer, and a rendered ``0.0%`` is a share a +#: reader acts on. And where no ribosomal gene was named, or none of the named ones is on the +#: annotation being counted against, NEITHER column is written — a plate-wide ``0.0%`` is precisely +#: the claim this measurement exists to keep anybody from making, and absence is the honest answer. +N_UMIS_RRNA = "n_umis_rrna" +RRNA_FRACTION = "rrna_fraction" + # -------------------------------------------------------------------------------------------- # The annotation, read once @@ -1027,8 +1049,42 @@ def _hit_counts(counted: Sequence[_Counted]) -> np.ndarray: ) +def _rrna_umis( + entries: Sequence[Mapping[str, Mapping[int, int]]], + annotation: Annotation, + rrna_gene_ids: Sequence[str] | None, +) -> list[int] | None: + """Each cell's combined molecules in the named genes, or ``None`` when there are none to name. + + The ids are matched against the gene axis of the annotation the plate was actually counted + against, and an id that is not on it is dropped rather than refused: the ribosomal set is a + property of the assembly's curation and the axis is a property of the GTF that was registered, so + the two overlap partially all the time and a partial overlap is still a measurement. + + An EMPTY overlap is not, and it is the one case that gets its own answer. Returning zeros there + would put a column on the object saying every cell of this plate is 0.0% ribosomal, which is + indistinguishable from a measurement and is the exact wrong claim — see :data:`RRNA_FRACTION`. + ``None`` reaches :func:`count_plate` as "write neither column", and a reader then finds the metric + absent, which is what an unmeasured thing looks like everywhere else on this path. + + The buckets are the ones the cells were already counted into, summed per cell over the matched + indices — no second pass over the matrices and no second traversal of a BAM. Genes are walked in + sorted order, like everything else here, so a plate summed twice adds in one order. + """ + if not rrna_gene_ids: + return None + axis = {gene_id: index for index, gene_id in enumerate(annotation.gene_ids)} + matched = sorted({axis[gene_id] for gene_id in rrna_gene_ids if gene_id in axis}) + if not matched: + return None + return [sum(entry["umi_combined"].get(gene, 0) for gene in matched) for entry in entries] + + def count_plate( - cells: Sequence[tuple[str, Path]], annotation: Annotation, workers: int = 1 + cells: Sequence[tuple[str, Path]], + annotation: Annotation, + workers: int = 1, + rrna_gene_ids: Sequence[str] | None = None, ) -> anndata.AnnData: """Every cell of a plate -> one AnnData, rows in the order the cells were given. @@ -1055,6 +1111,15 @@ def count_plate( the first twenty thousand fragments of the first cell and flat from there — so the width the memory arithmetic can afford is far wider than any node's core count. Asking for more workers than there are cells simply gets one per cell. + + ``rrna_gene_ids`` is which of this annotation's genes are ribosomal, and it arrives as a list of + ids rather than as a category name because the database that knows the categories is + ``liulab-genome`` and the import that asks it lives in the CLI verb — the same line this module + already holds for the annotation itself, and what keeps this function strictly typed and testable + against a synthetic annotation. Nothing is counted twice for it: the ids select columns of the + combined UMI buckets these cells were already counted into. Omit them, or hand over ids none of + which is on the gene axis, and the object carries no ribosomal column at all rather than a zero + one — see :data:`RRNA_FRACTION`. """ import anndata as ad @@ -1087,6 +1152,17 @@ def count_plate( [np.nan if cell.saturation is None else cell.saturation for cell in counted], dtype=np.float64, ) + ribosomal = _rrna_umis(entries, annotation, rrna_gene_ids) + if ribosomal is not None: + adata.obs[N_UMIS_RRNA] = np.array(ribosomal, dtype=np.int32) + # A cell with no molecules has no share of them, by the same rule as the ratio above it. + adata.obs[RRNA_FRACTION] = np.array( + [ + np.nan if cell.n_umis == 0 else molecules / cell.n_umis + for molecules, cell in zip(ribosomal, counted, strict=True) + ], + dtype=np.float64, + ) adata.obsm[MULTIMAPPING_HITS] = _hit_counts(counted) adata.uns["primary_matrix"] = PRIMARY_MATRIX adata.uns["multimapping_caveat"] = MULTIMAPPING_CAVEAT @@ -1094,7 +1170,11 @@ def count_plate( def write_umi_counts( - cells: Sequence[tuple[str, Path]], annotation_db: Path, out: Path, workers: int = 1 + cells: Sequence[tuple[str, Path]], + annotation_db: Path, + out: Path, + workers: int = 1, + rrna_gene_ids: Sequence[str] | None = None, ) -> Path: """The one entry point: N per-cell BAMs + the built annotation -> one ``.h5ad``. Returns ``out``. @@ -1105,9 +1185,14 @@ def write_umi_counts( **The annotation is read BEFORE the fan-out and never inside it**, which is what makes the width free: one gffutils read for the plate, and every worker forked from the process that already holds the result. + + ``rrna_gene_ids`` passes straight through to :func:`count_plate`, which is where what it costs + and what an empty overlap means are written down. It is a list of gene ids here for the same + reason ``annotation_db`` is a resolved path: whoever asked ``liulab-genome`` the question did it + on the far side of this signature. """ annotation = read_annotation(annotation_db) - adata = count_plate(cells, annotation, workers) + adata = count_plate(cells, annotation, workers, rrna_gene_ids) out.parent.mkdir(parents=True, exist_ok=True) adata.write_h5ad(out) return out @@ -1277,6 +1362,24 @@ def fate_metrics(row: Mapping[str, float], sample: str) -> SampleStats: built.append( sequencing_saturation(None if saturated is None or isnan(saturated) else saturated) ) + # The ribosomal share is a SHARE and never the count beside it: `n_umis_rrna` stays on the object + # for whoever wants to recompute, and the page carries the one figure that is comparable across + # cells whose depths differ by three orders of magnitude -- the same trade the fates make. It is + # ungraded on the argument every column here is ungraded on: how much rRNA a library SHOULD carry + # is a property of the prep and the organism, nobody has measured a bar for it, and an invented + # one is worse than a number a reader compares across the plate themselves. + ribosomal = row.get(RRNA_FRACTION) + built.append( + fraction( + RRNA_FRACTION, + "rRNA", + None if ribosomal is None or isnan(ribosomal) else ribosomal, + group="counts", + hint="Share of this cell's combined molecules that fell in ribosomal genes. Those genes " + "are counted as expression in the matrices beside this — nothing here removes them — so " + "read it before deciding what to do about them.", + ) + ) return SampleStats(sample_id=sample, metrics=[m for m in built if m is not None]) @@ -1295,7 +1398,15 @@ def _obs_columns(adata: anndata.AnnData) -> dict[str, list[float]]: frame: Any = adata.obs return { column: [float(value) for value in frame[column]] - for column in (*FATES, N_FRAGMENTS, SATURATION, N_UMIS, GENES_DETECTED) + for column in ( + *FATES, + N_FRAGMENTS, + SATURATION, + N_UMIS, + GENES_DETECTED, + N_UMIS_RRNA, + RRNA_FRACTION, + ) if column in frame.columns } @@ -1348,7 +1459,9 @@ def read_plate_stats(path: Path, samples: Sequence[str]) -> dict[str, SampleStat "MULTIMAPPING_LAYER", "N_FRAGMENTS", "N_UMIS", + "N_UMIS_RRNA", "PRIMARY_MATRIX", + "RRNA_FRACTION", "SATURATION", "UMI_TAG", "Annotation", diff --git a/tests/test_workflows.py b/tests/test_workflows.py index 3d81845..1e3a55f 100644 --- a/tests/test_workflows.py +++ b/tests/test_workflows.py @@ -174,7 +174,9 @@ MULTIMAPPING_LAYER, N_FRAGMENTS, N_UMIS, + N_UMIS_RRNA, PRIMARY_MATRIX, + RRNA_FRACTION, SATURATION, UmiCountError, _step_index, @@ -6014,6 +6016,92 @@ def test_saturation_is_the_molecules_over_the_gene_assigned_fragments_that_carri assert math.isnan(float(_frame(empty.obs).loc["untagged", SATURATION])) +def test_the_ribosomal_columns_are_the_named_genes_molecules_over_the_cells_own( + tmp_path: Path, +) -> None: + """A share of the yield, over genes this counter deliberately goes on counting as expression. + + A ribosomal fragment maps uniquely and reaches exactly one gene, so nothing here excludes it and + every matrix carries it like any other transcript. That is what makes the share worth writing + down rather than acting on: a cell can report an ordinary molecule total over an ordinary number + of genes with almost none of it the transcriptome anybody wanted, and this column is the only + thing on the object that would say so. + + Every figure is read off `_PLATE` by hand. `cell_a`'s combined UMI matrix holds GENE_A twice and + GENE_B once, so naming GENE_A ribosomal makes two of that cell's three molecules ribosomal and + none of `cell_b`'s one. The denominator is `n_umis` and the numerator comes off `umi_combined` — + the same population the molecule total and the saturation beside it are both read off, so all + three are readings of one number rather than three derivations nothing forces to agree. + """ + db, cells = _plate(tmp_path) + annotation = read_annotation(db) + + # A partial match proceeds over what did match, silently: the ribosomal set is curated per + # assembly and the gene axis is whatever GTF was registered, so the two overlap partially all the + # time and a partial overlap is still a measurement. + adata = count_plate(cells, annotation, rrna_gene_ids=["GENE_A", "GENE_NOBODY_DECLARED"]) + obs = _frame(adata.obs) + + assert int(obs.loc["cell_a", N_UMIS_RRNA]) == 2 + assert float(obs.loc["cell_a", RRNA_FRACTION]) == pytest.approx(2 / 3) + assert int(obs.loc["cell_a", N_UMIS]) == 3 # the denominator, on the object beside it + assert int(obs.loc["cell_b", N_UMIS_RRNA]) == 0 + assert float(obs.loc["cell_b", RRNA_FRACTION]) == 0.0 # measured, and a real zero + + # Over `umi_combined` and never over `X`, on the one cell that can tell the two apart: its GENE_A + # corrects to two molecules exonically and to ONE in the union, so a numerator taken off the + # primary matrix would report twice the cell's whole molecule total — a share above 1. + split = count_plate( + [("one_cell", _synthetic_bam(tmp_path / "split.bam", _SPLIT_NEIGHBOURS))], + annotation, + rrna_gene_ids=["GENE_A"], + ) + assert _row(split, "one_cell", "GENE_A") == 2 + assert int(_frame(split.obs).loc["one_cell", N_UMIS_RRNA]) == 1 + assert float(_frame(split.obs).loc["one_cell", RRNA_FRACTION]) == 1.0 + + # A cell with no molecules at all has no share of them: `nan` and never a zero, the rule the + # saturation above already follows and the page already reads back as absence. + untagged = _synthetic_bam( + tmp_path / "untagged.bam", (_Fragment("read_only", "chr1", 120, 180),) + ) + nothing = _frame( + count_plate([("untagged", untagged)], annotation, rrna_gene_ids=["GENE_A"]).obs + ) + assert int(nothing.loc["untagged", N_UMIS_RRNA]) == 0 + assert math.isnan(float(nothing.loc["untagged", RRNA_FRACTION])) + + +@pytest.mark.parametrize( + "rrna_gene_ids", + [None, (), ("GENE_NOBODY_DECLARED", "GENE_ALSO_ABSENT")], + ids=["unasked", "empty", "off-axis"], +) +def test_a_plate_with_no_ribosomal_gene_to_count_gets_neither_column_rather_than_a_zero( + tmp_path: Path, rrna_gene_ids: tuple[str, ...] | None +) -> None: + """The anti-zero clause, which is the reason this measurement exists at all. + + A plate whose every cell reads `0.0% rRNA` is indistinguishable from a plate that really is + ribosome-free, and it is the one answer nobody may act on. So an annotation nobody named genes + for, and one whose gene axis shares no id with the ones named — the assembly was curated + elsewhere, or the registered GTF uses another id namespace — write NEITHER column, and the reader + then omits the metric exactly as it does for a column an older object never carried. + + Absence rather than a refusal, because the plate itself counted fine: what is missing is a QC + number about it, and a fan-in that refused a whole plate over a gene list would cost every cell's + matrices for a column. + """ + db, cells = _plate(tmp_path) + + adata = count_plate(cells, read_annotation(db), rrna_gene_ids=rrna_gene_ids) + + assert not {N_UMIS_RRNA, RRNA_FRACTION} & set(adata.obs.columns) + assert ( + N_UMIS in adata.obs.columns + ) # ...and the yield it would have been a share of is untouched + + def test_counting_the_same_plate_twice_gives_a_byte_identical_h5ad(tmp_path: Path) -> None: """Determinism, asserted on the artifact rather than on the absence of `random` in the source. @@ -7551,6 +7639,23 @@ def test_a_cell_that_counted_nothing_has_no_rates_rather_than_four_zeroes() -> N SATURATION: 0.0 } + # The ribosomal share reads back under both rules at once. An object written against an + # annotation that could not name a ribosomal gene carries no column, and a cell with no molecules + # carries `nan` — and a page that rendered either as `0.0%` would be claiming a clean library. + assert RRNA_FRACTION not in {m.key for m in counted.metrics} + assert RRNA_FRACTION not in { + m.key for m in fate_metrics({RRNA_FRACTION: float("nan")}, "cell_v").metrics + } + measured = _by_key(fate_metrics({N_UMIS_RRNA: 1, RRNA_FRACTION: 0.25}, "cell_u")) + assert measured[RRNA_FRACTION].value == 0.25 + # Ungraded, like every other column this reader builds: nobody has measured how much rRNA a + # library SHOULD carry, and a bar invented at review is worse than a number read across cells. + assert measured[RRNA_FRACTION].level == "none" + # The COUNT stays on the object and gets no column of its own — a raw ribosomal molecule total is + # not comparable across cells whose depths differ by three orders of magnitude, which is the same + # trade the fates make against `n_fragments`. + assert N_UMIS_RRNA not in measured + # ---- the plate-assay UMI extractor --------------------------------------------------------------- # From 16eda8225f6d9606ef74ef94b0664c3d4cd6a261 Mon Sep 17 00:00:00 2001 From: hq <29302823+lhqing@users.noreply.github.com> Date: Sat, 22 Aug 2026 13:13:05 -0400 Subject: [PATCH 3/4] the verb asks the annotation which of its genes are ribosomal MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The counter has taken a ribosomal gene set since the commit before this one, and nothing handed it one, so no plate ever got the two `obs` columns it writes from them. `io umi-count` supplies them now: the annotation it has just resolved is asked for its rRNA gene list, and the ids go straight into `write_umi_counts`. The verb is where that question belongs for the same reason the annotation database is resolved here — the package that knows which genes are which is the one this file already imports, and the counting module stays strictly typed and testable against a synthetic annotation with no genome store in sight. The lookup sits AFTER the `--component` arm has rewritten the pair, which is load-bearing rather than incidental. Asked of a Chimera, `gene_list("rRNA")` concatenates the worm's ids and the bacterium's across its sources, and a fraction whose numerator pools two organisms is a share of a population nobody asked about. The counter already runs once per Component, so asking the component's own genome about the component's own registered annotation makes that impossible by construction — and seqforge handles no `.sources` at all. The two ways an annotation can decline the question are caught BY NAME: `NoGeneCategoriesError`, where no curated list ships for it, and `GeneCategoryNotDeclaredError`, where one ships and does not name ribosomes. Either costs the plate that one column and never a zero one, silently and with nothing on stdout, which is what an unmeasured thing looks like everywhere else on this path. Deliberately NOT `LookupError`: an unregistered annotation raises an `AnnotationNotRegisteredError`, which is a `KeyError` and therefore also a `LookupError`, and catching the base class would turn that refusal into a silently missing column. A curated file that is present and broken, or curated against another assembly, raises a `ValueError` and falls through to the exit that was already there — a defect in a shipped data file must not become an absent metric either. Both existing `io umi-count` tests stub `Genome`, so both stubs grow a `gene_list`; the first is now parametrized over what the annotation answers, with the two absences as rows, and reads the columns back off the written object rather than off the seam — a plate whose annotation could not answer carries neither. The `--component` test pins who was asked, since a lookup moved one line earlier would pool two organisms and produce a number that looks fine. Closes #472 Co-Authored-By: Claude Opus 5 (1M context) --- src/seqforge/cli/io.py | 37 ++++++++++++++-- tests/test_cli.py | 98 ++++++++++++++++++++++++++++++++++++------ 2 files changed, 118 insertions(+), 17 deletions(-) diff --git a/src/seqforge/cli/io.py b/src/seqforge/cli/io.py index 410f6b7..77a2bfe 100644 --- a/src/seqforge/cli/io.py +++ b/src/seqforge/cli/io.py @@ -847,6 +847,18 @@ def io_umi_count( rendered command line saying two contradictory things about which GTF was used, in a repo whose wiring gate reads rendered commands. + **Which of that annotation's genes are ribosomal is looked up here too, and deliberately after + the pair above is settled.** The counter is handed gene ids rather than a category name, for the + same reason it is handed a resolved database: the package that knows which genes are which is the + one this verb already imports, and the module below stays testable without it. Asking after the + `--component` rewrite is what keeps the answer one organism's — the same question put to a + Chimera concatenates every Component's ids, and a ribosomal share pooling two organisms is a + share of a population nobody asked about, while this verb is already invoked once per Component. + An annotation no curated list ships for, or one whose list declares no ribosomes, is not a + refusal: it costs the plate that one `obs` column and never a zero one, which is what an + unmeasured thing looks like everywhere else on this path. A curated file that is present and + broken is a different thing, and stays the refusal below. + **`--threads` counts cells at once, and defaults to one.** The counting rule asks the scheduler for threads and hands them over here; a hand invocation that says nothing gets the single-core plate it used to get. The cells are independent and the annotation is read once before the @@ -889,9 +901,9 @@ def io_umi_count( typer.echo(json.dumps({"error": str(exc)}), err=True) raise typer.Exit(2) from exc try: - from genome import ( - Genome, # untyped lab package; resolved here, off the strict workflow path - ) + # untyped lab package; resolved here, off the strict workflow path. The two absences are + # imported by name because they have to be CAUGHT by name -- see the lookup below. + from genome import GeneCategoryNotDeclaredError, Genome, NoGeneCategoriesError except ImportError as exc: # pragma: no cover - depends on the host typer.echo(json.dumps({"error": f"liulab-genome is not importable: {exc}"}), err=True) raise typer.Exit(3) from exc @@ -921,7 +933,24 @@ def io_umi_count( # The gate above leaves exactly one of the two forms and the arm above has filled the # annotation in; the checker cannot see a relation between two options, so it is said here. gtf = Path(str(Genome(assembly).annotations.path(cast(str, annotation)))) - written = write_umi_counts(plate, gtf.with_suffix(".db"), out, workers=threads) + # Which of this annotation's genes are ribosomal, asked AFTER the pair above is settled, so + # a Component asks its OWN genome about its OWN annotation. The same question put to a + # Chimera concatenates worm and bacterial ids across its sources, and a share pooling two + # organisms measures neither; the counter already runs once per Component, so asking the + # component directly makes that impossible and needs no handling of sources here at all. + try: + rrna_gene_ids = Genome(assembly).gene_list("rRNA", annotation).gene_ids + except (GeneCategoryNotDeclaredError, NoGeneCategoriesError): + # No curated list ships for this annotation, or one does and does not declare + # ribosomes: either way the metric is unmeasurable here, which costs the plate that + # column and never a zero one. Caught by NAME rather than as `LookupError`, which an + # unregistered annotation also raises -- swallowing that would turn the refusal below + # into a silently missing column, and a defect in a shipped curated file lands there + # too, as the `ValueError` it is. + rrna_gene_ids = None + written = write_umi_counts( + plate, gtf.with_suffix(".db"), out, workers=threads, rrna_gene_ids=rrna_gene_ids + ) except UmiCountError as exc: typer.echo(json.dumps({"error": str(exc)}), err=True) raise typer.Exit(3) from exc diff --git a/tests/test_cli.py b/tests/test_cli.py index 4f6e044..e9ed5d3 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -5,10 +5,16 @@ import json import random from pathlib import Path +from types import SimpleNamespace from typing import Any, cast import pytest import yaml + +# The two ways liulab-genome says an annotation cannot answer which of its genes are ribosomal. At +# module scope because they are parametrize rows below, and constructed rather than named because +# what `io umi-count` must catch is the class — a stand-in raised in their place proves nothing. +from genome import GeneCategoryNotDeclaredError, NoGeneCategoriesError from typer.testing import CliRunner, Result from conftest import ( @@ -123,8 +129,30 @@ def test_the_cli_surface_exits_and_answers_as_documented( assert needle in result.stdout +@pytest.mark.parametrize( + "answer, ribosomal", + [ + pytest.param(["GENE_A"], (1, 1.0), id="the-genomes-ribosomal-ids-reach-the-counter"), + # Two different absences, and the verb owes both the same answer. Neither is a defect: no + # curated list ships for most annotations, and one that ships may simply not name ribosomes. + # The real exceptions rather than stand-ins for them, so what is caught is the class. + pytest.param( + NoGeneCategoriesError("synthetic", "mm10", ()), + None, + id="an-annotation-nothing-is-curated-for-costs-the-column", + ), + pytest.param( + GeneCategoryNotDeclaredError("synthetic", "mm10", "rRNA", ("protein_coding",)), + None, + id="an-annotation-declaring-no-ribosomes-costs-the-column", + ), + ], +) def test_io_umi_count_finds_the_annotation_database_beside_the_gtf_liulab_genome_registered( - tmp_path: Path, monkeypatch: pytest.MonkeyPatch + answer: list[str] | Exception, + ribosomal: tuple[int, float] | None, + tmp_path: Path, + monkeypatch: pytest.MonkeyPatch, ) -> None: """The verb's whole job: marshal arguments, resolve the annotation, and answer on stdout. @@ -134,11 +162,18 @@ def test_io_umi_count_finds_the_annotation_database_beside_the_gtf_liulab_genome a real assembly needs a genome store this box may not have; what is under test is the derivation and the wiring, not liulab-genome. - **`--threads` is the other thing it marshals, and it is the whole of what #397 changed here.** - The counting rule asks the scheduler for threads and renders them into this command; a verb that - accepted the option and dropped it would leave the plate on one core at exit 0, which is the - state this ticket found. So the number is caught on its way into the counter rather than - inferred from how fast one cell counted. + **`--threads` is the other thing it marshals.** The counting rule asks the scheduler for threads + and renders them into this command; a verb that accepted the option and dropped it would leave + the plate on one core at exit 0, which is exactly the state this claim was added to catch. So the + number is caught on its way into the counter rather than inferred from how fast one cell counted. + + **The ribosomal gene ids are the third**, and they are the parameter above. Handed over, they + put both `obs` columns on the written object — the whole point of asking, proved end to end + rather than at the seam. Where the annotation cannot answer, the verb hands over nothing and the + object carries NEITHER column: the rows are the two ways liulab-genome says so, and a `0.0` + share would be a measurement claim about a plate nobody measured. An annotation that is not + registered at all is deliberately not a row here — it raises through the same base class as + these two and must stay the refusal it is, which is why the verb catches these by name. """ import anndata as ad import genome as liulab_genome @@ -188,16 +223,24 @@ def __init__(self, assembly: str) -> None: self.assembly = assembly self.annotations = _StubRegistry() + def gene_list(self, category: str, annotation: str | None = None) -> Any: + assert category == "rRNA", "the ribosomal category is what this metric is of" + if isinstance(answer, Exception): + raise answer + return SimpleNamespace(gene_ids=list(answer)) + monkeypatch.setattr(liulab_genome, "Genome", _StubGenome) - asked_for: list[int] = [] + marshalled: list[tuple[int, Any]] = [] write_counts = counter.write_umi_counts - def note_the_width(cells: Any, db: Path, out: Path, workers: int = 1) -> Path: - asked_for.append(workers) - return write_counts(cells, db, out, workers) + def note_what_was_handed_over( + cells: Any, db: Path, out: Path, workers: int = 1, rrna_gene_ids: Any = None + ) -> Path: + marshalled.append((workers, rrna_gene_ids)) + return write_counts(cells, db, out, workers, rrna_gene_ids) - monkeypatch.setattr(counter, "write_umi_counts", note_the_width) + monkeypatch.setattr(counter, "write_umi_counts", note_what_was_handed_over) written = tmp_path / "plate.h5ad" result = runner.invoke( @@ -207,7 +250,10 @@ def note_the_width(cells: Any, db: Path, out: Path, workers: int = 1) -> Path: ) # fmt: skip assert result.exit_code == 0, result.stdout - assert asked_for == [3], "the verb took a thread count and counted the plate on one core" + assert marshalled == [(3, ["GENE_A"] if ribosomal else None)], ( + "the verb marshals the thread count and the ids the genome answered with, and hands over " + "no ids at all where it could not answer" + ) assert json.loads(result.stdout)["written"] == str(written) adata = ad.read_h5ad(written) assert list(adata.obs_names) == ["cell_a"] @@ -215,6 +261,15 @@ def note_the_width(cells: Any, db: Path, out: Path, workers: int = 1) -> Path: # it is the sparse matrix that was written, and only the cast says so to the checker. counts = cast("Any", adata.X) assert int(counts[0, adata.var_names.get_loc("GENE_A")]) == 1 + if ribosomal is None: + assert not {"n_umis_rrna", "rrna_fraction"} & set(adata.obs), ( + "an annotation that cannot say which genes are ribosomal costs the column, and a plate " + "reading 0.0% ribosomal is the one answer nobody may act on" + ) + else: + molecules, share = ribosomal + assert int(adata.obs["n_umis_rrna"].iloc[0]) == molecules + assert float(adata.obs["rrna_fraction"].iloc[0]) == share def test_io_umi_count_reads_a_components_annotation_off_the_chimeras_completion_record( @@ -233,6 +288,12 @@ def test_io_umi_count_reads_a_components_annotation_off_the_chimeras_completion_ Component that contributed nothing is named rather than counted against nothing — `tinyEcDub` ships no GTF, and neither would a spike-in or a plasmid. + **Which genes are ribosomal is asked of that same rewritten pair**, and that is the second claim + here. Asked one line earlier it would go to the Chimera, whose answer concatenates every + Component's ids — a share of worm molecules over a set holding bacterial genes too, reported as + this Component's. Nothing about the number would look wrong, which is why the recipient of the + question is pinned rather than left to the counter to notice. + `Genome` is stubbed because a real Chimera needs a built genome store this box may not have, and the counter is stubbed out too: what a `.db` beside its `.gtf` actually counts is the neighbouring test's claim, proved there against a synthetic annotation. This one is about which. @@ -244,6 +305,7 @@ def test_io_umi_count_reads_a_components_annotation_off_the_chimeras_completion_ chimera = "tinyCe_tinyEcDub" record: dict[str, str | None] = {"tinyCe": "wormbase_ws298", "tinyEcDub": None} resolved: list[tuple[str, str]] = [] + asked_which_genes: list[tuple[str, str, str | None]] = [] class _StubRegistry: def __init__(self, assembly: str) -> None: @@ -268,8 +330,14 @@ def default_gtf(self) -> str: "a Component's default annotation now is not necessarily what went into the merge" ) + def gene_list(self, category: str, annotation: str | None = None) -> Any: + asked_which_genes.append((self.assembly, category, annotation)) + return SimpleNamespace(gene_ids=["WBGene00000001"]) + monkeypatch.setattr(liulab_genome, "Genome", _StubGenome) - monkeypatch.setattr(counter, "write_umi_counts", lambda cells, db, out, workers=1: out) + monkeypatch.setattr( + counter, "write_umi_counts", lambda cells, db, out, workers=1, rrna_gene_ids=None: out + ) written = tmp_path / "combined.tinyCe.h5ad" result = runner.invoke( @@ -282,6 +350,10 @@ def default_gtf(self) -> str: assert resolved == [("tinyCe", "wormbase_ws298")], ( "the registered name comes off the Chimera's record and the GTF off the Component itself" ) + assert asked_which_genes == [("tinyCe", "rRNA", "wormbase_ws298")], ( + "the ribosomal ids are asked of the Component's own genome and its own annotation — the " + "Chimera would answer with two organisms' genes, and the share would be of neither" + ) refused = runner.invoke( app, From faa595878fe266e045ead9a3230419817d7b0186 Mon Sep 17 00:00:00 2001 From: hq <29302823+lhqing@users.noreply.github.com> Date: Sat, 22 Aug 2026 13:25:32 -0400 Subject: [PATCH 4/4] the prose enumeration and the genome surface catch up with the ribosomal columns MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Review follow-up on the three commits that added the ribosomal share. All tidy-up: no behaviour moves, no artifact changes shape, and the composed pipeline renders the same bytes. The one real breach was `count.py`'s module docstring. It enumerates the `obs` surface in prose a few lines after warning that a surface spelled again in prose is the copy that goes stale — and two columns had landed since anyone reread it. The sentence now carries the ribosomal share as a clause, conditional on the annotation being able to name those genes, and defers the argument to the constant that holds it. `_GENOME_API` in the repo invariants gains `gene_list`. It is the hand-written set of every liulab-genome attribute we call, asserted against the real package so an upstream rename goes red; the verb had grown a call site the set did not know about, which is exactly the drift it exists to catch. `io umi-count` built `Genome(assembly)` twice in one `try` — once for the annotation path and once for the gene list. Bound once, after the `--component` rewrite has settled the pair, which is the ordering both lookups depend on and a test pins. `N_UMIS_RRNA` leaves the `_obs_columns` whitelist. `fate_metrics` never reads that key — the page carries the share and not the count beside it, which was settled — so listing it read a column back for nobody. The count stays on the written object, and the metric table's claim that it earns no report row is unchanged. Prose: the Chimera-concatenates-every-Component's-ids argument was made in full twice in `io.py`, and "neither column rather than a zero" three times over the one at `RRNA_FRACTION`. Each keeps one home and the rest point at it. Tests: `test_cli.py` names the two columns by the constants the counter exports rather than by string literal, and its happy-path row drops the arithmetic. What those columns hold is proved against a synthetic annotation in the counter's own tests, so the CLI row was a second test going red for one cause; what it claims on its own — that the ids the genome answered with reach `write_umi_counts` and that the columns land on the object — it still asserts. Co-Authored-By: Claude Opus 5 (1M context) --- src/seqforge/cli/io.py | 14 +++++------ src/seqforge/workflows/umite/count.py | 19 +++++++-------- tests/test_cli.py | 34 ++++++++++++++------------- tests/test_repo_invariants.py | 1 + 4 files changed, 35 insertions(+), 33 deletions(-) diff --git a/src/seqforge/cli/io.py b/src/seqforge/cli/io.py index 77a2bfe..cfacd4a 100644 --- a/src/seqforge/cli/io.py +++ b/src/seqforge/cli/io.py @@ -930,16 +930,16 @@ def io_umi_count( # pair becomes the Component's own and the resolution below is the single one both # forms take -- and the refusal below names the Component whose GTF went missing. assembly, annotation = component, registered + # Bound once, AFTER the pair above is settled, so both lookups below are put to the same + # genome -- under `--component` that is the Component's own and never the Chimera's. + genome = Genome(assembly) # The gate above leaves exactly one of the two forms and the arm above has filled the # annotation in; the checker cannot see a relation between two options, so it is said here. - gtf = Path(str(Genome(assembly).annotations.path(cast(str, annotation)))) - # Which of this annotation's genes are ribosomal, asked AFTER the pair above is settled, so - # a Component asks its OWN genome about its OWN annotation. The same question put to a - # Chimera concatenates worm and bacterial ids across its sources, and a share pooling two - # organisms measures neither; the counter already runs once per Component, so asking the - # component directly makes that impossible and needs no handling of sources here at all. + gtf = Path(str(genome.annotations.path(cast(str, annotation)))) + # Which of this annotation's genes are ribosomal, asked of that same settled pair -- what a + # Chimera would answer instead, and why it may not be asked, is in the docstring above. try: - rrna_gene_ids = Genome(assembly).gene_list("rRNA", annotation).gene_ids + rrna_gene_ids = genome.gene_list("rRNA", annotation).gene_ids except (GeneCategoryNotDeclaredError, NoGeneCategoriesError): # No curated list ships for this annotation, or one does and does not declare # ribosomes: either way the metric is unmeasurable here, which costs the plate that diff --git a/src/seqforge/workflows/umite/count.py b/src/seqforge/workflows/umite/count.py index dc52bdd..eccde6f 100644 --- a/src/seqforge/workflows/umite/count.py +++ b/src/seqforge/workflows/umite/count.py @@ -32,8 +32,10 @@ correction in its output shape. As columns they need no leading underscore either: the underscore was there to keep them out of the gene id namespace, and they are not in it any more. Sequencing saturation is a column beside them, what the cell YIELDED — its molecules and how many genes hold -one — is two more, and how many loci each multiply-placed fragment had is an ``obsm`` array: a -per-cell vector rather than a scalar, so it is the one figure ``obs`` cannot hold. +one — is two more, how much of that yield was ribosomal joins them wherever the annotation can say +which genes those are (:data:`RRNA_FRACTION`), and how many loci each multiply-placed fragment had +is an ``obsm`` array: a per-cell vector rather than a scalar, so it is the one figure ``obs`` cannot +hold. **A matrix is materialised when it cannot be derived from the others, and only then.** That is why the combined UMI matrix is here and a combined *read* matrix is not: reads carry no UMI and are @@ -1061,11 +1063,9 @@ def _rrna_umis( property of the assembly's curation and the axis is a property of the GTF that was registered, so the two overlap partially all the time and a partial overlap is still a measurement. - An EMPTY overlap is not, and it is the one case that gets its own answer. Returning zeros there - would put a column on the object saying every cell of this plate is 0.0% ribosomal, which is - indistinguishable from a measurement and is the exact wrong claim — see :data:`RRNA_FRACTION`. - ``None`` reaches :func:`count_plate` as "write neither column", and a reader then finds the metric - absent, which is what an unmeasured thing looks like everywhere else on this path. + An EMPTY overlap is not, and it is the one case that gets its own answer: ``None`` reaches + :func:`count_plate` as "write neither column" rather than as a plate of zeros. The argument for + that is at :data:`RRNA_FRACTION`. The buckets are the ones the cells were already counted into, summed per cell over the matched indices — no second pass over the matrices and no second traversal of a BAM. Genes are walked in @@ -1118,8 +1118,8 @@ def count_plate( already holds for the annotation itself, and what keeps this function strictly typed and testable against a synthetic annotation. Nothing is counted twice for it: the ids select columns of the combined UMI buckets these cells were already counted into. Omit them, or hand over ids none of - which is on the gene axis, and the object carries no ribosomal column at all rather than a zero - one — see :data:`RRNA_FRACTION`. + which is on the gene axis, and the object carries no ribosomal column at all — see + :data:`RRNA_FRACTION`. """ import anndata as ad @@ -1404,7 +1404,6 @@ def _obs_columns(adata: anndata.AnnData) -> dict[str, list[float]]: SATURATION, N_UMIS, GENES_DETECTED, - N_UMIS_RRNA, RRNA_FRACTION, ) if column in frame.columns diff --git a/tests/test_cli.py b/tests/test_cli.py index e9ed5d3..d9002dd 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -28,6 +28,7 @@ ) from seqforge import __version__, kb from seqforge.cli import app +from seqforge.workflows.umite.count import N_UMIS_RRNA, RRNA_FRACTION runner = CliRunner() @@ -132,25 +133,25 @@ def test_the_cli_surface_exits_and_answers_as_documented( @pytest.mark.parametrize( "answer, ribosomal", [ - pytest.param(["GENE_A"], (1, 1.0), id="the-genomes-ribosomal-ids-reach-the-counter"), + pytest.param(["GENE_A"], True, id="the-genomes-ribosomal-ids-reach-the-counter"), # Two different absences, and the verb owes both the same answer. Neither is a defect: no # curated list ships for most annotations, and one that ships may simply not name ribosomes. # The real exceptions rather than stand-ins for them, so what is caught is the class. pytest.param( NoGeneCategoriesError("synthetic", "mm10", ()), - None, + False, id="an-annotation-nothing-is-curated-for-costs-the-column", ), pytest.param( GeneCategoryNotDeclaredError("synthetic", "mm10", "rRNA", ("protein_coding",)), - None, + False, id="an-annotation-declaring-no-ribosomes-costs-the-column", ), ], ) def test_io_umi_count_finds_the_annotation_database_beside_the_gtf_liulab_genome_registered( answer: list[str] | Exception, - ribosomal: tuple[int, float] | None, + ribosomal: bool, tmp_path: Path, monkeypatch: pytest.MonkeyPatch, ) -> None: @@ -168,12 +169,14 @@ def test_io_umi_count_finds_the_annotation_database_beside_the_gtf_liulab_genome number is caught on its way into the counter rather than inferred from how fast one cell counted. **The ribosomal gene ids are the third**, and they are the parameter above. Handed over, they - put both `obs` columns on the written object — the whole point of asking, proved end to end - rather than at the seam. Where the annotation cannot answer, the verb hands over nothing and the - object carries NEITHER column: the rows are the two ways liulab-genome says so, and a `0.0` - share would be a measurement claim about a plate nobody measured. An annotation that is not - registered at all is deliberately not a row here — it raises through the same base class as - these two and must stay the refusal it is, which is why the verb catches these by name. + put both `obs` columns on the written object — that the ids travel the whole way rather than + only as far as the seam. What those columns then HOLD is arithmetic against a known annotation + and is proved once, in the counter's own tests. Where the annotation cannot answer, the verb + hands over nothing and the object carries NEITHER column: the rows are the two ways + liulab-genome says so, and a `0.0` share would be a measurement claim about a plate nobody + measured. An annotation that is not registered at all is deliberately not a row here — it raises + through the same base class as these two and must stay the refusal it is, which is why the verb + catches these by name. """ import anndata as ad import genome as liulab_genome @@ -261,15 +264,14 @@ def note_what_was_handed_over( # it is the sparse matrix that was written, and only the cast says so to the checker. counts = cast("Any", adata.X) assert int(counts[0, adata.var_names.get_loc("GENE_A")]) == 1 - if ribosomal is None: - assert not {"n_umis_rrna", "rrna_fraction"} & set(adata.obs), ( + columns = {N_UMIS_RRNA, RRNA_FRACTION} + if ribosomal: + assert columns <= set(adata.obs), "the ids the genome answered with reached the object" + else: + assert not columns & set(adata.obs), ( "an annotation that cannot say which genes are ribosomal costs the column, and a plate " "reading 0.0% ribosomal is the one answer nobody may act on" ) - else: - molecules, share = ribosomal - assert int(adata.obs["n_umis_rrna"].iloc[0]) == molecules - assert float(adata.obs["rrna_fraction"].iloc[0]) == share def test_io_umi_count_reads_a_components_annotation_off_the_chimeras_completion_record( diff --git a/tests/test_repo_invariants.py b/tests/test_repo_invariants.py index 6dd9246..88b6224 100644 --- a/tests/test_repo_invariants.py +++ b/tests/test_repo_invariants.py @@ -127,6 +127,7 @@ def test_seqforge_defines_no_genome_machinery(src_trees: SrcTrees) -> None: "components", # io split-chimera: which Components a Chimera holds, off the completion record "separator", # io split-chimera: the underscore run those Components' chromosome names carry "component_annotations", # io umi-count --component: what each Component gave the merged GTF + "gene_list", # io umi-count: which of that annotation's genes are the curated `rRNA` category }