Skip to content

feat(verab): methoxy-aromatic O-demethylation rule discovery + genome degradation screening (Pickaxe)#48

Open
VibhavSetlur wants to merge 13 commits into
cshenry:mainfrom
VibhavSetlur:feature/verab-methoxy-pickaxe
Open

feat(verab): methoxy-aromatic O-demethylation rule discovery + genome degradation screening (Pickaxe)#48
VibhavSetlur wants to merge 13 commits into
cshenry:mainfrom
VibhavSetlur:feature/verab-methoxy-pickaxe

Conversation

@VibhavSetlur

@VibhavSetlur VibhavSetlur commented Jul 22, 2026

Copy link
Copy Markdown

Summary

Implements the verAB use case from the AI meeting: discover the Pickaxe rule(s) that reproduce verAB aromatic O-demethylation, then screen methoxylated aromatic compounds against genome gene content to predict degradation capability (ADP1 and others). Wraps Pickaxe / MINE-Database (Shebek et al. 2023; mechanism-informed operators, Pate et al. 2026) in KBUtilLib under a single kbu.verab facade + kbu verab CLI, and emits reproducible KING coscientist input files + prompts.

⚠️ Stacks on #47. This branch is built on top of PR #47 (unified thermodynamics/cheminformatics), which is still open. Until #47 merges, this PR's diff also shows #47's commits. Merge #47 first, then this diff reduces to the verAB-only changes (23 files, ~6.6k lines). The verAB-only commits are 3fa9a63 → e235519.

The scientific workflow (Dr. Henry's ask)

Phase 1 — verAB rule discovery

  • Canonical SMILES for the 5 seeds baked in: vanillate, isovanillate, guaiacol, 4-methoxybenzoate, veratrate.
  • kbu verab discover runs Pickaxe expansion and identifies which operator(s) realize Ar-OCH3 → Ar-OH + HCHO, using the operator provenance Pickaxe records per reaction.
  • Rule matching is structural (RDKit reaction-pattern match), with a text-match fallback when RDKit is absent.
  • Emits KING coscientist artifacts (kbu verab emit-king): seeds.tsv/csv, discovered_rules.tsv, target_transformation.txt, prompt.md, manifest.json — ready to submit to KING.

Phase 2 — genome/degradation screening

  • kbu verab enumerate: RDKit substructure scan over the ModelSEED biochem DB to find ALL methoxylated aromatics (works around the exact-match-only compound search).
  • kbu verab screen: for each predicted product answers (a) is the reaction already in the DB? (b) is the product in the DB? (c) does the product have a downstream degradation pathway? (d) is that pathway in our models?
  • predict_genome_degradation: maps the operator's EC (default EC 1.14.13.82, vanillate monooxygenase) to gene content per genome (ADP1 etc.) to predict which genomes can degrade which compounds.

API surface

kbu.verab.seed_compounds()                 # the 5 canonical methoxy-aromatic seeds
kbu.verab.discover_rules(generations=1)     # -> VerabDiscoveryResult (firing operators + confidence + EC hint)
kbu.verab.enumerate_methoxy_aromatics()     # RDKit substructure scan of biochem DB
kbu.verab.screen(...)                        # cross-reference products vs DB / models
kbu.verab.emit_king_workflow(outdir=...)     # KING coscientist input files + prompt
kbu.verab.status()                           # rdkit / minedatabase / backend availability

CLI: kbu verab discover|enumerate|screen|emit-king (all support --json).

Architecture

  • New sub-package src/kbutillib/cheminformatics/verab/ (smarts, models, substructure, rule_discovery, screening, king_artifacts) — extends the existing real PickaxeBackend, never duplicates minedatabase calls.
  • src/kbutillib/verab_utils.py facade (VerabUtils/VerabUtilsImpl) composed from network_expansion + biochem + model + genome + annotation, wired as a lazy kbu.verab toolkit property.
  • Optional KING app bundle at src/kbutillib/king_app_verab/.
  • No top-level rdkit/minedatabase imports anywhere — full graceful degradation via BackendUnavailableError.

Dependencies / Python

  • RDKit is optional (pip install KBUtilLib[cheminformatics] or conda). minedatabase's python_requires<3.10 is a pip-metadata gate only — bypassed by KBUtilLib's DependencyManager sys.path injection (runs on 3.11/3.12). MINE-Database resolves via dependencies.yaml. Documented in docs/thermo-cheminformatics/INSTALL.md.

Tests

  • 158 passed / 13 skipped across the verAB suite (unit tests use fake backends; the 13 skips are the gated live-integration test + RDKit-absent-path tests).
  • tests/test_verab.py, test_verab_screening.py, test_verab_king.py, test_verab_cli.py, test_verab_king_bundle.py, and gated test_verab_live.py (real MINE-Database + RDKit; KBUTILLIB_LIVE_CHEM=1).
  • Zero new failures in the broad suite: verified against a pristine upstream/main worktree — baseline 339 failed / 1349 passed (pre-existing optional-dep gaps: cobra, requests_toolbelt, httpx, etc.); this branch 322 failed / 1572 passed (+223 passing = the new tests, −17 failures, same 59 pre-existing collection errors).

Follow-ups / deferred

  • Live empirical confirmation of the exact firing operator id (test_verab_live.py) requires a MINE-Database checkout + RDKit in CI.
  • The Pate-2026 mechanism-informed operator set isn't packaged in stock minedatabase; supported via a custom rule_list (config cheminformatics.verab.operator_rule_tsv); the operator is confirmed empirically, not hardcoded.
  • Phase-2 substructure scan depends on ModelSEED SMILES coverage; compounds lacking SMILES are skipped with counts.

Update: mechanism-informed operators + robustness hardening

Following a critical self-audit against Dr. Henry's exact asks, five fixes were added so the feature genuinely delivers the requested scientific workflow (not just infrastructure):

FIX 1 (the core ask): Pate-2026 mechanism-informed operators are now the default

  • Henry explicitly asked for the most accurate operators (Pate et al. 2026, JCIM). These are not bundled in stock Pickaxe. They ARE published in the paper's official data repo: stefanpate/coarse-grain-rxns (data/processed/rules/mechinferred_dt_019_rules_w_coreactants.tsv, ~14.5k operators).
  • verAB now defaults to rule_set="mechinformed", resolved at runtime via cheminformatics.verab.operator_rule_tsv config → KBUTILLIB_VERAB_OPERATOR_TSV env → DependencyManager (repo added to dependencies.yaml).
  • The operator TSV is NOT vendored — that repo has no LICENSE, so redistribution isn't permitted. If it isn't resolved, verAB falls back to the bundled metacyc_intermediate set (more mechanism-aware than metacyc_generalized) with an honest warning, and records the rule set actually used in the discovery result + KING manifest.
  • Coreactant roles (WATER, METHYL_DONOR_CoF, …) are resolved from a KBUtilLib-authored role→SMILES table (our own code, not repo data) or a user-supplied operator_coreactant_tsv.
  • The verAB O-demethylation operator is still confirmed empirically via RDKit at runtime (candidate ids 2783_1/7534_1/4304_1 recorded only as hints, never hardcoded as the answer).
  • Install steps documented in docs/thermo-cheminformatics/INSTALL.md.

FIX 2: firing-operator attribution preserved as a list

Pickaxe reactions can fire multiple operators; the ingest previously collapsed them to one ;-joined string. Now PredictedReaction.operators keeps the full list (scalar operator retained for back-compat), so per-operator attribution isn't lost through discovery/screening/KING.

FIX 3: phase-2 screening distinguishes lookup FAILURE from not-found

DB query errors were silently swallowed to {} (false negatives). ScreeningRecord now carries lookup_ok + lookup_errors, so a failed lookup is visible instead of being reported as "not in DB".

FIX 4: genome degradation prediction is multi-EC + explainable

Was single hardcoded EC (1.14.13.82) + a bare bool. Now checks multiple EC terms (derivable from operators), and reports which EC/MSRXN matched which gene/feature per genome, keeping 1.14.13.82 only as a fallback hint.

FIX 5: KING artifacts aligned + honest

Consolidated the installable KING bundle to src/kbutillib/king_app/verab/ (validates against king_install.load_bundle); the manifest records the actual rule set used; prompt.md honestly states the dir is reproducible inputs + a ready-to-run prompt (not itself a KING bundle) and names the mechanism-informed default + fallback.

Tests

255 passed / 20 skipped across the cheminformatics + thermo + verab surface; no regressions. New unit tests cover operator resolution precedence, TSV normalization, fallback + honest warning, operator-list preservation, lookup-failure vs not-found, and multi-EC explainable genome prediction. The gated live test (KBUTILLIB_LIVE_CHEM=1 + RDKit + a MINE-Database/operator checkout) empirically confirms the firing operator.

Honest limitations (not yet delivered)

  • Empirical confirmation of the exact mechanism-informed operator id requires running the live path (RDKit + minedatabase + the cloned operator TSV). Offline/CI runs use the metacyc_intermediate fallback and say so.
  • The coreactant role→SMILES table is best-effort; cofactor-sensitive expansions may benefit from a user-curated operator_coreactant_tsv.
  • Phase-2 substructure enumeration depends on ModelSEED SMILES coverage; compounds without SMILES are skipped (counted).
  • KING compatibility is validated against load_bundle's schema; end-to-end execution inside KING coscientist has not been run here.

Add a PredictiveThermoUtils facade that estimates standard transformed
Gibbs energies of formation (compounds) and reaction (delta_r G'^0) by
dispatching across a priority-ordered set of pluggable backends, each of
which degrades gracefully and never fabricates a value (unknown -> None
plus a warning).

Backends (src/kbutillib/thermo_predictors/):
- equilibrator: live component-contribution via equilibrator-api
  (optional dep). Condition-aware (pH, ionic strength, T, p_Mg).
- dgpredictor: Maranas/Tyo-lab group-contribution reaction predictor
  (Andrew Freiburger's fork). Runs out-of-process against a local clone
  to keep rdkit/openbabel and the ~150 MB model out of this process;
  KEGG stoichiometry in, kJ/mol transformed dG' out. Detects a
  missing/LFS-pointer model and reports unavailable.
- modelseed_db: baked transformed dG from the expanded-eQuilibrator
  ModelSEED Database fork, read directly from compounds.tsv/reactions.tsv
  with the stdlib csv reader (kcal->kJ, fixed pH7/I0.25/298.15 K,
  EQU/EQC/EQP-preferred provenance gate).
- modelseed: always-on fallback over the existing ThermoUtils lookups.
- molgpk: stub for the Tyo-lab molGPK/OPAM2 predictor (unreleased);
  reports unavailable until configured.

Dispatch order: equilibrator -> dgpredictor -> modelseed_db -> modelseed.
Nothing in the subpackage imports a heavy/optional dependency at import
time. Exposed as the lazy `predictive_thermo` property on KBUtilLib.

Adds [project.optional-dependencies] equilibrator / thermo extras and a
dependency-free test suite (tests/test_predictive_thermo.py) covering
graceful degradation, the no-fabrication invariant, dispatch fall-through,
and the kcal->kJ conversion against worked ATP / rxn00001 values.
…ules backends

Add a pluggable in-silico reaction-prediction / metabolic-network-expansion
layer mirroring the predictive-thermo facade pattern:

- cheminformatics/base.py: ExpansionBackend protocol + PredictedCompound/
  PredictedReaction/ExpansionResult dataclasses + BackendUnavailableError.
  Stdlib-only; importing it pulls no heavy dependency.
- cheminformatics/pickaxe_backend.py: optional-dep wrapper over minedatabase
  (Pickaxe) SMARTS-operator expansion; robust rule-data-dir resolution
  (config -> env -> wheel -> source checkout) since the wheel ships no rule
  data; graceful degradation when absent.
- cheminformatics/retrorules_backend.py: applies RetroRules SMARTS via RDKit
  (no-AddHs + partial-sanitize recipe), streaming the flat-rules TSV; filters
  by diameter/direction; graceful degradation when rdkit/TSV absent.
- network_expansion_utils.py: NetworkExpansionUtils facade + Impl wrapper with
  backend dispatch (pickaxe -> retrorules), backend_status diagnostics, and
  never-fabricate semantics (empty result + warnings when nothing expands).
- Wire NetworkExpansionUtils(Impl) into package __init__ and the toolkit
  (lazy 'network_expansion' property); add 'cheminformatics' optional-deps
  extra (rdkit).
- tests/test_network_expansion.py: import/degradation invariants, injected-fake
  dispatch success path, and opt-in live integration tests gated on env vars.

16/16 tests pass (incl. live runs against real RetroRules + MINE-Database
rule data); new modules ruff clean.
Replace the dead-endpoint molGPK stub with a working backend that wraps
Andrew Freiburger's OPAM2 (MolGpKa fine-tuned on the ModelSEED database).

- Subprocess-isolated so KBUtilLib never imports torch/torch-geometric;
  the worker is bundled in-package (_molgpk_worker.py) and run with
  OPAM2_SRC pointing at a cloned OPAM2 repo, so users only clone OPAM2
  and set thermo.molgpk.repo_path/python (or MOLGPK_REPO/MOLGPK_PYTHON).
- compound_dgf -> pKa values + major microspecies; dgf left None (molGPK
  feeds equilibrator, so ΔGf is never double-counted or fabricated).
- Probes worker + ModelSEED-tuned .pth checkpoints with a Git-LFS-pointer
  size guard; degrades gracefully (per-compound miss != backend down).
- Batch path: worker accepts {compounds:[...]} and loads the model once;
  compounds_dgf + facade.compounds_microspecies amortize the model load
  across a whole compound set.

Tested: hermetic unavailable/degrade/LFS-pointer/hard-failure/misalign
cases plus env-gated live single + batch runs against a real OPAM2
checkout (acetate, glycine zwitterion, phosphate).
dGPredictor backend rewritten to Andrew Freiburger's ModelSEED fork API
(ModelSEEDdGPredictor.predict_reaction / predict_from_equation): addresses
reactions by rxnNNNNN accession or ModelSEED cpd stoichiometry, fixes the
model filename to model/modelseed_M12_model_BR.pkl, adds a batch path, and
ships the previously-missing subprocess worker (_dgpredictor_worker.py).

molGPK backend confirmed against the real OPAM2 API (predict_pka.predict /
protonate.Opam_protonate_mol + modelseed .pth weights).

Both validated end-to-end through the PredictiveThermoUtils facade with
config-driven repo paths (thermo.{dgpredictor,molgpk}.{repo_path,python}):
  molgpk  acetic acid pKa 2.01 -> CC(=O)[O-]; glycine -> [NH3+]CC(=O)[O-]
  dgpredictor  rxn00001 -15.76+/-3.63 kJ/mol; ATP hydrolysis -25.46 kJ/mol

Adds gated live integration tests (skip unless *_REPO/*_PYTHON set) and a
documented thermo: config block. 31 unit tests pass; 15 pass with repos
configured.
…-007

Andrew (relaying the Tyo lab) confirmed github.com/tyo-nu/MINE-Database as the
pickaxe version to use. It IS the minedatabase package PickaxeBackend already
wraps, with the bundled rule TSVs — so no separate tyo_pickaxe_backend.py is
needed; point cheminformatics.pickaxe.data_dir at the clone's minedatabase/data.

Validated end-to-end through the NetworkExpansionUtils (kbu.chem) facade:
  D-glucose, 1 generation, metacyc_generalized -> 901 new compounds /
  1696 new reactions. Live integration test test_pickaxe_live_expansion passes.

Resolves ADR-007 (wrap, don't reimplement). Documents the pickaxe/retrorules
config block. Fixes test_backend_status_reports_all_backends to also clear the
DGPREDICTOR_* env vars so it tests the not-configured contract deterministically.
45 unit tests pass; live thermo + pickaxe tests pass with repos configured.
Make the dGPredictor / molGPK / pickaxe backends self-contained the way the rest
of KBUtilLib's external deps are: register the three research repos in
dependencies.yaml (dGPredictor, OPAM2, MINE-Database) with path + git URL, and
have each backend fall back to DependencyManager.get_dependency_path() when no
explicit config key / env var is set. The pickaxe backend also puts a declared
MINE-Database checkout on sys.path itself (no pip install, no .pth shim needed --
the Tyo repo pins python <3.10 which blocks pip install -e on 3.11).

Net result: with zero env vars and zero config paths, all three backends resolve
and compute purely from dependencies.yaml -- the same mechanism used for
ModelSEEDpy / ModelSEEDDatabase / cobrakbase. Adds base.dependency_repo_path()
helper, an INSTALL.md documenting the common footprint, and tightens the
hermetic 'not configured' tests to also neutralize the DependencyManager so they
stay deterministic. 45 unit tests pass; live thermo + pickaxe tests pass.
…eature

Finalization pass to make the equilibrator + dGPredictor + molGPK/OPAM2
unified facade PR-ready.

- docs/config: resolve API-name drift; document real attributes
  (kbu.predictive_thermo, kbu.network_expansion) and clarify legacy kbu.thermo;
  add kbu.chem alias on the toolkit facade
- INSTALL.md: fix broken model-rebuild reference to dGPredictor's in-repo
  retrain_modelseed.py so the documented install path is real
- thermo_predictors: add machine-readable `units` field (kJ/mol) to
  CompoundThermoEstimate/ReactionThermoEstimate and to_dict()
- thermo_utils: fix stored-reaction deltaG bug (read reaction_obj.delta_g to
  match its guard instead of reaction_obj.deltag)
- tests: rewrite stale test_ms_biochem_deltag.py against the relocated deltaG
  API in ThermoUtils; add gated molGPK predominant-ion pH-sweep test and
  units-field assertions
- environment.yml: pinned conda env covering rdkit/openbabel/torch and backend
  runtime deps for a reproducible one-env install
@VibhavSetlur
VibhavSetlur force-pushed the feature/verab-methoxy-pickaxe branch from e235519 to 50fcf36 Compare July 22, 2026 16:29
…ness fixes (FIX1-5)

FIX1: verab defaults to mechanism-informed operators (Pate 2026), runtime-resolved via config/env/DependencyManager from stefanpate/coarse-grain-rxns (not vendored, no license); honest fallback to bundled metacyc_intermediate; records rule set actually used. FIX2: preserve per-reaction firing-operator attribution as a list. FIX3: screening distinguishes lookup failure from not-found. FIX4: genome prediction is multi-EC + explainable. FIX5: consolidated installable KING bundle to king_app/verab/; manifest records rule_set_used.
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