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 diff --git a/src/seqforge/cli/io.py b/src/seqforge/cli/io.py index 410f6b7..cfacd4a 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 @@ -918,10 +930,27 @@ 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)))) - written = write_umi_counts(plate, gtf.with_suffix(".db"), out, workers=threads) + 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.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/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..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 @@ -229,6 +231,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 +1051,40 @@ 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: ``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 + 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 — 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,14 @@ 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, + RRNA_FRACTION, + ) if column in frame.columns } @@ -1348,7 +1458,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_cli.py b/tests/test_cli.py index 4f6e044..d9002dd 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 ( @@ -22,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() @@ -123,8 +130,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"], 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", ()), + False, + id="an-annotation-nothing-is-curated-for-costs-the-column", + ), + pytest.param( + GeneCategoryNotDeclaredError("synthetic", "mm10", "rRNA", ("protein_coding",)), + 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( - tmp_path: Path, monkeypatch: pytest.MonkeyPatch + answer: list[str] | Exception, + ribosomal: bool, + tmp_path: Path, + monkeypatch: pytest.MonkeyPatch, ) -> None: """The verb's whole job: marshal arguments, resolve the annotation, and answer on stdout. @@ -134,11 +163,20 @@ 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 — 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 @@ -188,16 +226,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 +253,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 +264,14 @@ 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 + 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" + ) def test_io_umi_count_reads_a_components_annotation_off_the_chimeras_completion_record( @@ -233,6 +290,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 +307,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 +332,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 +352,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, 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 } 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 --------------------------------------------------------------- #