diff --git a/CHANGELOG.md b/CHANGELOG.md index e6c94f5..51bdc78 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,14 @@ Milestones in the raven-toolbox port. For function-level status see differed only in operand order. This brings `diff_models` in line with MATLAB RAVEN's `diffModels` ([RAVEN #686](https://github.com/SysBioChalmers/RAVEN/pull/686)); a rule cobra cannot parse falls back to the old string comparison, so malformed rules are still compared rather than silently equated. +* **Fix `load_delta_g_csv` recording the ΔG side-car tables' "missing" sentinel as a real measurement.** The + side-car tables encode "no valid ΔG" as `10000000`, and the loader — written for exactly these files, down to its + `Var1`/`Var2` defaults — stamped it verbatim, presenting a physically impossible 10⁷ kJ/mol as a + measurement on **777 of yeast-GEM's 4102 reactions (19.5%)**. yeast-GEM's own `checkrxnDirection.m` gates + on the same value (`if ~isequal(seed_rxnInfo{...},'10000000.0') %check if database contains valid deltaG + value`). The sentinel is now treated as missing, as NaN already was, recognised whichever dtype the CSV + round-trip produces; the new keyword-only `missing_value` (default `DELTA_G_MISSING`) tunes or + disables it. Real ΔG coverage of yeast-GEM is 78.2%, not the 97.1% the loader previously implied. ## 0.3.0 — 2026-07-16 diff --git a/docs/studies/confidence_tracking.md b/docs/studies/confidence_tracking.md index 3570d52..8d69da4 100644 --- a/docs/studies/confidence_tracking.md +++ b/docs/studies/confidence_tracking.md @@ -1,9 +1,9 @@ # Per-reaction confidence tracking -**Status: P1 + P2 shipped.** The data model, the notes round-trip (YAML + SBML), and the `localization`, -`equation` and `gene_association` scorers live in +**Status: the three facets are shipped.** The data model, the notes round-trip (YAML + SBML), and the +`localization`, `equation` and `gene_association` scorers live in [`raven_toolbox/confidence.py`](../../src/raven_toolbox/confidence.py) (tests in `tests/test_confidence.py`). -The `reversibility` facet (P3) and the standards mapping (P4) below remain planned. +The facet set is **closed** — see [§10](#10-what-is-left) for the remaining work. Every reaction carries a small structured record scoring how well-supported each of its **facets** is, persisted in the model file, computed from evidence, updated by curation, and consumed by the raven-toolbox @@ -43,7 +43,6 @@ no facet is left. Those two facts, taken together, force the design: | `localization` | the compartment assignment | DeepLoc support at the assigned compartment + FBA certification | | `equation` | mass & charge balance, formula completeness | `get_elemental_balance` + a recomputed charge sum | | `gene_association` | is there gene evidence, and is it corroborated | GPR presence + a `pubmed` annotation | -| `reversibility` *(planned)* | are the bounds thermodynamically justified | ΔG hook / FVA-attainable direction | Each facet is scored independently, so a model can be annotated one facet at a time and the record grows incrementally. A `ConfidenceEntry` is a continuous `score` in [0, 1] plus optional provenance: a categorical @@ -82,8 +81,8 @@ and must balance. Detecting biomass by *name* is likewise refused — `\bgrowth\ model carries **no** reaction SBO terms at all, the scorers warn: they then cannot tell a pseudo-reaction from a defect. -Note also that `reaction.boundary` is `len(metabolites) == 1`, independent of id, bounds and reversibility — -so it catches an exchange reaction that is not named `EX_`, and a blocked `(0, 0)` one. +Note also that `reaction.boundary` is `len(metabolites) == 1`, independent of id and of bounds — so it +catches an exchange reaction that is not named `EX_`, and a blocked `(0, 0)` one. ### Bands @@ -207,9 +206,9 @@ not a guarantee about another. `confidence_report`, `facet_summary`. Storage lives in `reaction.notes["raven_confidence"]`; there is no separate save step — the record serialises with the model. -**Planned (P3):** `score_reversibility_confidence` and an umbrella `annotate_confidence(model, types=[...])`. +**Planned:** an umbrella `annotate_confidence(model, facets=[...])` — see [§10](#10-what-is-left). -## 9. Standards alignment (P4, for the paper) +## 9. Standards alignment (for the paper) Map the categorical `level` onto the established **Thiele & Palsson reconstruction confidence score (0–4)** so it is familiar to modellers and reviewers, and reference **ECO** (Evidence & Conclusion Ontology) terms where a @@ -223,11 +222,35 @@ facet maps to an evidence class. Two cautions carried forward from the design wo biomass production, `SBO:0000395` encapsulating process, `SBO:0000630` ATP maintenance, `SBO:0000672` spontaneous reaction, `SBO:0000655` transport reaction. -## 10. Phasing +## 10. What is left -- **P1 — foundation + the facet we already had** *(shipped)*: data model, storage/round-trip helpers (YAML - **and** SBML), the `localization` scorer, `confidence_report`. -- **P2 — cheap structural facets** *(shipped)*: `equation` (mass/charge/formula) + `gene_association`, the - abstain-vs-zero discipline, `facet_summary`, and the exemption predicates. -- **P3 — reversibility:** the bounds-vs-FVA heuristic, with a ΔG hook. -- **P4 — paper:** Thiele-Palsson / ECO / SBO mapping and documentation. +The **facet set is closed**: `localization`, `equation` and `gene_association` are shipped, and no further +facet is planned. What remains is finishing the work *around* them. + +### 10.1 Wire the facets together + +- **`annotate_confidence(model, facets=[...])`** — one call that runs every applicable scorer, instead of + making a caller know the three scorer names and their argument shapes. `localization` needs a proposal and + a score table while the other two need only the model, so the umbrella must skip a facet whose inputs are + absent rather than fail — the same abstain-rather-than-guess rule the scores themselves follow. +- **Let `curation_priority` read the record.** Today it re-derives localisation evidence from scratch and + cannot see that a curator already settled a placement, so a `mark_curated` reaction keeps surfacing in the + review queue. Skipping facets at `level == "curated"` closes the loop between the two tools: score → review + → curate → *stop being asked about it*. This is the single change that makes curation feel finished. +- **Point the SBO precondition at its remedy.** The scorers warn on a model with no SBO terms, but do not say + that `raven_toolbox.annotation.add_sbo_terms(model)` is the fix. The warning should name it. + +### 10.2 Validate beyond one model + +Every number in §6 comes from yeast-GEM, and §6 already flags that three bands (`formula-unparseable`, +`formula-generic`, `charge-unknown`) never fire there — they are covered only by synthetic fixtures. A +distribution measured on one model is not a guarantee about another. Running +`scripts/measure_confidence_facets.py` over Human-GEM and a non-curated draft would show whether the bands +are calibrated or merely yeast-shaped, and would exercise the branches yeast cannot reach. The +gene-rubric-vs-`Confidence Level` check only replicates on a model that records that note, so its absence +elsewhere is itself worth reporting. + +### 10.3 Standards alignment for the paper + +§9 above: map `level` onto Thiele & Palsson 0–4, and attach ECO terms where a facet maps to an evidence +class — with the two cautions recorded there. The SBO half is already done and verified. diff --git a/src/raven_toolbox/annotation/delta_g.py b/src/raven_toolbox/annotation/delta_g.py index e0a90e0..aef343d 100644 --- a/src/raven_toolbox/annotation/delta_g.py +++ b/src/raven_toolbox/annotation/delta_g.py @@ -18,6 +18,25 @@ import cobra import pandas as pd +#: The "no valid ΔG" sentinel used by the ΔG side-car tables (e.g. yeast-GEM's +#: ``model_rxnDeltaG.csv``). yeast-GEM's ``checkrxnDirection.m`` gates on it verbatim: +#: ``if ~isequal(seed_rxnInfo{...},'10000000.0') %check if database contains valid deltaG +#: value``. Stamping it would present a physically impossible 10⁷ kJ/mol as a measurement, +#: so it is treated as missing. +DELTA_G_MISSING = 1e7 + + +def _is_missing(value, sentinel: float) -> bool: + """True when ``value`` is the sentinel, whether it arrived as a number or as text. + + The CSV round-trips through MATLAB and pandas, so the same sentinel shows up as ``10000000``, + ``10000000.0`` or ``"10000000.0"`` depending on the writer and the column's inferred dtype. + """ + try: + return math.isclose(float(value), sentinel, rel_tol=1e-9) + except (TypeError, ValueError): + return False + def load_delta_g_csv( entities: Iterable, @@ -26,6 +45,7 @@ def load_delta_g_csv( id_column: str = "Var1", value_column: str = "Var2", note_key: str = "deltaG", + missing_value: float | None = DELTA_G_MISSING, verbose: bool = False, ) -> int: """Stamp ``note_key`` on each entity from a CSV of ``id → value``. @@ -42,12 +62,19 @@ def load_delta_g_csv( note_key Key under which the value is stored on ``entity.notes``. Default ``"deltaG"``. + missing_value + A sentinel standing for "no value", left unstamped rather than + recorded as a measurement. Defaults to + :data:`DELTA_G_MISSING` (10⁷), which covers 777 of + yeast-GEM's 4102 reaction rows. Pass ``None`` to stamp every + value verbatim. verbose Print a summary of unmatched entity ids. Returns ------- - The number of entities that were stamped (i.e. matched the CSV). + The number of entities that were stamped (i.e. matched the CSV and + carried a real value). """ df = pd.read_csv(path) if id_column not in df.columns or value_column not in df.columns: @@ -64,6 +91,9 @@ def load_delta_g_csv( if value is None or (isinstance(value, float) and math.isnan(value)): missing.append(entity.id) continue + if missing_value is not None and _is_missing(value, missing_value): + missing.append(entity.id) + continue entity.notes[note_key] = str(value) stamped += 1 @@ -111,5 +141,5 @@ def save_delta_g_csv( # Re-export the cobra Model type for type-checker friendliness; helps # IDEs surface the right hints to callers that hand us model.metabolites # / model.reactions directly. -__all__ = ["load_delta_g_csv", "save_delta_g_csv"] +__all__ = ["DELTA_G_MISSING", "load_delta_g_csv", "save_delta_g_csv"] _ = cobra # silence "imported but unused" — used for typing context above diff --git a/src/raven_toolbox/confidence.py b/src/raven_toolbox/confidence.py index f544c0b..f6bdf82 100644 --- a/src/raven_toolbox/confidence.py +++ b/src/raven_toolbox/confidence.py @@ -1,10 +1,10 @@ """Per-reaction, multi-facet confidence — persisted in the model, ignored by plain cobra. -Attaches a small structured record to a reaction scoring how well-supported each *facet* of it is -(``localization``, ``equation``, ``gene_association``; ``reversibility`` follows the same shape). Each -facet is a :class:`ConfidenceEntry` — a continuous 0-1 ``score`` plus optional provenance (a categorical -``level``, the ``basis`` evidence, ``method``/``source``/``note``). A reaction carries a -:class:`ReactionConfidence` (facet → entry) whose ``overall`` is the weakest facet. +Attaches a small structured record to a reaction scoring how well-supported each *facet* of it is: +``localization``, ``equation`` and ``gene_association``. Each facet is a :class:`ConfidenceEntry` — a +continuous 0-1 ``score`` plus optional provenance (a categorical ``level``, the ``basis`` evidence, +``method``/``source``/``note``). A reaction carries a :class:`ReactionConfidence` (facet → entry) whose +``overall`` is the weakest facet. **Two rules govern every score**, because ``overall = min(facets)`` and :func:`_write` drops the record when no facet remains: @@ -30,9 +30,9 @@ defect. Detecting biomass by name instead is deliberately *not* done: ``\\bgrowth\\b`` matches "non-growth associated maintenance reaction", and a name regex must never silence a chemistry check. -The design and roadmap (the ``reversibility`` facet, ECO/SBO and Thiele-Palsson mapping) are in -``docs/studies/confidence_tracking.md``. Wire it in by calling :func:`score_localization_confidence` on -an :class:`~raven_toolbox.localization.AssignmentProposal`, :func:`score_equation_confidence` and +The design and the measured yeast-GEM distributions are in ``docs/studies/confidence_tracking.md``; the +facet set above is closed. Wire it in by calling :func:`score_localization_confidence` on an +:class:`~raven_toolbox.localization.AssignmentProposal`, :func:`score_equation_confidence` and :func:`score_gene_association_confidence` on any model, and :func:`mark_curated` when a curator firmly fixes a facet (e.g. after :func:`~raven_toolbox.localization.relocate_reactions`). """ diff --git a/tests/test_annotation.py b/tests/test_annotation.py index d61a651..98378d9 100644 --- a/tests/test_annotation.py +++ b/tests/test_annotation.py @@ -185,6 +185,43 @@ def test_load_skips_nan_rows(tmp_path): assert m.metabolites.get_by_id("glc_e").notes["deltaG"] == "1.0" +def test_load_skips_the_missing_sentinel(tmp_path): + """The ΔG side-car tables write 10000000 for "no valid ΔG", and yeast-GEM's side-car carries it on + 777 of its 4102 reaction rows. Stamping it would record a physically impossible 10^7 kJ/mol as a + measurement. yeast-GEM's own checkrxnDirection.m gates on the same value.""" + m = _toy_model() + m.metabolites.get_by_id("atp_c").notes["deltaG"] = "preserved" + + csv = tmp_path / "met_dg.csv" + pd.DataFrame({"Var1": ["atp_c", "glc_e"], "Var2": [10000000.0, 1.0]}).to_csv(csv, index=False) + + stamped = load_delta_g_csv(m.metabolites, csv) + assert stamped == 1 # only the real value + assert m.metabolites.get_by_id("atp_c").notes["deltaG"] == "preserved" + assert m.metabolites.get_by_id("glc_e").notes["deltaG"] == "1.0" + + +def test_sentinel_skipping_can_be_opted_out_of(tmp_path): + m = _toy_model() + csv = tmp_path / "met_dg.csv" + pd.DataFrame({"Var1": ["atp_c"], "Var2": [10000000.0]}).to_csv(csv, index=False) + + assert load_delta_g_csv(m.metabolites, csv, missing_value=None) == 1 + assert m.metabolites.get_by_id("atp_c").notes["deltaG"] == "10000000.0" + + +def test_sentinel_recognised_whatever_dtype_the_csv_round_trip_produces(tmp_path): + """The same sentinel arrives as 10000000, 10000000.0 or "10000000.0" depending on the writer and + the column's inferred dtype -- a string column appears as soon as one row holds text.""" + m = _toy_model() + csv = tmp_path / "met_dg.csv" + pd.DataFrame({"Var1": ["atp_c", "glc_e"], "Var2": ["10000000.0", "-2.5"]}).to_csv(csv, index=False) + + assert load_delta_g_csv(m.metabolites, csv) == 1 + assert "deltaG" not in m.metabolites.get_by_id("atp_c").notes + assert m.metabolites.get_by_id("glc_e").notes["deltaG"] == "-2.5" + + def test_custom_columns_and_note_key(tmp_path): m = _toy_model() m.metabolites.get_by_id("atp_c").notes["dG_kJ"] = "-30.5"