Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
53 changes: 38 additions & 15 deletions docs/studies/confidence_tracking.md
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand All @@ -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.
34 changes: 32 additions & 2 deletions src/raven_toolbox/annotation/delta_g.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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``.
Expand All @@ -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:
Expand All @@ -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

Expand Down Expand Up @@ -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
16 changes: 8 additions & 8 deletions src/raven_toolbox/confidence.py
Original file line number Diff line number Diff line change
@@ -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:
Expand All @@ -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`).
"""
Expand Down
37 changes: 37 additions & 0 deletions tests/test_annotation.py
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down