Skip to content

Additive per-source thermodynamic descriptions (reaction free energies + compound pKas) - #263

Merged
samseaver merged 1 commit into
ModelSEED:devfrom
freiburgermsu:additive-thermodynamics-and-pka-descriptions
Jun 4, 2026
Merged

Additive per-source thermodynamic descriptions (reaction free energies + compound pKas)#263
samseaver merged 1 commit into
ModelSEED:devfrom
freiburgermsu:additive-thermodynamics-and-pka-descriptions

Conversation

@freiburgermsu

Copy link
Copy Markdown
Member

Summary

Records every thermodynamic / pKa source side-by-side instead of collapsing them into a single value. Purely additive — no canonical free-energy or served-pKa value changes, and no TSV changes (thermodynamics and pkas are JSON-only fields).

Reaction free energies — thermodynamics dict

Each method now carries an [energy, error, operator] triple. The operator (>, <, =, ?) is that estimate's own thermodynamic direction, computed with the same heuristic as the canonical reversibility applied to that method's dG:

"thermodynamics": {
    "Group contribution": [4.15, 1.22, "="],
    "eQuilibrator":       [-3.46, 0.05, ">"],
    "dGPredictor":        [-3.82, 0.02, ">"]
}
  • dGPredictor (Wang et al. 2021; predictions already staged under Biochemistry/Thermodynamics/dGPredictor/, kJ→kcal /4.184) is recorded additively for every reaction it predicts (27,715), next to GC/eQ. It does not alter the canonical deltag/deltagerr/reversibility.
  • Estimate_Reaction_Reversibility.py is refactored to expose a reusable reversibility_from_energy(). The dev coercion fix (float(rxn_dg)) is preserved, and the refactor reproduces dev's reversibility output exactly on all 56,012 reactions (validated). New Update_Reaction_dGPredictor_Energies.py (additive-only) and Add_Reaction_Thermodynamics_Operators.py; the GC/eQ writers attach the operator.

Compound pKas — pkas dict

New per-tool pkas dict keyed by predicting tool, each {pKa, pKb} in the existing <microstate>:<atom>:<value> encoding, next to the served pka/pkb:

"pkas": {
    "Marvin":  {"pKa": "1:14:12.60;1:22:3.33;...", "pKb": "1:9:-3.03;1:14:-3.85;..."},
    "MolGpKa": {"pKa": "1:14:10.12;1:29:0.98;...", "pKb": "1:6:3.66;1:9:3.42;..."}
}
  • Marvin = the cascade-winning per-database Marvin value; MolGpKa = the OPAM2/MolGpKa prediction (new source file Biochemistry/Structures/ModelSEED/pkas/opam2_molgpka.tsv, 22,399 compounds).
  • The served pka/pkb are unchanged (the existing cascade-picked Marvin value).
  • New Scripts/Structures/Add_Compound_pKa_Sources.py builds the dict (idempotent, never touches served pka/pkb); added to the manifest cascade as the AddPkaSources stage.

Data changed

  • 27,715 reactions gained a dGPredictor record; 75,497 GC/eQ records gained an operator.
  • 35,632 compounds gained a Marvin record; 22,395 a MolGpKa record.
  • Verified: zero changes to canonical deltag/deltagerr/reversibility/notes, served pka/pkb, any other field, or any .tsv.

Docs updated in Scripts/Thermodynamics/README.md and Biochemistry/Structures/sources.yaml.

🤖 Generated with Claude Code

…und pKas)

Record every thermodynamic/pKa source side-by-side rather than collapsing them
into a single value. Purely additive: no canonical free-energy or served-pKa
value changes, and no TSV changes (both `thermodynamics` and `pkas` are
JSON-only fields).

Reaction free energies (`thermodynamics` dict):
- Each method now carries an [energy, error, operator] triple; the operator
  ('>', '<', '=', '?') is that estimate's own thermodynamic direction, computed
  with the same heuristic as the canonical reversibility applied to that
  method's dG, e.g. rxn00001:
    "Group contribution": [4.15, 1.22, "="]
    "eQuilibrator":       [-3.46, 0.05, ">"]
    "dGPredictor":        [-3.82, 0.02, ">"]
- dGPredictor (Wang et al. 2021; staged predictions already in
  Biochemistry/Thermodynamics/dGPredictor/, kJ->kcal /4.184) is recorded
  additively for every reaction it predicts (27,715), next to GC/eQ. It does
  NOT alter the canonical deltag/deltagerr/reversibility.
- Estimate_Reaction_Reversibility.py is refactored to expose a reusable
  reversibility_from_energy(); the dev coercion fix (float(rxn_dg)) is preserved
  and the refactor reproduces dev's reversibility output exactly on all 56,012
  reactions. New Update_Reaction_dGPredictor_Energies.py (additive-only) and
  Add_Reaction_Thermodynamics_Operators.py; GC/eQ writers attach the operator.

Compound pKas (`pkas` dict):
- New per-tool `pkas` dict keyed by predicting tool, each {pKa, pKb} in the
  existing <microstate>:<atom>:<value> encoding, next to the served pka/pkb:
    "pkas": {
        "Marvin":  {"pKa": "...", "pKb": "..."},
        "MolGpKa": {"pKa": "...", "pKb": "..."}
    }
  Marvin = the cascade-winning per-database Marvin value; MolGpKa = the
  OPAM2/MolGpKa prediction (new source file
  Biochemistry/Structures/ModelSEED/pkas/opam2_molgpka.tsv, 22,399 compounds).
  The served pka/pkb are unchanged (the dev cascade-picked Marvin value).
- New Scripts/Structures/Add_Compound_pKa_Sources.py builds the dict (idempotent,
  never touches served pka/pkb); added to the manifest cascade as AddPkaSources.

Data: 27,715 reactions gained a dGPredictor record and 75,497 GC/eQ records
gained an operator; 35,632 compounds gained a Marvin record and 22,395 a MolGpKa
record. Verified zero changes to canonical deltag/deltagerr/reversibility/notes,
served pka/pkb, any other field, or any TSV. Docs: Scripts/Thermodynamics/
README.md and Biochemistry/Structures/sources.yaml.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
@samseaver
samseaver merged commit 8282481 into ModelSEED:dev Jun 4, 2026
Cooper-Taylor added a commit to Cooper-Taylor/ModelSEEDDatabase that referenced this pull request Jun 4, 2026
… + adapted refactor)

Brings PR ModelSEED#263 (origin/dev c887ede, "Add additive per-source thermodynamic
descriptions") onto claude-changes:

DATA (force-copied from origin/dev, unmodified):
- 61 Biochemistry/reaction_*.json shards (each thermodynamics sublist now
  carries [energy, error, operator] triples; new dGPredictor source on
  27,715 reactions)
- 50 Biochemistry/compound_*.json shards (new per-tool `pkas` dict beside
  served pka/pkb; 35,632 Marvin + 22,395 MolGpKa records)
- New source file Biochemistry/Structures/ModelSEED/pkas/opam2_molgpka.tsv
- Updated Biochemistry/Structures/sources.yaml

NEW + NON-CONFLICT SCRIPTS (taken from origin/dev as-is):
- Scripts/Thermodynamics/Add_Reaction_Thermodynamics_Operators.py (NEW)
- Scripts/Thermodynamics/Update_Reaction_dGPredictor_Energies.py (NEW)
- Scripts/Structures/Add_Compound_pKa_Sources.py (NEW)
- Scripts/Structures/Update_Compound_pKas.py (updated to write pkas dict)
- Scripts/Thermodynamics/Rerun_Thermodynamics.sh (chains new steps)
- Scripts/Thermodynamics/README.md
- Scripts/Updates/cascade.py

ADAPTED claude-changes REFACTOR (kept structured around _thermo_helpers.py):
- _thermo_helpers.py: added fmt_dg_dge_op + DEFAULT_OP placeholder. The
  reaction-energy write paths (sum_reaction_energy, the sentinel branch
  in run_reaction_aggregation_update, run_reaction_lookup_update) now
  emit 3-element [dg, dge, '?'] triples; the '?' is later replaced by
  Estimate_Reaction_Reversibility (cascade winner) and
  Add_Reaction_Thermodynamics_Operators (comprehensive backfill).
  Compound-side writes stay 2-element (operators are reaction-only).
- Estimate_Reaction_Reversibility.py: added 'DGP' to DB_LEVEL_LABEL /
  DB_LEVEL_PRIORITY / _parse_db_level. _is_source_eligible falls
  through to the sublist-only check when no legacy note is mapped
  (dGPredictor has no legacy note). Extracted the cascade body into a
  _cascade helper so the new module-level reversibility_from_energy()
  shim — used by the per-source updaters and the operator backfill —
  can run the same heuristic on per-source (rather than top-level)
  energy values.

The two thin per-source updaters (Update_Reaction_GroupContribution_Energies
and Update_Reaction_eQuilibrator_Energies) needed no edits — they
delegate to the helper module, so the helper change is sufficient for
their writes to emit triples.

Verified: claude-changes' Estimate_Reaction_Reversibility (GC, then EQ
pass) against the new-format data reproduces dev's reversibility output
for all 56,012 reactions; 3-element thermodynamics entries with valid
operators are emitted for every populated sublist (55,999 Group
contribution, 19,498 eQuilibrator, 27,715 dGPredictor). The EQ-pass
status/reason column differs cosmetically for 3,792 reactions where
claude-changes' more permissive eligibility logic (sublist OR legacy
note) lets the cascade run instead of falling back to "Incomplete
(GCC)" — the load-bearing operator is identical in every case.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Cooper-Taylor added a commit to Cooper-Taylor/ModelSEEDDatabase that referenced this pull request Jun 4, 2026
…ll 3 sources

PR ModelSEED#263 stores the per-source thermodynamic-direction operator as the
third element of each ``thermodynamics[source]`` triple, computed from
THAT source's own dG via ``reversibility_from_energy``. The first-pass
adaptation (commit 0189e06) wrote a ``'?'`` placeholder and relied on
``Estimate_Reaction_Reversibility._append_direction`` and
``Add_Reaction_Thermodynamics_Operators`` to backfill the real operator
downstream.

That didn't work for two compounding reasons:
  1. ``_append_direction`` writes only the cascade-winning source's
     sublist, using the top-level ``deltag`` value (which is whatever
     ``Update_Reaction_*`` most recently set it to) — for non-winning
     sources the placeholder stayed ``'?'``.
  2. ``Add_Reaction_Thermodynamics_Operators`` is idempotent on
     length-3 entries (it only backfills length-2 pairs). The
     placeholder triples we wrote already had length 3, so the safety
     net silently skipped them — observed live as
     ``Per-method thermodynamics entries seen: 103215; Entries given an
     operator: 0``.

Fix: compute the per-source operator at write time in
``_thermo_helpers.run_reaction_aggregation_update`` and
``_thermo_helpers.run_reaction_lookup_update`` via a lazy import of
``Estimate_Reaction_Reversibility.reversibility_from_energy``. Each
source's stored operator is now a function of THAT source's own dG,
independent of the cascade's choice of top-level ``deltag``. With that
in place, ``_append_direction`` and the dead ``DEFAULT_OP`` /
``fmt_dg_dge_op`` helpers are removed; the upstream backfill script
remains as a no-op safety net for legacy length-2 entries.

Also extends ``Scripts/Tests/test_reaction_direction.py`` so the
regression test exercises the new format end-to-end:
  - Default ``--baseline-ref`` switched from ``dev`` (stale local
    branch) to ``origin/dev`` so the baseline naturally tracks the
    remote tip.
  - ``PIPELINE`` constant grows from 6 to 8 commands to mirror the
    updated ``Rerun_Thermodynamics.sh`` (adds
    ``Update_Reaction_dGPredictor_Energies.py`` and
    ``Add_Reaction_Thermodynamics_Operators.py``).
  - ``compare()`` now reports per-source operator agreement for each
    of Group contribution / eQuilibrator / dGPredictor in addition to
    the cascade-winning top-level ``reversibility``. The PASS-driving
    columns are ``Operator mismatch``, ``Regressed (missing key)``,
    and ``Added in current``. ``Baseline stale (sentinel)`` is shown
    for visibility but is upstream data drift (independently verified
    by applying dev's eligibility logic to the current compound state).

Test result after the fix (refresh-baselined against origin/dev):

  Reactions in baseline : 56012
  Reactions in current  : 56012
  Shared reactions      : 56012
  Top-level reversibility mismatches : 0

  | Source             | Operator match | Op mismatch | Regressed | Added | Baseline stale |
  |--------------------|----------------|-------------|-----------|-------|----------------|
  | Group contribution | 25536/25536    | 0           | 0         | 0     | 276            |
  | eQuilibrator       | 19498/19498    | 0           | 0         | 0     | 0              |
  | dGPredictor        | 27715/27715    | 0           | 0         | 0     | 0              |

  RESULT: PASS

So: zero operator regressions across all three sources and the
cascade-winning reversibility. The 276 ``Baseline stale (sentinel)``
entries on GC are reactions whose committed baseline carries a real GC
sum from an earlier compound-data snapshot; the current compound state
has no GC for cpd00074 / cpd02656 (and similar), so any reagent-sum
aggregator — ours or dev's — produces a sentinel. Verified by walking
the 276 with dev's own ``complete_gc_rxns_dict`` / ``incomplete_gc_rxns_dict``
logic on the current compounds: 276/276 land in ``incomplete`` (= sentinel).

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Cooper-Taylor added a commit to Cooper-Taylor/ModelSEEDDatabase that referenced this pull request Jul 6, 2026
Integrate 36 upstream commits from origin/dev (through 4b1f4c7). Conflict
resolution:

* The 4 conflicting Thermodynamics scripts (Estimate_Reaction_Reversibility,
  Update_Reaction_{GroupContribution,eQuilibrator}_Energies,
  Add_Reaction_Thermodynamics_Operators) are kept from claude-changes: our
  refactor already implements upstream PR ModelSEED#263's additive
  [energy, error, operator] format (via _thermo_helpers), plus our cascade
  bug-fixes and most-specific-first reordering. Taking ours loses no upstream
  script behavior.

* The 106 conflicting reaction data shards are taken from dev to preserve
  upstream curation (formula-conflict resolutions, SMILES canonicalization,
  equation/definition reprints). Our branch's only reaction-shard changes were
  thermo/reversibility fields, which are regenerated in the following commit.

Also brings in upstream's new tooling (Check_Database_Consistency.py,
Recanonicalize_SMILES.py, Report_Formula_Change_Impact.py, Compute_ACP_Overrides.py,
Curate_Structure_Pick.py) and curation data.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants