Skip to content

Document dGPredictor KEGG carry-over defect in the NAR 2026 Methods drafts - #286

Open
freiburgermsu wants to merge 2 commits into
ModelSEED:devfrom
freiburgermsu:methods-dgpredictor-kegg-carryover
Open

Document dGPredictor KEGG carry-over defect in the NAR 2026 Methods drafts#286
freiburgermsu wants to merge 2 commits into
ModelSEED:devfrom
freiburgermsu:methods-dgpredictor-kegg-carryover

Conversation

@freiburgermsu

Copy link
Copy Markdown
Member

Summary

Expands the multi-source thermodynamics Methods draft (Papers/NAR_Update_2026/drafts/methods_multi_source_thermodynamics.md) and its inlined copy in MANUSCRIPT.md with the discovery and correction of a systematic identifier carry-over defect in the original (KEGG-mediated) dGPredictor integration.

The finding being documented

  • Of the 27,715 reactions carrying a stored dGPredictor record, 17,271 (62%) hold an estimate computed for a KEGG reaction id that is not among the reaction's aliases (the reaction records and Unique_ModelSEED_Reaction_Aliases.txt agree on all 17,271 cases).
  • Root cause: a stale-variable bug in the 2023 staging notebook — the KEGG-alias extraction loop never reset its holding variable between reactions, so every KEGG-less reaction silently inherited the id of the nearest preceding KEGG-bearing reaction (e.g. rxn00019, a nitroalkane oxidase, carries RuBisCO's R00024 energy of −8.63 ± 0.04 kcal/mol; rxn13478–rxn13481 all inherited R09245 from rxn13477).
  • A database-wide fingerprint test explains 100.0% of the 17,784 contaminated staged entries by carry-over, with zero exceptions.
  • The contaminated subset shows null correlation with eQuilibrator and Group Contribution (r = −0.001 / −0.006) versus r = 0.61 / 0.74 for the correctly-mapped subset — explaining the method's previously puzzling low agreement with other sources.
  • Because Update_Reaction_dGPredictor_Energies.py never validates staged KEGG ids against aliases, every Rerun_Thermodynamics.sh pass re-propagated the contamination, and the promotion policy's lowest-uncertainty rule favored the deceptively tight copied errors: 7,466 contaminated values were promoted into canonical deltag/deltagerr/reversibility.

What this PR changes

Documentation only — no data or script changes:

  • New Methods subsection (defect finding → forensic confirmation → why it persisted → correction path), styled to mirror the existing GC Convention-A subsection.
  • The retrained dGPredictor-ModelSEED source is noted as structurally immune (keyed directly by ModelSEED rxn id, no KEGG detour) and unaffected.
  • [TBD] markers on the correction pieces not yet landed: regenerated staging JSONs (corrected extraction script lives in the dGPredictor repo), an alias-validation guard in the ingest script, and the re-promotion pass repairing the 7,466 canonical triples. The draft's "Open loose ends" list tracks these, plus a suggested regression test pinning known carry-over cases.
  • Results [NUMBERS PENDING] stub now flags that the head-to-head ΔrG′ comparison must be computed from post-repair records.

🤖 Generated with Claude Code

https://claude.ai/code/session_01XoJyv1ukSyPfME8WaHewJh

freiburgermsu and others added 2 commits August 20, 2026 20:00
Expand the multi-source thermodynamics Methods draft (and its inlined
copy in MANUSCRIPT.md) with the audit finding that 17,271 of 27,715
stored dGPredictor records (62%) carry an estimate for a KEGG reaction
id absent from the reaction's aliases, caused by a stale-variable bug
in the 2023 staging notebook: the alias-extraction loop never reset its
holding variable, so KEGG-less reactions inherited the previous
reaction's id. A database-wide fingerprint test explains 100.0% of
contaminated entries by carry-over (zero exceptions); the contaminated
subset shows null correlation with eQuilibrator/GC (r ~ 0.00 vs
0.61/0.74 for correct mappings), and 7,466 contaminated values were
promoted into canonical deltag/deltagerr/reversibility.

The new subsection documents the defect, forensic confirmation, why
reruns re-propagated it (no alias validation in
Update_Reaction_dGPredictor_Energies.py), and the correction path,
with [TBD] markers on the pieces not yet landed (regenerated staging
JSONs, alias-validation guard, re-promotion pass). The retrained
dGPredictor-ModelSEED source is noted as structurally immune (keyed
directly by rxn id). Results stub now flags that the head-to-head
comparison must use post-repair records.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01XoJyv1ukSyPfME8WaHewJh
samseaver pushed a commit that referenced this pull request Sep 4, 2026
…redictions

Three changes, all to thermodynamics["dGPredictor"].

1. RETRAINED, because the model was broken.

The predecessor dGPredictor-ModelSEED record had a median uncertainty of 88.6
kJ/mol (259.5 for quinones), 44.5% of predictions past 100 kJ/mol, worst error
8,531. The cause was RDKit canonicalization drift: the Wang lab shipped a
fragment vocabulary written by an older RDKit (CC(C)=O) and current RDKit emits
CC(=O)C for the same fragment, so reindexing silently dropped most fragment
counts. BayesianRidge then extrapolated from its prior, and the prior
covariance is what those error bars were.

The fix is a guarantee rather than a patch: the vocabulary is rebuilt from
current-RDKit output, and training and prediction call the same
count_substructures, so they can no longer disagree about a fragment's name.

2. FITTED TO THE SAME MEASUREMENTS AS eQUILIBRATOR.

Training moves from the raw Zenodo TECRDB.tsv (4,544 rows) to TECRDB_dedup.csv
(4,456) -- the table component-contribution's cc_params_dedup.npz was fitted
to. The 88 removed rows are redundant copies of 78 measurement groups counted
2-6 times; one copy of each is kept, so no measurement is lost. Fitting the two
methods to different row sets of the same measurements is not defensible when
the paper compares them head to head.

Training RMSE rises 5.4 -> 5.85 kJ/mol, which is what removing rows the fit got
for free must do. The question needs held-out data, so scripts/08 runs the
comparison eQuilibrator ran for component-contribution: 10 folds, both arms
predicting an identical test set, every copy of a tested measurement withheld
from both arms so the raw arm cannot win on leakage. Over 3,633 held-out
measurements, median |error| 1.7125 -> 1.6952, closer on 53.3%, paired Wilcoxon
p=1.5e-07. Highly significant and small -- the point is that consistency with
eQuilibrator costs nothing measurable, not that de-duplication is a large win.
Per-row errors in crossval_training_sets.tsv.

3. NO COVERAGE FLOOR. Every prediction ships with its own uncertainty.

The earlier policy installed only coverage >= 0.9 (8,279 reactions), on the
reasoning that below that a prediction is prior-dominated and is "we don't
know" wearing a wide error bar. Measured against eQuilibrator on the 20,691
reactions both methods score, that does not hold: the fraction within one
combined sigma is 67-79% in EVERY coverage stratum against 68% nominal. At
coverage < 0.5 the model says +/-111 kJ/mol and is off by 87. The bars are
honest, and mildly conservative if anything.

Coverage also turned out to be the wrong axis to filter on. Spearman against
|deviation from eQuilibrator| is -0.349 for coverage but +0.650 for sigma, and
the two correlate only -0.498 with each other. The 0.9 floor was admitting
1,947 predictions with high sigma (median deviation 43.3 kJ/mol) while
rejecting 5,038 with low sigma (median deviation 12.3).

So all 29,617 reactions and 30,185 compounds are installed, and consumers
filter on sigma. Coverage is stored as provenance -- it explains why sigma is
large -- but the README says plainly that it is not the filter, and gives a
sigma-threshold table instead. DGP_COVERAGE_FLOOR reinstates a floor if one is
ever wanted; it defaults to 0.

This does not smuggle false confidence into the direction field. DGP_HEURISTICS
(the Noor 2012 reversibility index at one sigma, merged in #285) reads the wide
bars and declines: of 29,617 records, 20,723 are "=", 6,870 ">" and 2,024 "<".

Record shape gains coverage, appended so positional readers are unaffected:
  reaction   [dg, err, operator, coverage]
  compound   [dg, err, coverage]         no operator; formation energy has no direction

Both ingest scripts, and the swap, now pass source= to
reversibility_from_energy. Without it a dGPredictor energy is silently scored
with Group-Contribution rules, and the stored operator disagrees with what
Estimate_Reaction_Reversibility.py computes for the same record. This was
latent in the pre-#285 code and would have been wrong the moment #285 landed.

Coverage against the other sources, on 56,012 reactions:
  Group contribution  56,002      eQuilibrator  21,853      dGPredictor  29,617
Compounds gain dGPredictor as a third source: 30,185 records.

dGPredictor and eQuilibrator overlap on 7,247 reactions at the old floor and
20,691 without it -- that overlap is the honest denominator for any head-to-head
comparison. dGPredictor reaches 1,032 reactions eQuilibrator cannot (518
unbalanced, 384 outside the CC span), so neither method dominates.

Also lands:
  dGPredictor-2020.tsv       27,715 archived KEGG-trained predictions, the
                             paper's baseline. Read with care: 17,271 of them
                             (62%) hold an energy computed for a KEGG reaction
                             that is not among that reaction's aliases -- the
                             2023 staging notebook's stale-variable carry-over
                             (PR #286). Confirmed here independently: only
                             10,444 rows carry a matching KEGG id.
  retrained_dG*.json         staged inputs, all predictions
  crossval_training_sets.tsv per-row held-out errors, both training sets
  README.md                  provenance, the calibration evidence, and what to
                             filter on

Removed: Update_Reaction_dGPredictor_ModelSEED_Energies.py, the 61 json_files
and 61 pickle_files staging shards, and modelseed_retrained_dG.json.

07_atomic_swap.py now refuses to overwrite an existing dGPredictor-2020.tsv.
Re-running it after a swap would otherwise capture the RETRAINED values as the
"2020 baseline" and destroy the only copy of the pre-retrain record.

Canonical deltag/deltagerr/reversibility are untouched and are now stale with
respect to every source beneath them. Re-promotion is deliberately out of scope:
the field is being retired.

Co-Authored-By: Claude Opus 5 (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.

1 participant