Add cofactor-EC coverage as a second, independent cognate-ligand source - #11
Draft
m-crown wants to merge 6 commits into
Draft
Add cofactor-EC coverage as a second, independent cognate-ligand source#11m-crown wants to merge 6 commits into
m-crown wants to merge 6 commits into
Conversation
Combines CoFactor DB 2010 + BRENDA (vendored from RelLig's cofactor_ec.csv and cofactors_details.json) with UniProt's structured COFACTOR annotations into cofactor_ligands_df.pkl (preprocess_cofactors.py, mirrors preprocess_rhea.py's pattern), then concats it into get_ec_information.py's cognate_ligands_df alongside the existing reaction-derived sources, tagged via a new ligand_source column. Partial ECs are handled per the plan doc's broadcast rule: subsubclass-level wildcards (N.N.N.-) are expanded to matching terminal ECs, class/subclass-level wildcards are dropped rather than broadcast (verified unsafe - would match every terminal EC in that class). See docs/cofactor_coverage_plan.md for the full design and the coverage numbers behind it. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…ssing Found while running preprocess_rhea.py end-to-end for the first time against real reference data (unrelated to the cofactor-coverage work, but this script is a prerequisite for it). None of these are cofactor- specific: - rdfreader's .metadata property fixed-column-parses a mol block's whole header, including the program/timestamp line - Rhea's CDK-written molfiles don't conform to that fixed-width format (e.g. " CDK 2/12/10,15:27" instead of MDL's fixed MMDDYY field), which broke metadata parsing on ~97% of real Rhea rd/ files. Read the molecule name directly from the mol block's raw first line instead, since that's all this script actually needed from .metadata. - Assigning int -1 into a pandas StringDtype column raises under modern pandas (silently allowed under older/object-dtype behavior). Assign the string "-1" instead - the column gets cast to int on the next line regardless. - ChEBI's names.tsv.gz moved to lowercase columns upstream (already documented in download_reference_data.py's derive_chebi_results docstring, but never applied here) - rename after reading. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…Chem fix Found and fixed while running the full cognate-ligand build end-to-end for the first time: - None of the module's live requests.get/post calls (KEGG, PubChem, GlyTouCan) had a timeout, so a stalled/dead connection (e.g. the host machine sleeping mid-request) hung forever. Added request_with_retry: a (connect, read) timeout tuple - not a single float, since a single combined timeout was observed to not reliably bound a stuck TCP handshake (socket stuck in SYN_SENT) - plus a small retry-with-backoff, wired into all 6 call sites. - Every live-fetch loop (PubChem block, KEGG-direct fallback records and SMILES, glycan compound records, GlyTouCan resolution) was completely silent between a "Getting X records" print and the next checkpoint, sometimes covering thousands of calls - a genuine stall and normal-but- slow progress looked identical from the log. Added fetch_all_with_progress() (timestamped, with an ETA) and used it everywhere; added timestamps to the two chunk-fetch loops that already printed. - The PubChem block was making one live KEGG call per code across the *entire* ~26k-code bulk PubChem<->KEGG cross-reference file, not filtered to the ~8.2k codes actually relevant to this run's reactions - the dominant, previously-invisible cost of that phase. Filtered to relevant codes before the live fetch. - PubChem's PUG-REST API has renamed/deprecated the "CanonicalSMILES" property - confirmed live that requesting it now returns the value under a "ConnectivitySMILES" JSON key instead. Request and read the current key name throughout. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Found running the full cognate-ligand build end-to-end for the first time against real, current data: - ChEBI's relation.tsv moved from a string TYPE column to a numeric relation_type_id foreign key (INIT_ID/FINAL_ID -> lowercase init_id/final_id too). Resolve "has_role" via the small relation_type.tsv vocabulary file (new --chebi_relation_types arg, new chebi_relation_types manifest entry) instead of hardcoding a number. Also found and fixed the relation direction itself: init_id is the specific compound and final_id is the role class (glucose --has_role--> nutrient), not the reverse an earlier version of this code assumed - confirmed empirically, real data has zero rows the old direction would have matched. - Two more instances of the same pandas strict-dtype-coercion issue already fixed once today in preprocess_rhea.py: assigning string values into a column pandas inferred as float64/int64 (the glycan offline-fallback "smiles" column, and the isCofactor "TYPE" column built from relation_type_id) raises under this pandas version. Fixed by casting to object dtype before the assignment; added a regression test reproducing the real-world float64-inference case (a literal Python list of str/NaN infers as object dtype and never reproduced this - only .apply() returning all-NaN does). - PubChem's PUG-REST API deprecated the "CanonicalSMILES" property in favour of "ConnectivitySMILES" - updated the request and all downstream column references. Also request PubChem's Title property in the same already-batched call and use it as the compound_name for PubChem-sourced ligands, replacing a separate per-KEGG-code live name lookup that turned out to be the single most expensive call in the whole build (~2.5h for ~6.5k codes) for a purely cosmetic field. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Ran the full cognate-ligand pipeline end-to-end and benchmarked the result against the last pre-cofactor-coverage cognate_ligands_df.pkl (2024-07-18 - checksum-identical to the file behind the published Zenodo v1-0-2 flat files and the copy vendored into AlphaCognate): +87% EC-ligand pairs (41,785 -> 78,163), +1,102 net ECs covered (6,214 -> 7,316), with the "strict addition" design intent confirmed (40,282 pairs unchanged) and a small (3.6%), fully-categorised set of differences from the old build. - nextflow/bin/benchmark_cognate_ligands.py: repeatable comparison script (row/EC/ligand counts, pair-level diff, lost-pair categorisation, isCofactor distribution) - re-run after any future cognate-ligand pipeline change instead of assuming its effect. - docs/v2_cofactor_coverage.md: first-class, standalone writeup of this work (motivation, method, full results, categorised limitations, reproduction instructions) - written for eventual reuse in a v2 paper/docs section, not just session notes. - docs/cofactor_coverage_plan.md: trimmed the results section down to a pointer at the new doc, avoiding two copies of the same numbers drifting out of sync. - README.md: filled in the long-standing "XXX cognate ligands" placeholder with real v2 numbers, split the Sources list into the reaction-derived (v1) and cofactor-derived (v2) routes, linked to the new doc. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Ran a small, hand-picked (2CPP/1MTY/1JB0) end-to-end validation of the
cofactor-coverage work against real PDB structures, standalone (no
Nextflow/conda - a Python 3.9 venv with pip-installed openbabel-wheel +
pdbe-arpeggio was enough). Confirmed the mechanism works on a real case:
P450cam's bound heme scores 0.717 PARITY similarity against the new
cofactor-sourced HEA cognate ligand (clears the 0.40 threshold; this
ligand exists in cognate_ligands_df.pkl only because of this session's
cofactor work). Also confirmed a real, specific limitation: MMO's bare
Fe3+ ion doesn't match anything, since its real cofactor (a
carboxylate-bridged diiron cluster) isn't captured by either combined
source as this structure deposits it.
Two more bugs found and fixed along the way, same classes as earlier
today:
- process_pdb_contacts.py: db_df["xref_db"] assigning a list into a
strictly-string-typed column (same pandas dtype-coercion issue fixed
repeatedly elsewhere today) - cast to object first.
- process_pdb_contacts.py: np.where(..., "minor", np.nan) - numpy 2.x's
stricter dtype promotion now rejects mixing a string result with
np.nan ("could not be promoted"). Use None instead, which pandas still
treats as missing downstream.
Documented in docs/v2_cofactor_coverage.md ("Small-scale PARITY
validation") and docs/cofactor_coverage_plan.md's next-steps list. A
proper, larger-scale sampled benchmark (not just 2-3 hand-picked
structures) is tracked as a standing memory to-do
(todo_v2_parity_benchmark) rather than left only in these docs, so it
survives across sessions.
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
docs/cofactor_coverage_plan.md: cognate ligands were previously derived only from reaction equations, which structurally misses catalytic/prosthetic cofactors that never appear as a reactant/product (e.g. NAD, FAD, PLP, metal centers).cofactor_ec.csv/cofactors_details.json) combined with UniProt's structuredCOFACTORannotations, built intocofactor_ligands_df.pklby newnextflow/bin/preprocess_cofactors.py(mirrorspreprocess_rhea.py's pattern).get_ec_information.pyconcats this alongside the existing reaction-derived sources intocognate_ligands_df, tagged via a newligand_sourcecolumn (reaction/cofactor, unioned when both apply). Strict addition - existing reaction-derived rows are unaffected (verified by test).N.N.N.-) are expanded to matching terminal ECs; class/subclass-level wildcards are dropped rather than broadcast (confirmed during planning that broadcasting those produces nonsense - e.g. a class-level wildcard matching literally every terminal EC in that class).get_chem_comp_descriptorsmoved fromprocess_all_pdb_contacts.pyintoutils.pyso both scripts can share the CCD→SMILES resolution logic.reference_data_manifest.yamlupdated: removed the stale RelLig bulk-TSV entries from an earlier version of the plan (checked against the real files - no EC column in either), added the three sources actually used, all optional and URL-verified.Test plan
nextflow/bin/tests/test_preprocess_cofactors.py- 21 tests covering UniProt parsing, the EC-completeness/broadcast safety rule (including a regression test for the class-level-wildcard bug caught during planning), CoFactor DB table loading, ChEBI resolution, and end-to-endbuild_cofactor_ligands_df.nextflow/bin/tests/test_utils_ec_broadcast.py- 10 tests forclassify_ec_completeness,broadcast_subsubclass_ec, and the movedget_chem_comp_descriptors.nextflow/bin/tests/test_get_ec_information_cofactor_merge.py- 4 tests replicating thecognate_ligands_dfconcat/groupby block to verifyligand_sourceplumbing and the strict-addition invariant.data_dirsetup).Known, explicitly out-of-scope gap: chlorophyll and other free-text-only UniProt
COFACTORannotations (~0.7% of the bulk pull) aren't mechanically resolvable by this path - deferred to a separate LLM-extraction follow-on per the plan doc.🤖 Generated with Claude Code