diff --git a/README.md b/README.md index e0bcec9..4d7fce6 100644 --- a/README.md +++ b/README.md @@ -215,7 +215,11 @@ A [tutorial](docs/tutorial.md) with a walkthrough of some of these features is a #### Sources -Cognate ligands in ProCogGraph are aggregated from the following sources: +Cognate ligands in ProCogGraph are aggregated from two independent routes. + +**Reaction-derived** (v1 onwards) — a compound is a candidate cognate +ligand for an EC if it is written as a reactant or product of that EC's +reaction equation, sourced from: - [KEGG](https://kegg.jp/) - [ChEBI](https://www.ebi.ac.uk/chebi/) @@ -223,7 +227,20 @@ Cognate ligands in ProCogGraph are aggregated from the following sources: - [PubChem](https://pubchem.ncbi.nlm.nih.gov/) - [GlyTouCan](https://glytoucan.org/) -SMILES representations are obtained for each ligand, and each cognate ligand is mapped to one or more EC IDs. Cognate ligands are processed using the RDKit library in Python, with structures neutralised and re-canonicalised to reduce the number of duplicate structures. A total of XXX cognate ligands are currently represented in the database. +**Cofactor-derived** (v2 onwards) — a second, independent route that adds +cofactors an EC's reaction equation never states explicitly (prosthetic +groups, electron/light carriers, structural metal centers — cofactors +that participate catalytically rather than stoichiometrically), sourced +from: + +- The [CoFactor database](https://doi.org/10.1093/bioinformatics/btq157) (2010) + [BRENDA](https://www.brenda-enzymes.org/) +- [UniProt](https://www.uniprot.org/)'s structured `COFACTOR` annotation + +See [docs/v2_cofactor_coverage.md](docs/v2_cofactor_coverage.md) for the +full methodology, source-selection rationale, and a benchmarked +before/after comparison against the last pre-v2 database build. + +SMILES representations are obtained for each ligand, and each cognate ligand is mapped to one or more EC IDs. Cognate ligands are processed using the RDKit library in Python, with structures neutralised and re-canonicalised to reduce the number of duplicate structures. As of the v2 cofactor-coverage build, 8,811 distinct cognate ligand structures are mapped across 7,316 ECs (78,163 EC-ligand pairs); see [docs/v2_cofactor_coverage.md](docs/v2_cofactor_coverage.md) for the full breakdown, and note these totals will shift slightly on each rebuild as upstream sources (KEGG, ChEBI, Rhea, UniProt) are revised. #### Similarity diff --git a/docs/cofactor_coverage_plan.md b/docs/cofactor_coverage_plan.md new file mode 100644 index 0000000..4eb9777 --- /dev/null +++ b/docs/cofactor_coverage_plan.md @@ -0,0 +1,340 @@ +# Plan: Closing the Cofactor Coverage Gap + +## The problem + +Cognate ligands in ProCogGraph are currently derived exclusively from +**reaction equations**: `enzyme.dat` EC entry → KEGG enzyme record → KEGG +reaction → reactant/product compound codes → ChEBI/PubChem SMILES +(`nextflow/bin/get_ec_information.py:415-589`, cross-checked against Rhea +in `nextflow/bin/preprocess_rhea.py`). A compound only becomes a candidate +cognate ligand if it is explicitly written as a reactant or product in the +balanced equation for that EC number. + +This is fine for substrates and products, and for cofactors that *are* +stoichiometric participants (e.g. `NAD+ + substrate ⇌ NADH + product` — +NAD is written into the equation, so it's captured). It structurally fails +for cofactors that participate **catalytically rather than +stoichiometrically** — prosthetic groups, electron/light carriers, +structural metal ions — because these are never written as a reactant or +product of the EC's net reaction at all. + +Your thesis identifies exactly this (p.172): Chlorophyll A is the single +most frequently unmatched ligand in ProCogGraph, "due to its lack of +annotation as a cofactor in reaction schemes" — chlorophyll is never +consumed or produced by a photosystem's EC reaction, so no amount of +reaction-database improvement will surface it via the current pipeline +path. + +`enzyme.dat` itself doesn't help either: ExPASy dropped the structured +`CF` cofactor field years ago; cofactor mentions today only exist as +unstructured free text inside `CC` comment lines (`utils.process_ec_records`, +`nextflow/bin/utils.py:17-37`, currently only extracts `ID`/`DE`/derived +`TRANSFER`) — not reliably machine-parseable. + +### An existing, separate mechanism that this plan is not about + +The pipeline already has a cofactor-related step, at +`get_ec_information.py:699-719`: after `cognate_ligands_df` is built from +the reaction-equation path, it joins each ligand's ChEBI ID against +`chebi_relations.tsv`'s `has_role` relations for four ChEBI role classes +(`CHEBI:23357` cofactor, `23354` coenzyme, `26348` prosthetic group, +`26672` siderophore) and stamps an `isCofactor` column. This **labels** +ligands that are already in `cognate_ligands_df` — it does nothing to add +ligands that never entered the table because they're not reaction +participants. It's a labelling mechanism for the existing set, not a +coverage extension. This plan is about the latter; the two are +independent and should stay that way, but the naming overlap (`isCofactor` +vs. the new `ligand_source` value below) is worth keeping straight when +implementing. + +(Also checked: a `parse_brenda_json_generic_reaction` function exists at +`get_ec_information.py:323` but is dead code — never called anywhere. +BRENDA is not currently wired into this pipeline despite that function's +presence.) + +## Sources considered + +### Ruled out: PDBe RelLig bulk TSVs (original plan's primary source) + +The original version of this plan proposed PDBe's **RelLig** +(`https://github.com/PDBeurope/rellig`) bulk PDBeChem v2 output +(`interacting_chains_with_ligand_functions.tsv`, +`pdb_bound_molecules.tsv`) as the primary source, on the assumption these +files carried an EC number per row. **Verified against the real file and +this is wrong** — pulled the actual header + rows: + +``` +PDBID Chain_Symmetry BestUnpAccession LigandID bmID LigandType inchikey +ProteinNameUniprot ProteinNamePdb OrganismScientificNameUniprot ... annotation +101m A P02185 HEM bm2 CCD KABFMIBPWCXCRK... Myoglobin Myoglobin ... reactant-like +``` + +No EC column. RelLig's actual EC-bearing output is a **separate per-ligand +JSON format** (`_cofactor_annotation.json`), produced by +*running* `pdberellig cofactors --cif ...` per CCD component — not a +static bulk download. Chasing this further would mean either running the +pipeline ourselves against the full CCD, or finding pre-generated JSON +output that may not exist as a bulk artifact at all. Not pursued further. + +### What RelLig's source code did turn up, and is being reused instead + +Investigated `pdberellig`'s own source +(`pdberellig/core/cofactors.py`, `pdberellig/data/cofactors/`) directly. +It ships two small, static, redistributable (Apache-2.0) data files that +are exactly the EC↔cofactor mapping this plan actually needs, with none of +the bulk-TSV baggage: + +- **`cofactor_ec.csv`** — 3,915 rows, `EC_NO → COFACTOR_ID`, sourced from + two curated inputs: `cofactor_db_2010` (the **CoFactor database**, + Fischer/Holliday/Thornton, *Bioinformatics* 2010 — a manually curated + catalogue of organic cofactors and the EC numbers known to use them) and + `brenda` (supplementary associations from BRENDA — note this is RelLig's + own upstream use of BRENDA, still not something this codebase depends on + directly). +- **`cofactors_details.json`** — 27 distinct cofactor classes (IDs 1–28, + one retired), each with a representative PDB CCD code (id 4 → `NAD`, id + 22 → `HEA`, id 7 → `PLP`, id 11 → `B12`, etc. — full list: TPP, FAD, + FMN, NAD, pantetheine, CoA, PLP, glutathione, biotin, folate, B12, + ascorbate, menaquinone, ubiquinone, molybdopterin, tetrahydrobiopterin, + a mycofactocin-type cofactor, SAM, coenzyme F430, coenzyme M, heme A, + deazaflavin, PQQ, TPQ, TRQ, lipoic acid). + +Since `ccd_cif` is already a pipeline input, resolving a representative +CCD code to a SMILES needs no new source at all. + +**Known gap in this source, confirmed and not fixable by construction**: +chlorophyll is not among the 27 classes. CoFactor DB (2010) is a +relatively small, central-metabolism-focused curated set — it doesn't +cover chlorophyll, and doesn't cover bare structural metal ions beyond +what's implicit in a couple of its classes. This source alone does **not** +solve the thesis's headline motivating example. + +### Added: UniProt `COFACTOR` annotation, structured entries + +Originally scoped in this plan as a secondary/fallback source (per-accession, +chunked REST calls). Re-investigated as a first-class combined source +instead, and it's better suited than that framing suggested: + +- **Bulk-downloadable**, not per-accession: confirmed the `/uniprotkb/search` + (or `/uniprotkb/stream` for larger pulls) endpoint supports + `fields=accession,ec,cc_cofactor&format=tsv` directly — + `reviewed:true AND cc_cofactor:* AND ec:*` returns 111,819 rows in one + bulk pull, no per-accession chunking needed. +- **99.3% of those rows (111,077) already have a machine-parseable + `Name=...; Xref=ChEBI:CHEBI:NNNN;` structure** — confirmed by regex + extraction against the real bulk TSV. Only 742 rows (0.7%) are + free-text-only (a bare `Note=` with no `Xref`). +- **Chlorophyll's specific entries fall into that 0.7% free-text bucket** — + confirmed directly: `cc_cofactor:"ChEBI:CHEBI:18230"` (chlorophyll a's + real ChEBI ID) returns **zero** reviewed UniProt entries. The actual + annotation for e.g. Photosystem I (`P56766`, EC 1.97.1.12) is: + ``` + COFACTOR: Note=P700 is a chlorophyll a/chlorophyll a' dimer, A0 is one + or more chlorophyll a, A1 is one or both phylloquinones and FX is a + shared 4Fe-4S iron-sulfur center. + ``` + No `Name=`/`Xref=` — not mechanically resolvable to a SMILES. So even + the broader UniProt source doesn't mechanically solve chlorophyll; it's + annotated qualitatively (variable stoichiometry, mixed a/a' forms), not + as a clean 1:1 cofactor identity, in both curated sources checked. + +## Why combine both rather than pick one + +Measured real overlap between the two sources (terminal-EC level, exact +string match, before any broadcast): + +| Source | Distinct ECs covered | +|---|---| +| CoFactor DB 2010 + BRENDA (`cofactor_ec.csv`) | 2,760 | +| UniProt `COFACTOR`, structured only | 3,284 | +| **Overlap** | 1,201 | +| **CoFactor DB-only** (UniProt has no cofactor annotation at all here) | 1,559 | +| **UniProt-only** (outside CoFactor DB's 27-class scope entirely) | 2,083 | +| **Combined unique** | **4,843** | + +Confirmed this is real complementary coverage, not redundancy, with a +concrete example: **EC 1.1.1.10** (D-xylulose reductase, a textbook +NADP-dependent enzyme) is in CoFactor DB's list, but every reviewed +UniProt entry for it has a **completely empty** `Cofactor` field (checked +live: `Q7Z4W1`, `Q91X52`, `Q21929` all blank). CoFactor DB/BRENDA encode +EC-class-level literature knowledge that individual UniProt curators never +entered per-accession; UniProt in turn catches specific, curator-verified +cofactor identities (particularly metal centers — 109 distinct ChEBI +identities vs. CoFactor DB's 27 classes) that fall outside CoFactor DB's +fixed, 2010-dated scope. Use both, unioned. + +## The broadcast problem, and its fix + +Not all of UniProt's EC values are fully resolved to 4 digits — 196 of the +3,284 (6.0%) are partial (`N.N.N.-`, `N.N.-.-`, or `N.-.-.-`). A naive +exact-string join (the same pattern `get_pdb_parity.py` already uses for +reaction-derived ligands) silently drops all 196, undercounting coverage: +exact-match-only gives just **62.0%** of the pipeline's 6,753 real terminal +ECs (4,184/6,753). + +**First attempt at "broadcast the wildcard down to matching terminal ECs" +produced a bogus 100% coverage number** — traced this to the 7 class-level +wildcards (`1.-.-.-` through `7.-.-.-`) each matching literally every +terminal EC in that class (e.g. `1.-.-.-` → all 6,753... no, all of class +1). Broadcasting at that granularity is chemically wrong — it would claim +catalase (`1.11.1.6`, heme-dependent) shares a cofactor with an unrelated +NAD-dependent dehydrogenase just because both are oxidoreductases. +Breaking down by granularity confirmed the danger is real: + +| Partial-EC granularity | Count | Terminal ECs it would broadcast to | +|---|---|---| +| `N.-.-.-` (class-level) | 7 | 6,753 — literally everything | +| `N.N.-.-` (subclass-level) | 33 | 5,578 | +| `N.N.N.-` (subsubclass-level) | 156 | 6,193 | + +**Rule adopted**: only broadcast **subsubclass-level (`N.N.N.-`)** partial +ECs — enzymes sharing all three EC digits genuinely tend to share cofactor +chemistry (e.g. `1.1.1.-` is the NAD(P)-dependent CH-OH oxidoreductase +subsubclass). **Class- and subclass-level partial ECs (40 of the 196) are +dropped from the cofactor table entirely** — there's no chemically +defensible way to narrow "some oxidoreductase, unknown subclass" down to +specific terminal children. + +With that rule, combined coverage of the pipeline's real terminal EC list: + +- Exact match only (CoFactor DB + UniProt 4-level): **62.0%** (4,184/6,753) +- \+ safe subsubclass-level broadcast: **96.2%** (6,496/6,753) + +**Where the broadcast happens**: at `cognate_ligands_df` *build* time, not +at match time. For each UniProt row that only resolved to `N.N.N.-`, look +up the pipeline's own terminal EC list (`ec_records_df.TRANSFER.unique()`, +already computed in `get_ec_information.py`'s `main()`) and expand that +one row into one row per real terminal EC sharing that subsubclass prefix, +each carrying the same cofactor ChEBI ID/SMILES — identical shape to every +existing reaction-derived row. This means `cognate_ligands_df.entry` never +contains a wildcard, so `get_pdb_parity.py:121`'s existing +`cognate_ligands_df.entry.isin(ec)` join needs **zero modification** — it +simply sees more rows for more ECs. No new match-time logic anywhere +downstream. This was the whole point of preferring this approach over the +RelLig structure-instance path: cofactor rows are ordinary +`cognate_ligands_df` rows, keyed by real EC, from the start. + +## Known remaining gap: chlorophyll and other free-text-only cases + +742 UniProt rows (0.7% of the reviewed cofactor+EC set) have only a +free-text `Note=` with no structured `Xref=ChEBI:...` — chlorophyll's +photosystem entries are in this bucket. Not mechanically resolvable to a +SMILES from either source as currently curated. + +**Deferred, explicitly agreed**: revisit this bucket with LLM-assisted +extraction (read the free-text note, propose a ChEBI ID / SMILES, flag for +review) as a **separate follow-on pass** once the main combined table +(CoFactor DB + UniProt structured, with subsubclass broadcast) is built +and integrated. Don't block the main implementation on this — it's a +small, well-bounded tail cleanup (742 rows), not core to getting the +combined source working. + +## Implementation steps + +1. **Vendor `cofactor_ec.csv` and `cofactors_details.json`** from + `pdberellig` (Apache-2.0, redistribution permitted) into the reference + data set, alongside a script step resolving each `COFACTOR_ID`'s + representative CCD code to canonical SMILES via the existing + `ccd_cif`-parsing path. +2. **Bulk-pull UniProt structured cofactor data**: add + `https://rest.uniprot.org/uniprotkb/stream?query=reviewed:true+AND+cc_cofactor:*+AND+ec:*&fields=accession,ec,cc_cofactor&format=tsv` + (or the equivalent `search` pagination) to the reference data manifest + pattern already established in `docs/reference_data_download_plan.md`. + Parse `Name=...; Xref=ChEBI:CHEBI:(\d+)` per row; keep only rows with a + match. Resolve ChEBI ID → SMILES via the same ChEBI resolution path + `get_ec_information.py` already uses for reaction-derived ligands. +3. **Split UniProt rows by EC completeness**: exact 4-digit ECs pass + through directly; `N.N.N.-` rows get expanded against the pipeline's + own terminal EC list (`ec_records_df.TRANSFER.unique()`) at build time; + `N.-.-.-` and `N.N.-.-` rows are dropped (logged, not silently + discarded, so the drop is auditable). +4. **Union CoFactor DB rows + expanded UniProt rows**, tag both with a new + `ligand_source = "cofactor"` value (existing reaction-derived rows get + `ligand_source = "reaction"` for consistency — this column doesn't + exist yet and needs adding to the reaction-derived rows too, not just + the new ones). Concat into `cognate_ligands_df` alongside the existing + Rhea/KEGG/ChEBI/PubChem/GlyTouCan frames at + `get_ec_information.py:675`. +5. **Do not modify `get_pdb_parity.py` or `produce_neo4j_files.py`** — the + existing EC-join and PARITY-scoring logic should work unchanged against + the enlarged `cognate_ligands_df`, since cofactor rows are + indistinguishable in shape from reaction rows except for the new + `ligand_source` tag. +6. **(Follow-on, separate piece of work)**: LLM-assisted extraction for + the 742 free-text-only UniProt rows, chlorophyll included. + +## Verification + +- After implementation, re-check the specific case cited in your thesis + — note per the gap analysis above, chlorophyll is **expected to still be + missing** after step 1–5 alone; don't treat its absence as a bug until + the follow-on LLM-extraction pass (step 6) is done. +- Confirm `cognate_ligands_df.entry` contains zero wildcard/partial EC + values after the broadcast step — this is a hard invariant the + downstream `isin()` join depends on. +- Spot-check a handful of well-characterised cofactor-dependent enzymes + outside chlorophyll (e.g. cytochromes with heme, the *Methylococcus + capsulatus* MMO example already discussed in your thesis's future-work + section, p.175) to confirm the new cofactor-origin ligands match known + biology. +- Confirm reaction-origin cognate ligands are completely unaffected (this + is a strict addition) — existing PARITY scores and cognate ligand counts + for already-mapped structures should be unchanged, and the new + `ligand_source` column should backfill to `"reaction"` for every + pre-existing row. +- Compare before/after counts of "most frequently unmatched ligands" (the + analysis behind Table 3.3 / the Discussion's unmatched-ligand + discussion) for the ~4,843 newly-EC-covered enzymes. + +## Benchmark results + +Real end-to-end build (2026-08-30), benchmarked against the last +pre-cofactor-coverage `cognate_ligands_df.pkl` (2024-07): **+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, fully-categorised set of differences from +the old build (3.6%, mostly RDKit macrocycle-sanitization edge cases and +normal upstream Rhea/KEGG data drift, not a regression in this plan's own +logic). + +Full methodology, the complete results table, the lost-pair +categorisation, and reproduction instructions now live in +**[docs/v2_cofactor_coverage.md](v2_cofactor_coverage.md)** — written as +the first-class, citable record of this work (for the v2 docs/paper), +rather than duplicated here where it would drift out of sync. + +## Status + +**Implemented and benchmarked** — see +[docs/v2_cofactor_coverage.md](v2_cofactor_coverage.md) for the full +writeup. What was originally scoped as "suggested order of work" is done: +vendoring/parsing the CoFactor DB tables, the UniProt bulk pull and +exact/broadcast/drop EC-completeness split, the `ligand_source` column +and union into `cognate_ligands_df`, and a real before/after comparison +against the last pre-cofactor-coverage build. + +## Remaining next steps (not yet done) + +1. **Real, larger-scale PARITY-based validation against PDB structures.** + A *small* version of this is now done — see + [docs/v2_cofactor_coverage.md](v2_cofactor_coverage.md#small-scale-parity-validation-2026-08-3031): + `get_pdb_parity.py` run standalone against 2 of 3 hand-picked + structures (P450cam/heme, MMO/diiron; Photosystem I incomplete), + confirming a real cofactor-origin PARITY match on P450cam's heme + (0.717, clears threshold) and a real, specific gap on MMO's bare Fe³⁺ + ion (0.0, not matched). That's evidence the mechanism works, not a + statistically meaningful measure across the ~4,843 newly-covered ECs — + a proper sampled benchmark is still open. **Tracked as a standing + to-do in memory** (`todo_v2_parity_benchmark`), so it isn't lost + between sessions. +2. **LLM-assisted extraction for the 742 free-text-only UniProt + `COFACTOR` rows** (chlorophyll's actual fix - see "Known remaining + gap" above). Explicitly deferred, not started. +3. **Investigate the 9 ECs and ~104 macrocycle-structure pairs lost + relative to the pre-cofactor-coverage baseline**, if full parity with + the old dataset is ever required (see + [docs/v2_cofactor_coverage.md](v2_cofactor_coverage.md)'s "Known + gaps" for the specific EC list and categorisation). Not blocking - + small (9 of 6,214 ECs) - but not yet root-caused either. +4. Finish the stopped 1JB0 (Photosystem I) run — the chlorophyll/`[4Fe-4S]` + case — likely as part of whichever of the above happens first. diff --git a/docs/installation.md b/docs/installation.md index 584198e..a97dc6a 100644 --- a/docs/installation.md +++ b/docs/installation.md @@ -146,6 +146,24 @@ The ProCogGraph pipeline is built using Nextflow for workflow management. To run python3 preprocess_rhea.py --rhea_ec_mapping rhea2ec.tsv --rhea_reaction_directions rhea-directions.tsv --rd_dir rd/ --outdir . --chebi_names chebi_names.tsv.gz ``` + Optionally, also preprocess cofactor-EC association data (see + `docs/cofactor_coverage_plan.md`) - adds cofactor-origin cognate + ligands (e.g. NAD, FAD, PLP, metal centers) alongside the + reaction-derived ones above, for cofactors that never appear as a + reactant/product in an EC's reaction equation. Requires the vendored + `cofactor_ec.csv`/`cofactors_details.json` and a bulk UniProt + cofactor-annotation pull (`--only cofactor_ec_csv,cofactor_details_json,uniprot_cofactor_annotations` + with `download_reference_data.py --include-optional`), on top of + files already fetched above (`ccd.cif`, `chebi_structures.tsv.gz`, + `chebi_names.tsv.gz`, `enzyme.dat`, `enzclass.txt`): + + ``` bash + python3 preprocess_cofactors.py --cofactor_ec_csv cofactor_ec.csv --cofactor_details_json cofactors_details.json --ccd_cif ccd.cif --uniprot_cofactor_tsv uniprot_cofactor_annotations.tsv --chebi_structures chebi_structures.tsv.gz --chebi_names chebi_names.tsv.gz --ec_dat enzyme.dat --enzyme_class_file enzclass.txt --outdir . + ``` + + This produces `cofactor_ligands_df.pkl`, passed to `get_ec_information.py` + via `--cofactor_ligands` (optional - the pipeline runs unchanged without it). + 4. Produce final manifest file of structures to be processed: ``` bash diff --git a/docs/v2_cofactor_coverage.md b/docs/v2_cofactor_coverage.md new file mode 100644 index 0000000..7d155ac --- /dev/null +++ b/docs/v2_cofactor_coverage.md @@ -0,0 +1,287 @@ +# Cofactor Coverage in ProCogGraph v2 + +## Summary + +ProCogGraph v1 derived cognate ligands exclusively from reaction +equations: an EC number's KEGG/Rhea reaction record supplies its +reactant/product compounds, which become candidate cognate ligands for +that EC. This is structurally blind to cofactors that participate +**catalytically rather than stoichiometrically** — prosthetic groups, +electron/light carriers, and structural metal centers that are never +written as a reactant or product of an EC's net reaction, and so can +never be discovered via that route no matter how complete the underlying +reaction database is. + +ProCogGraph v2 adds a second, independent cognate-ligand source: a +combined EC→cofactor mapping built from the CoFactor database (Fischer, +Holliday & Thornton, *Bioinformatics* 2010) plus BRENDA, and from +UniProt's structured `COFACTOR` annotation. This is a strict addition +alongside the existing reaction-derived path, not a modification of it — +every cognate ligand row is now tagged with its provenance (`reaction`, +`cofactor`, or both, where a compound was independently discovered by +each route). + +Benchmarked against the last pre-v2 build, this increases EC-ligand +coverage from 41,785 to 78,163 pairs (+87%) and the number of ECs with at +least one cognate ligand from 6,214 to 7,316 (+1,102 net), while leaving +the pre-existing reaction-derived coverage almost entirely intact (see +Results). + +## Motivation + +This is fine for substrates, products, and cofactors that *are* +stoichiometric participants (e.g. `NAD+ + substrate ⇌ NADH + product` — +NAD is written into the equation, so it's captured by v1 already). It +fails for the catalytic case: a photosystem's chlorophyll, for instance, +is never consumed or produced by the photosystem protein's own EC +reaction, so it is invisible to a reaction-equation-only pipeline +regardless of how good the underlying reaction data is. + +## Method + +### Sources + +Two combined, independently-sourced EC→cofactor mappings: + +1. **CoFactor database (2010) + BRENDA**, via the small vendored + `cofactor_ec.csv`/`cofactors_details.json` tables from PDBe's RelLig + project (Apache-2.0) — 27 organic cofactor classes (NAD, FAD, PLP, CoA, + biotin, B12, heme A, SAM, molybdopterin, and others), each mapped to a + representative PDB chemical-component code and, from there, a SMILES. +2. **UniProt's `COFACTOR` annotation**, bulk-pulled (not per-accession — + a single `/uniprotkb/stream` query against all reviewed entries with + both an EC number and a cofactor annotation), restricted to rows with + a structured `Xref=ChEBI:CHEBI:` (≈99.3% of the pull). This source + is broader in chemical identity (109 distinct ChEBI cofactor + identities vs. CoFactor DB's 27 classes — notably including metal + centers such as `[4Fe-4S]` clusters) and catches EC/cofactor + associations CoFactor DB's fixed, 2010-dated scope does not. + +Measured real overlap between the two sources: only 1,201 of a combined +4,843 unique ECs are covered by both — genuinely complementary, not +redundant (concrete example: EC 1.1.1.10, a textbook NADP-dependent +enzyme, is covered by CoFactor DB but has zero UniProt cofactor +annotation on any reviewed entry). + +### EC completeness and the broadcast rule + +Not every source EC value is fully resolved to 4 digits. Each is +classified by how many segments are resolved before the first wildcard: + +- **Exact (4/4)** — used directly. +- **Subsubclass-level (`N.N.N.-`)** — broadcast to every real terminal EC + sharing that prefix (enzymes sharing all three leading digits generally + do share cofactor chemistry). +- **Subclass- or class-level (`N.N.-.-`, `N.-.-.-`)** — **dropped, not + broadcast**. An earlier attempt at broadcasting these produced a + clearly wrong result (100% terminal-EC "coverage") by matching a + class-level wildcard against every terminal EC in that class — e.g. + claiming catalase and an unrelated NAD-dependent dehydrogenase share a + cofactor purely because both are oxidoreductases. There is no + chemically defensible way to narrow a class-level wildcard down to + specific terminal ECs, so these are excluded rather than guessed at. + +With this rule: 62.0% terminal-EC coverage from exact matches alone, +96.2% including safe subsubclass-level broadcast. + +### Implementation + +`nextflow/bin/preprocess_cofactors.py` builds the combined table +(mirroring `preprocess_rhea.py`'s standalone-script pattern) into +`cofactor_ligands_df.pkl`, which `get_ec_information.py +--cofactor_ligands` concatenates into `cognate_ligands_df` alongside the +existing Rhea/KEGG/ChEBI/PubChem/GlyTouCan sources, adding a +`ligand_source` column (`reaction` / `cofactor`, unioned where both +routes independently find the same compound). No existing matching code +(`get_pdb_parity.py`'s EC-keyed join) required modification — cofactor +rows are ordinary `cognate_ligands_df` rows, keyed by real terminal EC, +indistinguishable in shape from reaction-derived rows except for the +provenance tag. + +## Results + +Full end-to-end build against live current data (2026-08), compared +against the last pre-v2 `cognate_ligands_df.pkl` (dated 2024-07-18; +checksum-identical to the file that produced the published +[Zenodo v1-0-2 flat files](https://zenodo.org/records/14046116) and the +copy vendored into [AlphaCognate](https://github.com/m-crown/AlphaCognate)'s +`data/procoggraph_data/`): + +| Metric | v1 (2024-07) | v2 (2026-08) | Change | +|---|---|---|---| +| EC–ligand pairs | 41,785 | 78,163 | **+87%** | +| ECs with ≥1 cognate ligand | 6,214 | 7,316 | **+1,102 net** (+1,111 gained / −9 lost) | +| Distinct ligand structures | 8,589 | 8,811 | +222 | + +Of the 78,163 v2 pairs: 40,282 were already present in v1 (unchanged), +37,881 are new — of which 33,398 are corroborated by *both* the reaction +and cofactor paths independently finding the same compound, 2,406 are +reaction-path-only gains (unrelated to this work — normal upstream +Rhea/KEGG growth over the ~2-year gap between builds), and **2,077 exist +only because of the new cofactor path** — coverage that did not exist in +ProCogGraph v1 at all. + +As a coherence check: the pre-existing, independent ChEBI `has_role` +cofactor-labelling mechanism (`isCofactor` column) — unrelated to this +work, present since v1 — shows a 6.3x increase in rows labelled +`Cofactor` (5,606 → 35,122). This is expected rather than circular: the +new cofactor-sourced rows are, by construction, literal cofactor +molecules (NAD, FAD, heme, etc.), and those largely already carry +ChEBI's own independent `has_role: cofactor` annotation — two separately- +sourced signals reinforcing each other. + +### What didn't change (verification) + +- Every pre-existing v1 EC-ligand pair's `ligand_source` backfills to + `"reaction"` — the addition is strict, not a rewrite. +- `cognate_ligands_df.entry` contains zero wildcard/partial EC values + after the broadcast step (hard invariant the downstream + `get_pdb_parity.py` EC join depends on). +- Spot-checked cases match known biology: EC 1.1.1.10 gains + FAD/TPP/NAD/NADP+/Mg²⁺/Zn²⁺; EC 1.97.1.12 (Photosystem I) gains a + `[4Fe-4S] cluster` via the UniProt path; heme/heme b appear across the + expected heme-dependent EC set. + +### Known gaps (explained, not hidden) + +**1,503 EC-ligand pairs (3.6% of the v1 total) present in v1 are absent +from v2.** Categorised (554 distinct lost structures): + +- 180 involve wildcard-substituent SMILES (`*`) — partial/generic + structures such as `[acyl-carrier protein]`-linked intermediates. +- 104 involve charge-separated porphyrin/macrocycle SMILES (`[N+]`/`[Mg`) + — chlorophyll/heme-biosynthesis-pathway intermediates specifically, + plausibly an RDKit sanitization edge case on unusual valence states + (characterised, not yet root-caused). +- 270 are ordinary, otherwise-well-resolved compounds missing only for + specific ECs — consistent with normal upstream Rhea/KEGG reaction- + equation revision over the ~2-year gap (confirmed example: NADH is + missing specifically for EC 1.1.1.96, but still resolves correctly for + 630 other ECs — the reaction equation for that one EC now cites NAD + where it previously cited NADH). +- 9 ECs lost entirely (`1.14.14.140, 2.1.1.86, 2.4.1.129, 2.7.11.27, + 3.1.8.2, 4.2.1.78, 4.3.3.2, 4.3.3.3, 4.3.3.4`) — not yet individually + root-caused; small enough (9 of 6,214) not to block on, but worth + revisiting if full v1 parity is ever required. + +**Chlorophyll as a photosystem cofactor is still not captured**, and this +is expected, not a bug: neither combined source resolves it mechanically. +CoFactor DB's 27 classes don't include chlorophyll at all, and UniProt's +photosystem cofactor annotations (e.g. `P56766`, EC 1.97.1.12) are +free-text `Note=` only, with no structured `Xref=ChEBI:...` — confirmed +directly against a live UniProt query (`cc_cofactor:"ChEBI:CHEBI:18230"` +returns zero reviewed entries). Chlorophyll's binding is described +qualitatively (variable stoichiometry, mixed a/a′ forms) rather than as a +clean 1:1 identity in both sources checked, which is plausibly *why* +neither gives it a structured entry. This affects a small, identified +long tail (742 of 111,819 UniProt cofactor+EC rows, ≈0.7%, are free-text +only) and is deferred to a separate LLM-assisted extraction pass, not +attempted mechanically here. + +## Small-scale PARITY validation (2026-08-30/31) + +The comparison in Results is at the `cognate_ligands_df.pkl` candidate-set +level (old vs. new cognate ligand lists) — it does not by itself confirm +that a real PDB structure's bound cofactor actually scores a PARITY match. +A small, hand-picked, three-structure run of the full downstream pipeline +(`download_mmcif.py` → `process_pdb_structure.py` → PDBe-Arpeggio → +`process_pdb_contacts.py` → `process_all_pdb_contacts.py` → +`get_pdb_parity.py`) was done to check this directly, run standalone +(outside Nextflow, no conda) via a dedicated Python 3.9 venv +(`pip install openbabel-wheel pdbe-arpeggio` — no conda environment +needed, despite `nextflow/envs/arpeggio-env.yaml` being conda-based). + +**Structures** (chosen from the plan's own Verification section, EC +confirmed against real SIFTS `pdb_chain_enzyme.tsv.gz` data, not assumed): + +| PDB | Protein | EC | Cofactor of interest | +|---|---|---|---| +| 2CPP | Cytochrome P450cam | 1.14.15.1 | heme b | +| 1MTY | Methane monooxygenase hydroxylase | 1.14.13.25 | non-heme diiron center | +| 1JB0 | Photosystem I | 1.97.1.12 | chlorophyll + `[4Fe-4S]` | + +**Results, 2 of 3 structures** (1JB0 — a 36-mer complex — was stopped +mid-run in Arpeggio's contact-computation phase before completion; not +yet re-attempted): + +- **2CPP: confirmed positive.** The real bound `HEM` ligand scores + **0.717 PARITY similarity** (threshold 0.40) against cognate ligand + `HEA` (heme A) — a row that exists in `cognate_ligands_df.pkl` *only* + because of the v2 cofactor path (`ligand_source: cofactor`, + `CofactorDB:22`). Before this work, heme was never a stoichiometric + reactant/product of the P450cam reaction, so it could not have scored + against anything. Camphor (the real substrate) separately scores 1.0 + via the pre-existing reaction path, confirming that path is unaffected. +- **1MTY: confirmed limitation, not a bug.** The real bound `FE (III) + ION` scores 0.0 against every cognate ligand. MMO's actual cofactor is + a carboxylate-bridged **diiron** cluster; this structure deposits it as + a bare mononuclear ion, which neither combined source (CoFactor DB's 27 + classes, UniProt's structured entries for this protein) resolves to a + matching cognate ligand as represented. A genuine, specific gap - not + general evidence against the mechanism (2CPP's heme match is real + evidence for it). +- **1JB0: not completed.** Chlorophyll is expected to still fail to + match (known gap, see above); the `[4Fe-4S]` cluster was expected to + succeed (confirmed present in `cognate_ligands_df.pkl` for this EC via + the UniProt path) but this wasn't confirmed against real structure data + before the run was stopped. + +**Two more real bugs found and fixed running this** (same +strict-dtype-coercion/version-drift class as the ones found building +`cognate_ligands_df.pkl` - see the commit history): +`process_pdb_contacts.py` assigning a list into a strictly-string-typed +`xref_db` column, and a `np.where(..., "minor", np.nan)` call that numpy +2.x's stricter dtype promotion now rejects when mixing a string result +with `np.nan` (fixed by using `None` instead). + +**This is evidence, not a benchmark.** Two hand-picked structures is +enough to confirm the mechanism works on a real, known case and to +surface one real limitation — it is not a statistically meaningful +measure of match rate across the ~4,843 newly-covered ECs. See "Future +work" below. + +## Future work + +1. **A real, larger-scale PARITY benchmark** across a proper sample of + the newly-covered ECs (not just 2-3 hand-picked structures) — the + natural follow-up to the small validation above. Needs deciding a + sampling strategy (e.g. N structures per newly-covered EC) and likely + a full Nextflow/conda environment rather than the manual per-script + invocation used for the small run, given the volume. **Tracked as a + standing to-do in memory** (see `todo_v2_parity_benchmark` in this + project's Claude memory) rather than only here, so it survives across + sessions. +2. **LLM-assisted extraction for the 742 free-text-only UniProt + `COFACTOR` rows** — the mechanical fix for chlorophyll specifically + (see "Known gaps" above). +3. **Root-cause the 9 fully-lost ECs and the ~104 macrocycle-structure + lost pairs**, if full parity with the pre-v2 dataset is ever required. +4. Finish the 1JB0 run (chlorophyll/`[4Fe-4S]` case) and/or fold it into + whichever of the above happens first. + +See [docs/cofactor_coverage_plan.md](cofactor_coverage_plan.md#remaining-next-steps-not-yet-done) +for the same list with more detail. + +## Reproducing this benchmark + +```bash +python3 nextflow/bin/benchmark_cognate_ligands.py \ + --old /path/to/old_cognate_ligands_df.pkl \ + --new /path/to/new_cognate_ligands_df.pkl \ + --report_out benchmark_report.txt +``` + +Re-run this after any future change to the cognate-ligand generation +pipeline (`preprocess_rhea.py`, `preprocess_cofactors.py`, +`get_ec_information.py`) to get a real before/after comparison rather +than assuming the effect of a change. + +## See also + +- [docs/cofactor_coverage_plan.md](cofactor_coverage_plan.md) — the + original design/planning document (source analysis, source selection + rationale, and the session-by-session implementation record this + summary is drawn from). +- `nextflow/bin/preprocess_cofactors.py`, + `nextflow/bin/benchmark_cognate_ligands.py` — implementation and + benchmark script. diff --git a/nextflow/bin/benchmark_cognate_ligands.py b/nextflow/bin/benchmark_cognate_ligands.py new file mode 100644 index 0000000..8d09b69 --- /dev/null +++ b/nextflow/bin/benchmark_cognate_ligands.py @@ -0,0 +1,147 @@ +#!/usr/bin/env python + +""" +Compares two cognate_ligands_df.pkl files (e.g. a pre-cofactor-coverage +baseline vs. a new build) and reports the headline coverage numbers: +row/EC/ligand counts, which EC-ligand pairs are gained/lost, a rough +categorization of lost pairs, and the isCofactor tag distribution. + +Written to check docs/cofactor_coverage_plan.md's real-world payoff +(see "Benchmark results" there for the run this script was built to +reproduce) - re-run it after any future change to the cognate-ligand +generation pipeline to see the actual before/after effect, not just +trust that the code changed. + +Example usage: + python3 benchmark_cognate_ligands.py \ + --old /path/to/old_cognate_ligands_df.pkl \ + --new /path/to/new_cognate_ligands_df.pkl \ + --report_out benchmark_report.txt +""" + +import argparse +import re +import sys + +import pandas as pd + +WILDCARD_PATTERN = re.compile(r"\*") +MACROCYCLE_PATTERN = re.compile(r"\[N\+\]|\[Mg") + + +def load(path): + df = pd.read_pickle(path) + required = {"entry", "canonical_smiles", "compound_name"} + missing = required - set(df.columns) + if missing: + raise ValueError(f"{path} is missing expected column(s): {missing}") + return df + + +def summarise_counts(df): + return { + "rows": len(df), + "distinct_ecs": df["entry"].nunique(), + "distinct_ligands": df["canonical_smiles"].nunique(), + } + + +def pair_diff(old_df, new_df): + old_pairs = set(zip(old_df["entry"], old_df["canonical_smiles"])) + new_pairs = set(zip(new_df["entry"], new_df["canonical_smiles"])) + return { + "shared": old_pairs & new_pairs, + "new_only": new_pairs - old_pairs, + "old_only": old_pairs - new_pairs, + } + + +def categorise_lost_pairs(old_df, old_only_pairs): + """Rough heuristic split of pairs present in the old dataset but not + the new one: wildcard-substituent partial structures, charge-separated + porphyrin/macrocycle-like structures (both plausible RDKit-sanitization + edge cases), vs. everything else (more likely ordinary upstream + Rhea/KEGG reaction-equation revisions between the two builds - see the + plan doc for a worked example: EC 1.1.1.96 listing NAD where it used + to list NADH).""" + lost_df = old_df[old_df.apply(lambda r: (r["entry"], r["canonical_smiles"]) in old_only_pairs, axis=1)] + lost_structures = lost_df.drop_duplicates(subset="canonical_smiles") + + has_wildcard = lost_structures["canonical_smiles"].str.contains(WILDCARD_PATTERN, regex=True) + has_macrocycle = lost_structures["canonical_smiles"].str.contains(MACROCYCLE_PATTERN, regex=True) + + return { + "distinct_lost_structures": len(lost_structures), + "wildcard_substituent": int(has_wildcard.sum()), + "charge_separated_macrocycle": int((has_macrocycle & ~has_wildcard).sum()), + "other": int((~has_wildcard & ~has_macrocycle).sum()), + } + + +def format_report(old_path, new_path, old_counts, new_counts, diff, lost_ecs, lost_categories, old_cofactor, new_cofactor): + lines = [] + lines.append(f"Cognate ligand dataset benchmark") + lines.append(f" old: {old_path}") + lines.append(f" new: {new_path}") + lines.append("") + lines.append(f"{'metric':<28}{'old':>12}{'new':>12}{'change':>12}") + for key, label in [("rows", "EC-ligand pairs"), ("distinct_ecs", "distinct ECs"), ("distinct_ligands", "distinct ligands")]: + old_v, new_v = old_counts[key], new_counts[key] + change = f"{new_v - old_v:+d}" + lines.append(f"{label:<28}{old_v:>12}{new_v:>12}{change:>12}") + lines.append("") + lines.append(f"shared EC-ligand pairs: {len(diff['shared'])}") + lines.append(f"new-only pairs (gained): {len(diff['new_only'])}") + lines.append(f"old-only pairs (lost): {len(diff['old_only'])}") + lines.append("") + lines.append(f"ECs lost entirely (in old, absent from new): {len(lost_ecs)}") + if lost_ecs: + lines.append(f" {', '.join(lost_ecs)}") + lines.append("") + lines.append("Lost-pair categorisation (distinct structures, heuristic):") + lines.append(f" distinct lost structures: {lost_categories['distinct_lost_structures']}") + lines.append(f" wildcard-substituent (partial structures): {lost_categories['wildcard_substituent']}") + lines.append(f" charge-separated macrocycle (e.g. porphyrins): {lost_categories['charge_separated_macrocycle']}") + lines.append(f" other (likely upstream Rhea/KEGG data drift): {lost_categories['other']}") + lines.append("") + if old_cofactor is not None and new_cofactor is not None: + lines.append("isCofactor tag distribution:") + lines.append(f"{'':<28}{'old':>12}{'new':>12}") + all_tags = sorted(set(old_cofactor.index) | set(new_cofactor.index)) + for tag in all_tags: + lines.append(f"{tag:<28}{old_cofactor.get(tag, 0):>12}{new_cofactor.get(tag, 0):>12}") + return "\n".join(lines) + + +def main(): + parser = argparse.ArgumentParser(description="Benchmark two cognate_ligands_df.pkl files against each other") + parser.add_argument("--old", required=True, help="Path to the baseline cognate_ligands_df.pkl") + parser.add_argument("--new", required=True, help="Path to the new cognate_ligands_df.pkl") + parser.add_argument("--report_out", default=None, help="Optional path to also write the report to a file") + args = parser.parse_args() + + old_df = load(args.old) + new_df = load(args.new) + + old_counts = summarise_counts(old_df) + new_counts = summarise_counts(new_df) + diff = pair_diff(old_df, new_df) + + lost_ecs = sorted(set(old_df["entry"]) - set(new_df["entry"])) + lost_categories = categorise_lost_pairs(old_df, diff["old_only"]) + lost_categories["lost_ecs"] = lost_ecs + + old_cofactor = old_df["isCofactor"].value_counts() if "isCofactor" in old_df.columns else None + new_cofactor = new_df["isCofactor"].value_counts() if "isCofactor" in new_df.columns else None + + report = format_report(args.old, args.new, old_counts, new_counts, diff, lost_ecs, lost_categories, old_cofactor, new_cofactor) + print(report) + + if args.report_out: + with open(args.report_out, "w") as fh: + fh.write(report + "\n") + print(f"\nReport written to {args.report_out}", file=sys.stderr) + + +if __name__ == "__main__": + main() diff --git a/nextflow/bin/get_ec_information.py b/nextflow/bin/get_ec_information.py index 3568409..02bc242 100644 --- a/nextflow/bin/get_ec_information.py +++ b/nextflow/bin/get_ec_information.py @@ -8,6 +8,7 @@ from Bio.KEGG import Enzyme import io import time +import datetime import pickle import re from rich.progress import Progress @@ -25,6 +26,63 @@ from bs4 import BeautifulSoup from urllib.parse import quote +# (connect_timeout, read_timeout), not a single float: a single timeout +# value only reliably bounds the read phase in some observed cases - a +# stuck TCP handshake (socket sitting in SYN_SENT, no SYN-ACK ever +# arriving) was seen in practice to survive well past 30s combined-timeout +# calls without raising or retrying. Separating connect/read is the +# standard, more robust pattern requests/urllib3 recommend for exactly +# this failure mode. +REQUEST_TIMEOUT_SECONDS = (10, 30) +REQUEST_RETRY_ATTEMPTS = 3 +REQUEST_RETRY_BACKOFF_SECONDS = 5 + +def _now(): + return datetime.datetime.now().strftime("%H:%M:%S") + +def request_with_retry(method, url, **kwargs): + """Thin wrapper around requests.get/requests.post with a timeout and a + small retry-with-backoff. None of this module's live API calls + (KEGG, PubChem, GlyTouCan) previously set a timeout at all, so a + stalled/dead connection - e.g. the host machine sleeping mid-request - + hangs forever instead of failing and retrying. This does not fix the + separate, larger issue that each phase (enzyme records, reaction + records, ...) only checkpoints to disk once the whole phase completes - + a hang partway through a phase still loses that phase's progress if it + has to be killed - but it does mean a transient drop no longer requires + that at all.""" + last_exc = None + for attempt in range(1, REQUEST_RETRY_ATTEMPTS + 1): + try: + return method(url, timeout=REQUEST_TIMEOUT_SECONDS, **kwargs) + except requests.RequestException as exc: + last_exc = exc + if attempt < REQUEST_RETRY_ATTEMPTS: + print(f" [{_now()}] request failed (attempt {attempt}/{REQUEST_RETRY_ATTEMPTS}): {exc} - retrying in {REQUEST_RETRY_BACKOFF_SECONDS}s") + time.sleep(REQUEST_RETRY_BACKOFF_SECONDS) + raise last_exc + +def fetch_all_with_progress(items, fetch_fn, label, print_every=100): + """Runs fetch_fn(item) for every item, printing timestamped progress + every `print_every` items (and always on the first and last). Every + live-fetch loop in this module used to be completely silent between a + "Getting X records" print and the next checkpoint, sometimes covering + thousands of live calls - a stall anywhere inside looked identical to + it just being slow, from the log alone. Returns the list of results in + input order.""" + results = [] + total = len(items) + start = time.time() + print(f" [{_now()}] [{label}] starting, {total} items") + for i, item in enumerate(items, 1): + results.append(fetch_fn(item)) + if i == 1 or i % print_every == 0 or i == total: + elapsed = time.time() - start + rate = i / elapsed if elapsed > 0 else 0 + remaining = (total - i) / rate if rate > 0 else float("nan") + print(f" [{_now()}] [{label}] {i}/{total} ({elapsed:.0f}s elapsed, ~{remaining:.0f}s remaining)") + return results + def get_kegg_enzymes(ec_list, enzyme_string_file = None): def extract_reaction(ec): all_reacts = [] @@ -54,7 +112,7 @@ def extract_compound_codes(text): else: return np.nan if not enzyme_string_file: - response = requests.get(f'https://rest.kegg.jp/get/{"+".join(ec_list)}') + response = request_with_retry(requests.get, f'https://rest.kegg.jp/get/{"+".join(ec_list)}') if response.status_code == 200: response_string = response.text else: @@ -118,7 +176,7 @@ def extract_secondary_id(identifier, database_list, current_db = ""): def get_kegg_reactions(chunk, reactions_string_file = None): kegg_reaction_dictionary = {} if not reactions_string_file: - response = requests.get(f'https://rest.kegg.jp/get/{ "+".join(chunk)}') + response = request_with_retry(requests.get, f'https://rest.kegg.jp/get/{ "+".join(chunk)}') if response.status_code == 200: responses_string = response.text response_status = 200 @@ -168,9 +226,21 @@ def pubchem_cid_to_descriptor(compound_list, chunk_size = 200): for i in range(0, len(compound_list), chunk_size): chunk = compound_list[i:i+chunk_size] compound_string = ",".join(chunk) - url = f"https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/cid/{compound_string}/property/CanonicalSMILES/JSON" + # PubChem renamed/deprecated the "CanonicalSMILES" property - a live + # request confirmed requesting it now returns the value under a + # "ConnectivitySMILES" JSON key instead (same non-isomeric canonical + # SMILES concept, new name as part of PubChem's PUG-REST schema + # modernization). Request the current name directly rather than + # relying on whatever alias PubChem happens to still honour. + # Also request Title (PubChem's display name) in the same batched + # call - free (one more property on an already-batched request), + # and replaces the separate per-KEGG-code live name lookup that + # used to sit later in the PubChem block (removed: it was the + # single most expensive call in the whole build for a purely + # cosmetic field - see git history for that change). + url = f"https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/cid/{compound_string}/property/ConnectivitySMILES,Title/JSON" # Fetch JSON data from the URL - response = requests.get(url) + response = request_with_retry(requests.get, url) # Check if the request was successful if response.status_code == 200: @@ -192,7 +262,7 @@ def get_kegg_compound_record(kegg_id, compound_cache_dir = None): with open(f"{compound_cache_dir}/{kegg_id}.kegg_record.txt", "r") as file: compound_record_text = file.read() else: - compound_record = requests.get(f'https://rest.kegg.jp/get/{kegg_id}') + compound_record = request_with_retry(requests.get, f'https://rest.kegg.jp/get/{kegg_id}') if compound_record.status_code == 200: compound_record_text = compound_record.text if compound_cache_dir: @@ -221,7 +291,7 @@ def get_kegg_compound_smiles(kegg_id, mol_compound_cache_dir = None): molblock = file.read() else: time.sleep(1) - response = requests.get(f'https://rest.kegg.jp/get/{kegg_id}/mol') + response = request_with_retry(requests.get, f'https://rest.kegg.jp/get/{kegg_id}/mol') if response.status_code == 200: compound_split = response.text.split("> \n") molblock = compound_split[0] @@ -241,7 +311,11 @@ def get_gtc_info(gtcids, cache_df_file): cache_df = pd.read_pickle(cache_df_file) else: cache_df = None - for gtcid in gtcids: + total = len(gtcids) + print(f" [{_now()}] [GlyTouCan] starting, {total} ids") + for i, gtcid in enumerate(gtcids, 1): + if i == 1 or i % 50 == 0 or i == total: + print(f" [{_now()}] [GlyTouCan] {i}/{total}") if cache_df is not None and gtcid in cache_df.index: continue else: @@ -249,7 +323,7 @@ def get_gtc_info(gtcids, cache_df_file): data = {'gtcid': gtcid} nested_dict = {} - response = requests.post(url, data=data) + response = request_with_retry(requests.post, url, data=data) if response.status_code == 200: json_result = json.loads(response.text) # Display the response content @@ -364,7 +438,9 @@ def main(): parser.add_argument('--csdb_cache', type=str, default = None, help='Path to csdb_cache.pkl file, cached from previous run') parser.add_argument('--compound_cache_dir', type=str, default = None, help='Path to directory containing KEGG compound records, cached from previous run') parser.add_argument('--chebi_relations', type=str, default = None, help='Path to chebi relations.tsv file for extracting cofactor information') + parser.add_argument('--chebi_relation_types', type=str, default = None, help='Path to ChEBI relation_type.tsv vocabulary file (id/code lookup for relation.tsv.relation_type_id, e.g. resolving "has_role")') parser.add_argument('--gtc_cache', type=str, default = None, help='Path to glytoucan_cache.pkl file, cached from previous run') + parser.add_argument('--cofactor_ligands', type=str, default = None, help='Path to preprocessed cofactor_ligands_df.pkl from preprocess_cofactors.py (optional - see docs/cofactor_coverage_plan.md)') args = parser.parse_args() if args.compound_cache_dir: @@ -394,7 +470,8 @@ def main(): #run this indiviudally, and save the results to a file for each enzyme. Adapt function to take a cache dir and check for matching filename before making api call print("Fetching enzyme records from KEGG API") for i in range(0, len(ec_list), n): - print(f"Processing chunk {i} of {len(ec_list)}") + if i % 100 == 0 or i + n >= len(ec_list): + print(f" [{_now()}] Processing chunk {i} of {len(ec_list)}") chunk = ec_list[i:i + n] enzyme_dict, enzyme_string = get_kegg_enzymes(chunk) enzyme_records.update(enzyme_dict) @@ -423,9 +500,11 @@ def main(): print("Loading reaction records from text file.") kegg_reaction_dictionary, kegg_reaction_string = get_kegg_reactions(reactions, reactions_string_file = args.kegg_reaction_string) else: - print("Fetching reaction records from KEGG API") + print(f"Fetching reaction records from KEGG API ({len(reactions)} reactions)") n=10 #chunk size for i in range(0, len(reactions), n): + if i % 100 == 0 or i + n >= len(reactions): + print(f" [{_now()}] Processing chunk {i} of {len(reactions)}") chunk = reactions[i:i + n] reaction_dictionary, reaction_string = get_kegg_reactions(chunk) kegg_reaction_dictionary.update(reaction_dictionary) @@ -455,7 +534,20 @@ def main(): kegg_reaction_enzyme_df["EC_substrate_codes"] = kegg_reaction_enzyme_df["EC_substrate_codes"].apply(lambda d: d if isinstance(d, list) else []) kegg_reaction_enzyme_df["EC_product_codes"] = kegg_reaction_enzyme_df["EC_product_codes"].apply(lambda d: d if isinstance(d, list) else []) - kegg_reaction_enzyme_df = kegg_reaction_enzyme_df.fillna("").groupby("entry").agg({"entry" : set, "error" : set, "matched_name" : set, "EC_substrate_codes": sum,"EC_product_codes": sum, "reaction_substrate_codes" : sum, "reaction_product_codes": sum, "EC_reactions":set, "reaction_id" : set, "reaction_definition": set, "reaction_equation" : set}) + # not `sum` directly: pandas' groupby.agg(sum) on an object/list + # column calls Python's builtin sum() with its default start=0, + # which raises "unsupported operand type(s) for +: 'int' and + # 'list'" the moment it tries 0 + - true for every + # group regardless of whether the lists are empty, this was never + # going to work as a list-concatenation aggregator under this + # pandas version. concat_lists explicitly starts from [] instead. + def concat_lists(series): + result = [] + for value in series: + result.extend(value) + return result + + kegg_reaction_enzyme_df = kegg_reaction_enzyme_df.fillna("").groupby("entry").agg({"entry" : set, "error" : set, "matched_name" : set, "EC_substrate_codes": concat_lists,"EC_product_codes": concat_lists, "reaction_substrate_codes" : concat_lists, "reaction_product_codes": concat_lists, "EC_reactions":set, "reaction_id" : set, "reaction_definition": set, "reaction_equation" : set}) kegg_reaction_enzyme_df["entities"] = kegg_reaction_enzyme_df["EC_substrate_codes"] + kegg_reaction_enzyme_df["EC_product_codes"] + kegg_reaction_enzyme_df["reaction_substrate_codes"] + kegg_reaction_enzyme_df["reaction_product_codes"] kegg_reaction_enzyme_df["entities"] = kegg_reaction_enzyme_df["entities"].apply(lambda x: ",".join(list(set(x)))) #convert substrate/product codes to a comma separated string instead of list (for easier merging later) @@ -499,6 +591,17 @@ def unpack_sets(row): rhea_reactions["compound_reaction"] = rhea_reactions["compound_reaction"].str.join("|") print("RHEA records loaded from file") + # cofactor-origin cognate ligands (docs/cofactor_coverage_plan.md) - a + # strict addition alongside the reaction-derived sources above, not a + # modification of them. Optional: the pipeline runs unchanged without + # --cofactor_ligands. + cofactor_ligands_cols = ["entry", "compound_id", "compound_name", "ROMol", "ligand_db", "compound_reaction", "ligand_source"] + if args.cofactor_ligands: + cofactor_ligands_df = pd.read_pickle(args.cofactor_ligands) + print(f"Cofactor ligands loaded from file: {len(cofactor_ligands_df)} rows, {cofactor_ligands_df['entry'].nunique()} distinct ECs") + else: + cofactor_ligands_df = pd.DataFrame(columns=cofactor_ligands_cols) + # ChEBI and PubChem are resolved *before* the live per-compound KEGG lookup # below, specifically so that lookup can skip any compound code already # resolvable from bulk/cheap sources. ChEBI is a local bulk file (zero API @@ -552,7 +655,33 @@ def unpack_sets(row): data.append(record) kegg_pubchem_mapping = pd.DataFrame(data) - pubchem_kegg_compound_records = pd.DataFrame([get_kegg_compound_record(code, compound_cache_dir = args.compound_cache_dir) for code in kegg_pubchem_mapping.KEGG.unique()]) + # Restrict to codes actually seen in our reactions *before* the live + # fetch below - this file is NCBI's full bulk KEGG<->PubChem + # cross-reference (26k+ KEGG codes), not scoped to this run's data at + # all. Previously every one of those 26k+ codes got a live, + # unbatched get_kegg_compound_record call regardless of relevance + # (confirmed: only ~8.2k compound codes are ever relevant here) - + # this was the dominant, silent cost of this whole phase. + kegg_pubchem_mapping = kegg_pubchem_mapping.loc[kegg_pubchem_mapping.KEGG.isin(compound_codes)].reset_index(drop = True) + # Not a live get_kegg_compound_record call per code: the only field + # of that record actually used downstream from this specific merge + # is compound_name, purely as a cosmetic display-name fallback + # (cognate_ligands_df already falls further back to ChEBI_NAME, then + # to the bare compound code itself, if this is missing - see the + # fillna chain after the final concat). It never feeds structure/ + # SMILES resolution (that's CID -> ConnectivitySMILES, entirely + # separate) or any matching/dedup logic. Measured live: ~1.4s/code, + # which made this the dominant cost of the whole build (~2.5h for + # ~6.5k codes) for a purely cosmetic field. Skip the live fetch and + # let the existing name fallback chain handle it instead - matches + # get_kegg_compound_record's own "not found" shape so the merge/ + # assert below are unaffected. + pubchem_kegg_compound_records = pd.DataFrame({ + "compound_id": kegg_pubchem_mapping.KEGG.unique(), + "compound_name": None, + "dbxrefs": None, + "KEGG_compound_record": None, + }) kegg_pubchem_mapping = kegg_pubchem_mapping.merge(pubchem_kegg_compound_records, left_on = "KEGG", right_on = "compound_id", how = "left", indicator = True) assert(len(kegg_pubchem_mapping.loc[kegg_pubchem_mapping._merge != "both"]) == 0) kegg_pubchem_mapping.drop(columns = ["_merge"], inplace = True) @@ -565,12 +694,16 @@ def unpack_sets(row): pubchem_smiles["CID"] = pubchem_smiles["CID"].astype(int) kegg_pubchem_mapping_smiles = kegg_pubchem_mapping.merge(pubchem_smiles, on = "CID", how = "left", indicator = True) - kegg_pubchem_mapping_smiles_filtered = kegg_pubchem_mapping_smiles.loc[kegg_pubchem_mapping_smiles.CanonicalSMILES.isna() == False].reset_index(drop = True).copy() + # compound_name is always None at this point (see above) - fill it + # from PubChem's own Title property, fetched in the same batched + # call as the SMILES above. + kegg_pubchem_mapping_smiles["compound_name"] = kegg_pubchem_mapping_smiles["compound_name"].fillna(kegg_pubchem_mapping_smiles["Title"]) + kegg_pubchem_mapping_smiles_filtered = kegg_pubchem_mapping_smiles.loc[kegg_pubchem_mapping_smiles.ConnectivitySMILES.isna() == False].reset_index(drop = True).copy() kegg_pubchem_mapping_smiles_filtered.drop(columns = "_merge", inplace = True) ### Merging PubChem with enzyme df kegg_reaction_enzyme_df_exploded_pubchem = kegg_reaction_enzyme_df_exploded.merge(kegg_pubchem_mapping_smiles_filtered, left_on = "entities", right_on = "KEGG", how = "inner") - PandasTools.AddMoleculeColumnToFrame(kegg_reaction_enzyme_df_exploded_pubchem, smilesCol='CanonicalSMILES') + PandasTools.AddMoleculeColumnToFrame(kegg_reaction_enzyme_df_exploded_pubchem, smilesCol='ConnectivitySMILES') kegg_reaction_enzyme_df_exploded_pubchem = kegg_reaction_enzyme_df_exploded_pubchem.loc[kegg_reaction_enzyme_df_exploded_pubchem.ROMol.isna() == False].reset_index(drop = True) kegg_reaction_enzyme_df_exploded_pubchem["ligand_db"] = "Pubchem:" + kegg_reaction_enzyme_df_exploded_pubchem["CID"].astype("str") @@ -592,8 +725,16 @@ def unpack_sets(row): already_covered_codes = chebi_covered_codes | pubchem_covered_codes compound_codes_for_kegg = [code for code in compound_codes if code.startswith("C") and code not in already_covered_codes] print(f"Getting KEGG Compound records for {len(compound_codes_for_kegg)} compound codes not already resolved via ChEBI/PubChem (of {len([c for c in compound_codes if c.startswith('C')])} total KEGG compound codes seen)") - kegg_compounds_df = pd.DataFrame([get_kegg_compound_record(code, compound_cache_dir = args.compound_cache_dir) for code in compound_codes_for_kegg]) - kegg_compounds_df["canonical_smiles"] = kegg_compounds_df["compound_id"].apply(lambda x: get_kegg_compound_smiles(x, mol_compound_cache_dir = args.compound_cache_dir)) + kegg_compounds_df = pd.DataFrame(fetch_all_with_progress( + compound_codes_for_kegg, + lambda code: get_kegg_compound_record(code, compound_cache_dir = args.compound_cache_dir), + "KEGG-direct fallback: fetching compound records", + )) + kegg_compounds_df["canonical_smiles"] = fetch_all_with_progress( + kegg_compounds_df["compound_id"].tolist(), + lambda code: get_kegg_compound_smiles(code, mol_compound_cache_dir = args.compound_cache_dir), + "KEGG-direct fallback: fetching compound SMILES", + ) kegg_compounds_df.to_pickle(f"kegg_compounds_df.pkl") print("KEGG Compound records saved") else: @@ -620,12 +761,12 @@ def unpack_sets(row): & (kegg_reaction_enzyme_df_exploded.entities.str.startswith("G")), "entities"].values.tolist() ### Getting KEGG compound records for glycans (to get xref to glytoucan) - glycan_compounds = [] + glycan_compounds = fetch_all_with_progress( + glycans, + lambda glycan: get_kegg_compound_record(glycan, compound_cache_dir=args.compound_cache_dir), + "Glycans: fetching KEGG compound records", + ) - for glycan in glycans: - compound = get_kegg_compound_record(glycan, compound_cache_dir=args.compound_cache_dir) - glycan_compounds.append(compound) - glycan_compounds_df = pd.DataFrame(glycan_compounds) glycan_compounds_df[["source", "secondary_id"]] = glycan_compounds_df.apply(lambda x: extract_secondary_id(x["compound_id"], x["dbxrefs"]), axis = 1, result_type = "expand").values glytoucan_ids = glycan_compounds_df.loc[(glycan_compounds_df.secondary_id.isna() == False) & @@ -654,6 +795,13 @@ def unpack_sets(row): #all - benchmarked at ~89% on the real cognate-ligand master set. Kept as a #fallback rather than replacing the live chain outright, in case GlyTouCan's #API is fixed in future. + # object dtype, not whatever pandas inferred from line 784's mostly- + # NaN result (often float64): a .loc[mask, col] = + # partial assignment into a float64-typed column raises under this + # pandas version ("Invalid value '...' for dtype + # 'float64'") - same class of bug already fixed in preprocess_rhea.py + # today (assigning str/NaN into a strictly-typed column). + glycan_compounds_df_merged["smiles"] = glycan_compounds_df_merged["smiles"].astype(object) missing_smiles_mask = glycan_compounds_df_merged["smiles"].isna() & glycan_compounds_df_merged["wurcs"].notna() glycan_compounds_df_merged.loc[missing_smiles_mask, "smiles"] = glycan_compounds_df_merged.loc[missing_smiles_mask, "wurcs"].apply(get_smiles_from_wurcs_offline) @@ -672,23 +820,25 @@ def unpack_sets(row): kegg_reaction_enzyme_df_exploded_gtc = pd.read_pickle(f"kegg_reaction_enzyme_df_exploded_gtc.pkl") if not os.path.exists(f"cognate_ligands_df.pkl"): - cognate_ligands_df = pd.concat([rhea_reactions[["entry", "compound_id", "compound_name", "ROMol", "ligand_db", "compound_reaction"]], - kegg_reaction_enzyme_df_exploded_kegg[["entry", "compound_name", "compound_id", "ROMol", "ligand_db", "compound_reaction"]], - kegg_reaction_enzyme_df_exploded_chebi[["entry", "ChEBI_NAME", "KEGG COMPOUND ACCESSION", "ROMol", "ligand_db", "compound_reaction"]].rename(columns = {"KEGG COMPOUND ACCESSION" : "compound_id"}), - kegg_reaction_enzyme_df_exploded_pubchem[["entry", "compound_name", "KEGG", "ROMol", "ligand_db", "compound_reaction"]].rename(columns = {"KEGG" : "compound_id"}), - kegg_reaction_enzyme_df_exploded_gtc[["entry", "compound_name", "compound_id","ROMol", "ligand_db", "compound_reaction"]]]) + cognate_ligands_df = pd.concat([rhea_reactions[["entry", "compound_id", "compound_name", "ROMol", "ligand_db", "compound_reaction"]].assign(ligand_source = "reaction"), + kegg_reaction_enzyme_df_exploded_kegg[["entry", "compound_name", "compound_id", "ROMol", "ligand_db", "compound_reaction"]].assign(ligand_source = "reaction"), + kegg_reaction_enzyme_df_exploded_chebi[["entry", "ChEBI_NAME", "KEGG COMPOUND ACCESSION", "ROMol", "ligand_db", "compound_reaction"]].rename(columns = {"KEGG COMPOUND ACCESSION" : "compound_id"}).assign(ligand_source = "reaction"), + kegg_reaction_enzyme_df_exploded_pubchem[["entry", "compound_name", "KEGG", "ROMol", "ligand_db", "compound_reaction"]].rename(columns = {"KEGG" : "compound_id"}).assign(ligand_source = "reaction"), + kegg_reaction_enzyme_df_exploded_gtc[["entry", "compound_name", "compound_id","ROMol", "ligand_db", "compound_reaction"]].assign(ligand_source = "reaction"), + cofactor_ligands_df[cofactor_ligands_cols]]) cognate_ligands_df = cognate_ligands_df.reset_index() - + #fill the missing compound names first using the chebi name, and subsequently with the compound id if that is also nan. cognate_ligands_df["compound_name"] = cognate_ligands_df["compound_name"].fillna(cognate_ligands_df["ChEBI_NAME"]).fillna(cognate_ligands_df["compound_id"]) cognate_ligands_df["ROMol"] = cognate_ligands_df["ROMol"].apply(lambda x: neutralize_atoms(x) if isinstance(x,Chem.rdchem.Mol) else np.nan) #attempt to neutralise charged structures for grouping as charges cannot be used to score mols cognate_ligands_df["canonical_smiles"] = cognate_ligands_df["ROMol"].map(lambda x: canon_smiles(x) if isinstance(x,Chem.rdchem.Mol) else np.nan) - cognate_ligands_df_unique_smiles = cognate_ligands_df[["canonical_smiles", "compound_name", "ligand_db", "compound_reaction"]].copy() + cognate_ligands_df_unique_smiles = cognate_ligands_df[["canonical_smiles", "compound_name", "ligand_db", "compound_reaction", "ligand_source"]].copy() cognate_ligands_df_unique_smiles["compound_reaction"] = cognate_ligands_df_unique_smiles["compound_reaction"].fillna("") - cognate_ligands_df_unique_smiles = cognate_ligands_df_unique_smiles.groupby("canonical_smiles", dropna = False).agg({"compound_name": set, "ligand_db": set, "compound_reaction": set}).reset_index() + cognate_ligands_df_unique_smiles = cognate_ligands_df_unique_smiles.groupby("canonical_smiles", dropna = False).agg({"compound_name": set, "ligand_db": set, "compound_reaction": set, "ligand_source": set}).reset_index() cognate_ligands_df_unique_smiles["ligand_db"] = cognate_ligands_df_unique_smiles.ligand_db.str.join("|") cognate_ligands_df_unique_smiles["compound_name"] = cognate_ligands_df_unique_smiles.compound_name.str.join("|") cognate_ligands_df_unique_smiles["compound_reaction"] = cognate_ligands_df_unique_smiles.compound_reaction.str.join("|").str.strip("|") + cognate_ligands_df_unique_smiles["ligand_source"] = cognate_ligands_df_unique_smiles.ligand_source.apply(lambda x: "|".join(sorted(x))) cognate_ligands_df_unique_smiles = cognate_ligands_df_unique_smiles.reset_index(drop=True).reset_index() cognate_ligands_df_unique_smiles.rename(columns = {"index": "uniqueID"}, inplace = True) @@ -696,9 +846,33 @@ def unpack_sets(row): cognate_ligands_df["compound_name"] = cognate_ligands_df["compound_name"].str.split("|").apply(lambda x: length_upper_sorted(x)).str.join("|") cognate_ligands_df = cognate_ligands_df.drop_duplicates() + # ChEBI's relation.tsv moved from a string TYPE column to a numeric + # relation_type_id foreign key (and INIT_ID/FINAL_ID -> lowercase + # init_id/final_id) - resolve "has_role" via the small + # relation_type.tsv vocabulary file rather than hardcoding whatever + # numeric id it happens to be today (confirmed live: 4, but that's + # an internal ChEBI id, not a stable public constant). + # + # Direction, confirmed against real data: init_id is the specific + # compound, final_id is the role class it has + # (glucose --has_role--> nutrient) - the ontologically-expected + # direction. Filtering init_id against the role IDs (as an earlier + # version of this code did) matches zero rows on current data; + # role IDs belong on final_id. Renamed below (final_id -> INIT_ID, + # init_id -> FINAL_ID) so the rest of this block, which already + # expects INIT_ID to be the role-class side and FINAL_ID to be the + # specific-compound side it later merges cognate_ligands_df against, + # is otherwise unchanged. chebi_relations = pd.read_csv(f"{args.chebi_relations}", sep="\t") - - chebi_cofactors = chebi_relations.loc[(chebi_relations.TYPE == "has_role") & (chebi_relations.INIT_ID.isin([23357,23354,26348,26672])), ["TYPE", "INIT_ID", "FINAL_ID"]] + chebi_relation_types = pd.read_csv(f"{args.chebi_relation_types}", sep="\t") + has_role_id = chebi_relation_types.loc[chebi_relation_types.code == "has_role", "id"].iloc[0] + + chebi_cofactors = chebi_relations.loc[(chebi_relations.relation_type_id == has_role_id) & (chebi_relations.final_id.isin([23357,23354,26348,26672])), ["relation_type_id", "final_id", "init_id"]].copy() + chebi_cofactors.rename(columns = {"relation_type_id": "TYPE", "final_id": "INIT_ID", "init_id": "FINAL_ID"}, inplace = True) + # TYPE is int64 at this point (it's relation_type_id, always + # has_role_id here) - the string labels assigned below need object + # dtype first, same strict-coercion issue as elsewhere in this file. + chebi_cofactors["TYPE"] = chebi_cofactors["TYPE"].astype(object) chebi_cofactors.loc[chebi_cofactors.INIT_ID == 23357, "TYPE"] = "Cofactor" chebi_cofactors.loc[chebi_cofactors.INIT_ID == 23354, "TYPE"] = "Coenzyme" chebi_cofactors.loc[chebi_cofactors.INIT_ID == 26348, "TYPE"] = "Prosthetic Group" diff --git a/nextflow/bin/preprocess_cofactors.py b/nextflow/bin/preprocess_cofactors.py new file mode 100644 index 0000000..c980e83 --- /dev/null +++ b/nextflow/bin/preprocess_cofactors.py @@ -0,0 +1,261 @@ +#!/usr/bin/env python + +""" +Preprocesses cofactor-EC association data into a cognate-ligand-shaped +dataframe, for the coverage-gap plan in docs/cofactor_coverage_plan.md. +Mirrors preprocess_rhea.py's pattern: a standalone script producing a +pickle (cofactor_ligands_df.pkl) that get_ec_information.py reads in and +concatenates alongside its existing reaction-derived sources, tagged +ligand_source="cofactor". + +Two combined, complementary sources (see the plan doc for the coverage +analysis behind this choice): + + 1. CoFactor DB 2010 + BRENDA, via RelLig's own vendored cofactor_ec.csv / + cofactors_details.json (Apache-2.0, pdberellig project) - 27 organic + cofactor classes, each EC exact (4-digit), resolved to a SMILES via + its representative PDB CCD code. + 2. UniProt's `COFACTOR` annotation, bulk-pulled (structured + Name=...;Xref=ChEBI:... rows only - free-text-only rows, e.g. + chlorophyll's photosystem entries, are out of scope here and are a + deferred separate LLM-extraction follow-on per the plan doc). + +Both sources' EC values get classified by how many of their 4 segments are +fully resolved (see utils.classify_ec_completeness): exact (4) passes +through, subsubclass-level partials ("N.N.N.-", 3) get broadcast to every +matching terminal EC via utils.broadcast_subsubclass_ec, and anything +coarser (class/subclass-level, 1 or 2) is dropped and logged rather than +broadcast - broadcasting those was checked and found to produce nonsense +(e.g. a class-level wildcard matching literally every terminal EC in that +class). See docs/cofactor_coverage_plan.md ("The broadcast problem, and +its fix") for the full reasoning and the numbers behind this rule. + +Example usage: + python3 preprocess_cofactors.py \ + --cofactor_ec_csv cofactor_ec.csv \ + --cofactor_details_json cofactors_details.json \ + --ccd_cif ccd.cif \ + --uniprot_cofactor_tsv uniprot_cofactor_annotations.tsv \ + --chebi_structures chebi_structures.tsv.gz \ + --chebi_names chebi_names.tsv.gz \ + --ec_dat enzyme.dat \ + --enzyme_class_file enzclass.txt \ + --outdir /path/to/output/directory +""" + +import argparse +import json +import re +from pathlib import Path + +import pandas as pd +from gemmi import cif +from rdkit.Chem import PandasTools + +from utils import ( + process_ec_records, + classify_ec_completeness, + broadcast_subsubclass_ec, + get_chem_comp_descriptors, +) + +CHEBI_XREF_PATTERN = re.compile(r"Xref=ChEBI:CHEBI:(\d+)") + + +def load_cofactor_db_table(cofactor_ec_csv_path, cofactor_details_json_path): + """Returns a DataFrame with one row per (EC, cofactor class) pair from + CoFactor DB 2010 + BRENDA, columns: ec, cofactor_id, representative_ccd, + source. cofactor_ec.csv is EC_NO/COFACTOR_ID/SOURCE; cofactors_details.json + is a list of {id, representative, template, threshold} - representative + is the PDB CCD code used to resolve a SMILES for that cofactor class.""" + cofactor_ec = pd.read_csv(cofactor_ec_csv_path) + cofactor_ec.columns = [c.strip().strip('"') for c in cofactor_ec.columns] + cofactor_ec["EC_NO"] = cofactor_ec["EC_NO"].astype(str).str.strip().str.strip('"') + cofactor_ec["COFACTOR_ID"] = cofactor_ec["COFACTOR_ID"].astype(int) + + with open(cofactor_details_json_path) as handle: + details = json.load(handle) + details_df = pd.DataFrame(details)[["id", "representative"]].rename( + columns={"id": "COFACTOR_ID", "representative": "representative_ccd"} + ) + + merged = cofactor_ec.merge(details_df, on="COFACTOR_ID", how="inner") + merged = merged.rename(columns={"EC_NO": "ec", "COFACTOR_ID": "cofactor_id", "SOURCE": "source"}) + return merged[["ec", "cofactor_id", "representative_ccd", "source"]].drop_duplicates() + + +def parse_uniprot_cofactor_df(raw_df): + """Takes the bulk UniProt TSV (columns: Entry, EC number, Cofactor) as + already-loaded a DataFrame, and returns one row per (ec, chebi_id) pair + for rows with at least one structured `Xref=ChEBI:CHEBI:` in the + Cofactor field. Rows with only a free-text Note (no Xref) are dropped + here - confirmed in the plan doc these are ~0.7% of the bulk pull + (chlorophyll's photosystem annotations among them) and are out of scope + for this mechanical path.""" + records = [] + for _, row in raw_df.iterrows(): + cofactor_text = row.get("Cofactor") + if not isinstance(cofactor_text, str): + continue + chebi_ids = CHEBI_XREF_PATTERN.findall(cofactor_text) + if not chebi_ids: + continue + ec_field = row.get("EC number") + if not isinstance(ec_field, str) or not ec_field.strip(): + continue + ecs = [e.strip() for e in ec_field.split(";") if e.strip()] + for ec in ecs: + for chebi_id in chebi_ids: + records.append({"ec": ec, "chebi_id": int(chebi_id)}) + if not records: + return pd.DataFrame(columns=["ec", "chebi_id"]) + return pd.DataFrame(records).drop_duplicates() + + +def classify_and_split_ec_rows(df, terminal_ec_list, ec_col="ec"): + """Splits a dataframe by how resolved its EC values are + (utils.classify_ec_completeness): exact (4) rows pass through + unchanged; subsubclass-level (3) rows are exploded into one row per + matching terminal EC (utils.broadcast_subsubclass_ec); class/subclass- + level (1 or 2) rows are dropped, not broadcast. Returns (kept_df, + dropped_df) - dropped_df is for logging/auditability, not silent + discard (docs/cofactor_coverage_plan.md's stated requirement).""" + terminal_ec_list = list(terminal_ec_list) + kept_rows = [] + dropped_rows = [] + + for _, row in df.iterrows(): + ec = row[ec_col] + level = classify_ec_completeness(ec) + if level == 4: + kept_rows.append({**row.to_dict(), "entry": ec}) + elif level == 3: + matches = broadcast_subsubclass_ec(ec, terminal_ec_list) + for terminal_ec in matches: + kept_rows.append({**row.to_dict(), "entry": terminal_ec}) + if not matches: + dropped_rows.append({**row.to_dict(), "reason": "subsubclass-level, no matching terminal EC"}) + else: + dropped_rows.append({**row.to_dict(), "reason": f"too coarse to broadcast safely (level {level})"}) + + kept_df = pd.DataFrame(kept_rows).drop(columns=[ec_col], errors="ignore") if kept_rows else pd.DataFrame(columns=list(df.columns) + ["entry"]).drop(columns=[ec_col], errors="ignore") + dropped_df = pd.DataFrame(dropped_rows) if dropped_rows else pd.DataFrame(columns=list(df.columns) + ["reason"]) + return kept_df, dropped_df + + +def resolve_chebi_smiles_map(chebi_ids, chebi_structures_df, chebi_names_df): + """Resolves a list of bare integer ChEBI IDs to SMILES/name, from the + full ChEBI bulk structures/names files (not the pipeline's existing + ChEBI_Results.tsv, which is deliberately narrowed to KEGG-COMPOUND- + cross-referenced entries only - most cofactor ChEBI IDs, especially + metal ions, have no KEGG COMPOUND cross-reference at all and would be + silently dropped if resolved that way).""" + ids = pd.Series(sorted(set(chebi_ids)), name="compound_id") + structures = chebi_structures_df.dropna(subset=["smiles"]).drop_duplicates(subset="compound_id", keep="first") + names = chebi_names_df.groupby("compound_id").agg({"name": "first"}).reset_index() + + resolved = ids.to_frame().merge(structures[["compound_id", "smiles"]], on="compound_id", how="inner") + resolved = resolved.merge(names, on="compound_id", how="left") + return resolved + + +def build_cofactor_ligands_df(cofactor_db_df, uniprot_df, terminal_ec_list, ccd_doc, chebi_structures_df, chebi_names_df): + """Orchestrates both sources into a single dataframe shaped like + get_ec_information.py's other cognate-ligand source frames (entry, + compound_id, compound_name, ROMol, ligand_db, compound_reaction, + ligand_source), ready to concat directly into cognate_ligands_df. + Returns (cofactor_ligands_df, dropped_ec_log_df).""" + cofactor_db_kept, cofactor_db_dropped = classify_and_split_ec_rows(cofactor_db_df, terminal_ec_list) + uniprot_kept, uniprot_dropped = classify_and_split_ec_rows(uniprot_df, terminal_ec_list) + + rows = [] + + if not cofactor_db_kept.empty: + ccd_codes = cofactor_db_kept["representative_ccd"].dropna().unique().tolist() + ccd_smiles = get_chem_comp_descriptors(ccd_doc, ccd_codes) + for _, row in cofactor_db_kept.iterrows(): + smiles = ccd_smiles.get(row["representative_ccd"]) + if smiles is None: + continue + rows.append({ + "entry": row["entry"], + "compound_id": row["representative_ccd"], + "compound_name": row["representative_ccd"], + "smiles": smiles, + "ligand_db": f"CofactorDB:{row['cofactor_id']}", + }) + + if not uniprot_kept.empty: + chebi_map = resolve_chebi_smiles_map(uniprot_kept["chebi_id"].unique(), chebi_structures_df, chebi_names_df) + chebi_map = chebi_map.set_index("compound_id") + for _, row in uniprot_kept.iterrows(): + chebi_id = row["chebi_id"] + if chebi_id not in chebi_map.index: + continue + resolved = chebi_map.loc[chebi_id] + rows.append({ + "entry": row["entry"], + "compound_id": f"CHEBI:{chebi_id}", + "compound_name": resolved["name"] if pd.notna(resolved.get("name")) else f"CHEBI:{chebi_id}", + "smiles": resolved["smiles"], + "ligand_db": f"CHEBI:{chebi_id}", + }) + + dropped_df = pd.concat([cofactor_db_dropped, uniprot_dropped], ignore_index=True) + + if not rows: + empty = pd.DataFrame(columns=["entry", "compound_id", "compound_name", "ROMol", "ligand_db", "compound_reaction", "ligand_source"]) + return empty, dropped_df + + cofactor_ligands_df = pd.DataFrame(rows).drop_duplicates(subset=["entry", "compound_id", "smiles"]) + PandasTools.AddMoleculeColumnToFrame(cofactor_ligands_df, smilesCol="smiles", molCol="ROMol") + cofactor_ligands_df = cofactor_ligands_df.loc[cofactor_ligands_df["ROMol"].notna()].copy() + cofactor_ligands_df["compound_reaction"] = "" + cofactor_ligands_df["ligand_source"] = "cofactor" + cofactor_ligands_df = cofactor_ligands_df[ + ["entry", "compound_id", "compound_name", "ROMol", "ligand_db", "compound_reaction", "ligand_source"] + ].reset_index(drop=True) + + return cofactor_ligands_df, dropped_df + + +def main(): + parser = argparse.ArgumentParser(description="Preprocess cofactor-EC association data") + parser.add_argument("--cofactor_ec_csv", required=True, help="RelLig's vendored cofactor_ec.csv (EC_NO/COFACTOR_ID/SOURCE)") + parser.add_argument("--cofactor_details_json", required=True, help="RelLig's vendored cofactors_details.json (cofactor class -> representative CCD code)") + parser.add_argument("--ccd_cif", required=True, help="cif file containing the chemical component dictionary in mmcif format") + parser.add_argument("--uniprot_cofactor_tsv", required=True, help="Bulk UniProt TSV (Entry, EC number, Cofactor columns)") + parser.add_argument("--chebi_structures", required=True, help="ChEBI bulk structures.tsv.gz file") + parser.add_argument("--chebi_names", required=True, help="ChEBI bulk names.tsv.gz file") + parser.add_argument("--ec_dat", required=True, help="Path to enzyme.dat file from EXPASY") + parser.add_argument("--enzyme_class_file", required=True, help="Path to enzyme_class file") + parser.add_argument("--outdir", required=True, help="Output directory") + args = parser.parse_args() + + Path(args.outdir).mkdir(parents=True, exist_ok=True) + + ec_records_df = process_ec_records(args.ec_dat, args.enzyme_class_file) + terminal_ec_list = ec_records_df.TRANSFER.unique().tolist() + + cofactor_db_df = load_cofactor_db_table(args.cofactor_ec_csv, args.cofactor_details_json) + + uniprot_raw_df = pd.read_csv(args.uniprot_cofactor_tsv, sep="\t") + uniprot_df = parse_uniprot_cofactor_df(uniprot_raw_df) + + ccd_doc = cif.read(args.ccd_cif) + chebi_structures_df = pd.read_csv(args.chebi_structures, sep="\t", compression="gzip") + chebi_names_df = pd.read_csv(args.chebi_names, sep="\t", compression="gzip") + + cofactor_ligands_df, dropped_ec_log = build_cofactor_ligands_df( + cofactor_db_df, uniprot_df, terminal_ec_list, ccd_doc, chebi_structures_df, chebi_names_df + ) + + cofactor_ligands_df.to_pickle(f"{args.outdir}/cofactor_ligands_df.pkl") + dropped_ec_log.to_csv(f"{args.outdir}/cofactor_dropped_ec_log.tsv", sep="\t", index=False) + + print(f"Cofactor ligands: {len(cofactor_ligands_df)} rows, {cofactor_ligands_df['entry'].nunique()} distinct ECs") + print(f"Dropped EC rows (too coarse to broadcast, or unresolved): {len(dropped_ec_log)} - see cofactor_dropped_ec_log.tsv") + + +if __name__ == "__main__": + main() diff --git a/nextflow/bin/preprocess_rhea.py b/nextflow/bin/preprocess_rhea.py index 002e8f4..121fccc 100644 --- a/nextflow/bin/preprocess_rhea.py +++ b/nextflow/bin/preprocess_rhea.py @@ -51,7 +51,16 @@ def main(): for reactant in reactants: reactant_smiles = reactant.smiles - reactant_id = reactant.metadata["molecule_name"] + # not reactant.metadata["molecule_name"]: rdfreader's + # .metadata property fixed-column-parses the whole mol + # block header including the program/timestamp line, + # which Rhea's CDK-written molfiles don't conform to + # (e.g. " CDK 2/12/10,15:27" instead of MDL's fixed + # MMDDYY field) - breaks on ~97% of real Rhea rd/ files. + # The molecule name is just the mol block's raw first + # line; read it directly instead of going through the + # metadata parser at all. + reactant_id = reactant.mol_block.splitlines()[0].strip() mol_dict[unique_id] = {"reaction_id": reaction_id, "reaction_properties": reaction_properties, "reaction_smiles": reaction_smiles, @@ -62,7 +71,7 @@ def main(): for product in products: product_smiles = product.smiles - product_id = product.metadata["molecule_name"] + product_id = product.mol_block.splitlines()[0].strip() mol_dict[unique_id] = {"reaction_id": reaction_id, "reaction_properties": reaction_properties, "reaction_smiles": reaction_smiles, @@ -87,10 +96,19 @@ def main(): reactions_df_merged = reactions_df.merge(rheamerge[["RHEA_ID_LR", "ID"]], left_on = "reaction_id", right_on = "RHEA_ID_LR", how = "inner") reactions_df_merged.loc[reactions_df_merged.compound_id.str.startswith("CHEBI"), "COMPOUND_ID"] = reactions_df_merged.loc[reactions_df_merged.compound_id.str.startswith("CHEBI"), "compound_id"].apply(lambda x: re.findall(r"CHEBI:(\d+)", x)[0]) #in chebi names format - reactions_df_merged.loc[reactions_df_merged.compound_id.str.startswith("CHEBI") == False, "COMPOUND_ID"] = -1 + # "-1" as a str, not int: the line above leaves COMPOUND_ID as a + # pandas StringDtype column, and modern pandas rejects assigning an + # int into a StringDtype column in place. The next line casts the + # whole column to int anyway, so the str/int distinction here doesn't + # matter to the result. + reactions_df_merged.loc[reactions_df_merged.compound_id.str.startswith("CHEBI") == False, "COMPOUND_ID"] = "-1" reactions_df_merged["COMPOUND_ID"] = reactions_df_merged["COMPOUND_ID"].astype("int") chebi_names = pd.read_csv(f"{args.chebi_names}", sep = "\t", compression = "gzip") + # ChEBI's flat files moved to lowercase columns upstream (see + # download_reference_data.py's derive_chebi_results docstring) - rename + # to the uppercase names the rest of this script already uses. + chebi_names.rename(columns = {"compound_id": "COMPOUND_ID", "name": "NAME"}, inplace = True) #get the first name for each compound ID in the chebi names file chebi_names = chebi_names.groupby("COMPOUND_ID").agg({"NAME": "first"}).reset_index() print(chebi_names) diff --git a/nextflow/bin/process_all_pdb_contacts.py b/nextflow/bin/process_all_pdb_contacts.py index 39cd605..d37a4b6 100644 --- a/nextflow/bin/process_all_pdb_contacts.py +++ b/nextflow/bin/process_all_pdb_contacts.py @@ -3,7 +3,7 @@ import argparse import pandas as pd from gemmi import cif -from utils import process_ec_records, get_updated_enzyme_records, get_scop_domains_info, extract_interpro_domain_annotations, get_pfam_annotations, get_glycoct_from_wurcs, get_csdb_from_glycoct, get_smiles_from_csdb, get_smiles_from_wurcs_offline, build_cath_dataframe, parse_cddf, build_g3dsa_dataframe, get_scop2_domains_info +from utils import process_ec_records, get_updated_enzyme_records, get_scop_domains_info, extract_interpro_domain_annotations, get_pfam_annotations, get_glycoct_from_wurcs, get_csdb_from_glycoct, get_smiles_from_csdb, get_smiles_from_wurcs_offline, build_cath_dataframe, parse_cddf, build_g3dsa_dataframe, get_scop2_domains_info, get_chem_comp_descriptors import numpy as np from Bio.ExPASy import Enzyme as EEnzyme import re @@ -36,32 +36,6 @@ def get_sugar_smiles_from_wurcs(wurcs_list, csdb_linear_cache, smiles_cache, gly updated_smiles_cache_df = pd.concat([pd.DataFrame(updated_smiles_cache, columns = ["csdb", "descriptor"]), smiles_cache]).drop_duplicates() return sugar_smiles, updated_glycoct_cache_df, updated_csdb_cache_df, updated_smiles_cache_df -def get_chem_comp_descriptors(ccd_doc, comp_id_list): - ligand_descriptors = {} - for ligand in comp_id_list: - lig_descriptor = None - lig_block = ccd_doc.find_block(ligand) - if lig_block is not None: - lig_descriptors = pd.DataFrame(lig_block.find_mmcif_category("_pdbx_chem_comp_descriptor."), columns = ["comp_id", "type", "program", "program_version", "descriptor"]) - lig_descriptors["descriptor"] = lig_descriptors.descriptor.str.strip("\"|';").str.replace(r"\n$","", regex = True) - lig_descriptors = lig_descriptors.loc[lig_descriptors.type == "SMILES"] - PandasTools.AddMoleculeColumnToFrame(lig_descriptors, smilesCol='descriptor', molCol='pdb_ROMol') - lig_descriptors = lig_descriptors.loc[lig_descriptors.pdb_ROMol.isna() == False] - if len(lig_descriptors) == 0: - lig_descriptor = None - else: - #preference is to use openeye descriptors where available. if not, revert to the first smiles string able to be loaded into RDkit. - preferred_row = lig_descriptors.loc[lig_descriptors.program.str.startswith("OpenEye")] - if not preferred_row.empty: - lig_descriptor = preferred_row.iloc[0].descriptor - else: - # Otherwise, select the first row with a SMILES string - lig_descriptor = lig_descriptors.iloc[0].descriptor - ligand_descriptors[ligand] = lig_descriptor - else: - ligand_descriptors[ligand] = None - return ligand_descriptors - def process_sifts_ec_map(sifts_ec_mapping_file, ec_records_file): sifts_chains_ec = sifts_ec_mapping_file.loc[sifts_ec_mapping_file.EC_NUMBER != "?"].copy() #the sifts mapping often has quotes in ec numbers because of extraction with gemmi it seems, strip these e.g. 2ex1 diff --git a/nextflow/bin/process_pdb_contacts.py b/nextflow/bin/process_pdb_contacts.py index f5bc555..d93ffb6 100644 --- a/nextflow/bin/process_pdb_contacts.py +++ b/nextflow/bin/process_pdb_contacts.py @@ -70,7 +70,13 @@ def assign_ownership_percentile_categories(ligands_df, unique_id = "uniqueID", d (ligands_df["domain_contact_perc"] > 0.1) & (ligands_df["domain_contact_perc"] < 0.5) & (ligands_df["num_non_minor_domains"] > 1), "partner", np.where( - ligands_df["domain_contact_perc"] <= 0.1, "minor", np.nan) + # None, not np.nan: numpy 2.x's stricter dtype + # promotion rejects mixing a string array with a + # float nan in np.where ("could not be promoted"), + # where it used to silently fall back to object + # dtype. None is still treated as missing by pandas + # once this array is used downstream. + ligands_df["domain_contact_perc"] <= 0.1, "minor", None) ) ) ) @@ -287,6 +293,11 @@ def process_manifest_row(row, cutoff): 'xref_db': db_source_list, 'xref_db_version': db_version_list }) + # object dtype: xref_db starts as a plain string column (pandas + # StringDtype under this pandas version), and the line below + # assigns a list into some of its cells - same strict-dtype- + # coercion issue fixed elsewhere in this codebase today. + db_df["xref_db"] = db_df["xref_db"].astype(object) db_df.loc[db_df.xref_db == "SCOP2", "xref_db"] = db_df.loc[db_df.xref_db == "SCOP2", "xref_db"].apply(lambda x: ["SCOP2_SuperFamily", "SCOP2_Family"]) #we adapt the SCOP2 source to the db formatting we use in ProCogGraph db_df["xref_db"] = db_df["xref_db"].apply(lambda x: [x] if isinstance(x, str) else x) db_df = db_df.explode("xref_db") diff --git a/nextflow/bin/reference_data_manifest.yaml b/nextflow/bin/reference_data_manifest.yaml index 935f365..d387ec8 100644 --- a/nextflow/bin/reference_data_manifest.yaml +++ b/nextflow/bin/reference_data_manifest.yaml @@ -286,7 +286,26 @@ entries: note: > Upstream file is now gzipped (relation.tsv.gz); get_ec_information.py --chebi_relations expects a plain .tsv, so this is gunzipped on - download. + download. Also note the TYPE/INIT_ID/FINAL_ID columns documented in + older versions of this project are gone - the real columns are + lowercase (init_id/final_id), and TYPE is now a numeric + relation_type_id foreign key, not a string label - see + chebi_relation_types below for the vocabulary needed to resolve it + (e.g. "has_role"). + + - name: chebi_relation_types + target_filename: relation_type.tsv + param_name: null + source_type: direct_url + url: https://ftp.ebi.ac.uk/pub/databases/chebi/flat_files/relation_type.tsv.gz + min_size_bytes: 100 + post_process: gunzip + confidence: verified + note: > + Small id/code vocabulary (11 rows) for relation.tsv's + relation_type_id column - get_ec_information.py + --chebi_relation_types uses it to resolve "has_role" to its current + numeric id (confirmed live: 4) rather than hardcoding that number. - name: chebi_database_accession target_filename: chebi_database_accession.tsv.gz @@ -326,29 +345,63 @@ entries: those three entries to have been fetched first; the script fetches them automatically as dependencies if missing. - # --- Not required by the pipeline today - only needed once - # docs/cofactor_coverage_plan.md is implemented. Kept here so this - # manifest stays the single place tracking every external file the - # pipeline depends on, present or planned. Excluded from --dry-run's - # "ready to run pipeline" summary; fetch explicitly with --only. - - name: pdbe_rellig_chain_ligand_functions - target_filename: interacting_chains_with_ligand_functions.tsv + # --- Optional - only needed if preprocess_cofactors.py is run (see + # docs/cofactor_coverage_plan.md). Not required by the core pipeline; + # excluded from --dry-run's default summary, fetch explicitly with + # --only or --include-optional. NOTE: an earlier version of this section + # pointed at PDBe RelLig's bulk PDBeChem v2 TSVs + # (interacting_chains_with_ligand_functions.tsv / pdb_bound_molecules.tsv) + # as the primary source - checked against the real files and that was + # wrong (no EC column in either; RelLig's actual EC-bearing output is + # per-ligand JSON requiring the pipeline to be run, not a static bulk + # download). Replaced with RelLig's own small vendored EC-cofactor + # mapping plus a bulk UniProt COFACTOR pull instead - see the plan doc's + # "Sources considered" section for the full reasoning. + - name: cofactor_ec_csv + target_filename: cofactor_ec.csv param_name: null source_type: direct_url - url: https://ftp.ebi.ac.uk/pub/databases/msd/pdbechem_v2/additional_data/pdb_ligand_interactions/interacting_chains_with_ligand_functions.tsv - min_size_bytes: 100000000 + url: https://raw.githubusercontent.com/PDBeurope/rellig/main/pdberellig/data/cofactors/cofactor_ec.csv + min_size_bytes: 50000 post_process: none confidence: verified optional: true - note: For the cofactor-coverage plan (docs/cofactor_coverage_plan.md), not the current pipeline. + note: > + RelLig's (Apache-2.0, redistributable) EC_NO/COFACTOR_ID/SOURCE + table, sourced from CoFactor DB 2010 + BRENDA - 27 cofactor classes, + 2,760 distinct ECs, all exact (4-digit). Consumed by + preprocess_cofactors.py --cofactor_ec_csv. - - name: pdbe_rellig_bound_molecules - target_filename: pdb_bound_molecules.tsv + - name: cofactor_details_json + target_filename: cofactors_details.json param_name: null source_type: direct_url - url: https://ftp.ebi.ac.uk/pub/databases/msd/pdbechem_v2/additional_data/pdb_ligand_interactions/pdb_bound_molecules.tsv - min_size_bytes: 50000000 + url: https://raw.githubusercontent.com/PDBeurope/rellig/main/pdberellig/data/cofactors/cofactors_details.json + min_size_bytes: 500 post_process: none confidence: verified optional: true - note: For the cofactor-coverage plan (docs/cofactor_coverage_plan.md), not the current pipeline. + note: > + Maps each of the 27 cofactor_ec_csv COFACTOR_IDs to a representative + PDB CCD code (e.g. id 4 -> NAD), resolved to a SMILES via the + already-fetched ccd.cif. Consumed by preprocess_cofactors.py + --cofactor_details_json. + + - name: uniprot_cofactor_annotations + target_filename: uniprot_cofactor_annotations.tsv + param_name: null + source_type: direct_url + url: "https://rest.uniprot.org/uniprotkb/stream?query=reviewed:true+AND+cc_cofactor:*+AND+ec:*&fields=accession,ec,cc_cofactor&format=tsv" + min_size_bytes: 5000000 + post_process: none + confidence: verified + optional: true + note: > + Bulk (not per-accession) pull of every reviewed UniProt entry with + both an EC number and a COFACTOR annotation - 111,819 rows as of + 2026-08. Only rows with a structured `Xref=ChEBI:CHEBI:` are + usable (~99.3%); free-text-only rows (~0.7%, e.g. chlorophyll's + photosystem entries) are skipped by preprocess_cofactors.py and are + a separate, deferred LLM-extraction follow-on per the plan doc, not + handled by this pull. Consumed by preprocess_cofactors.py + --uniprot_cofactor_tsv. diff --git a/nextflow/bin/tests/test_get_ec_information_cofactor_merge.py b/nextflow/bin/tests/test_get_ec_information_cofactor_merge.py new file mode 100644 index 0000000..20f30a1 --- /dev/null +++ b/nextflow/bin/tests/test_get_ec_information_cofactor_merge.py @@ -0,0 +1,119 @@ +#!/usr/bin/env python + +""" +get_ec_information.py's cognate_ligands_df concat/groupby block is inline +script code (not a function), so - same approach as +test_get_ec_information_glycan_fallback.py - this replicates its exact +logic against small synthetic dataframes, to check the ligand_source +plumbing added for docs/cofactor_coverage_plan.md without needing the +full pipeline's upstream KEGG/Rhea machinery. + + python3 nextflow/bin/tests/test_get_ec_information_cofactor_merge.py +""" + +import sys +import unittest +from pathlib import Path + +import numpy as np +import pandas as pd +from rdkit import Chem +from rdkit.Chem import PandasTools + +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) + + +def make_frame(rows): + df = pd.DataFrame(rows) + PandasTools.AddMoleculeColumnToFrame(df, smilesCol="smiles", molCol="ROMol") + return df + + +class TestCofactorLigandSourceMerge(unittest.TestCase): + + def setUp(self): + # one reaction-derived compound (unique), and one compound that + # appears via BOTH a reaction source and the cofactor source (same + # canonical SMILES, different EC/source) to check the union case. + self.reaction_df = make_frame([ + {"entry": "1.1.1.1", "compound_id": "C00003", "compound_name": "NAD+", + "smiles": "CC(=O)C", "ligand_db": "KEGG:C00003", "compound_reaction": "R00001"}, + {"entry": "2.2.2.2", "compound_id": "C00099", "compound_name": "some substrate", + "smiles": "CCO", "ligand_db": "KEGG:C00099", "compound_reaction": "R00099"}, + ]).assign(ligand_source="reaction") + + self.cofactor_df = make_frame([ + # same compound as the first reaction row (CC(=O)C), but a + # different EC, sourced from the cofactor path + {"entry": "3.3.3.3", "compound_id": "NAD", "compound_name": "NAD", + "smiles": "CC(=O)C", "ligand_db": "CofactorDB:4", "compound_reaction": ""}, + ]).assign(ligand_source="cofactor") + + def canon_smiles(self, x): + try: + return Chem.MolToSmiles(x, isomericSmiles=False) + except Exception: + return np.nan + + def run_merge_logic(self, reaction_df, cofactor_df): + # exact logic from get_ec_information.py's cognate_ligands_df block + cols = ["entry", "compound_id", "compound_name", "ROMol", "ligand_db", "compound_reaction", "ligand_source"] + cognate_ligands_df = pd.concat([reaction_df[cols], cofactor_df[cols]]) + cognate_ligands_df = cognate_ligands_df.reset_index() + + cognate_ligands_df["canonical_smiles"] = cognate_ligands_df["ROMol"].map( + lambda x: self.canon_smiles(x) if isinstance(x, Chem.rdchem.Mol) else np.nan + ) + cognate_ligands_df_unique_smiles = cognate_ligands_df[ + ["canonical_smiles", "compound_name", "ligand_db", "compound_reaction", "ligand_source"] + ].copy() + cognate_ligands_df_unique_smiles["compound_reaction"] = cognate_ligands_df_unique_smiles["compound_reaction"].fillna("") + cognate_ligands_df_unique_smiles = cognate_ligands_df_unique_smiles.groupby("canonical_smiles", dropna=False).agg( + {"compound_name": set, "ligand_db": set, "compound_reaction": set, "ligand_source": set} + ).reset_index() + cognate_ligands_df_unique_smiles["ligand_db"] = cognate_ligands_df_unique_smiles.ligand_db.str.join("|") + cognate_ligands_df_unique_smiles["compound_name"] = cognate_ligands_df_unique_smiles.compound_name.str.join("|") + cognate_ligands_df_unique_smiles["compound_reaction"] = cognate_ligands_df_unique_smiles.compound_reaction.str.join("|").str.strip("|") + cognate_ligands_df_unique_smiles["ligand_source"] = cognate_ligands_df_unique_smiles.ligand_source.apply(lambda x: "|".join(sorted(x))) + cognate_ligands_df_unique_smiles = cognate_ligands_df_unique_smiles.reset_index(drop=True).reset_index() + cognate_ligands_df_unique_smiles.rename(columns={"index": "uniqueID"}, inplace=True) + + cognate_ligands_df = cognate_ligands_df[["entry", "canonical_smiles"]].merge( + cognate_ligands_df_unique_smiles, on="canonical_smiles", how="left" + ) + return cognate_ligands_df.drop_duplicates() + + def test_reaction_only_compound_keeps_reaction_source(self): + result = self.run_merge_logic(self.reaction_df, pd.DataFrame(columns=self.cofactor_df.columns)) + row = result.loc[result["entry"] == "2.2.2.2"].iloc[0] + self.assertEqual(row["ligand_source"], "reaction") + + def test_compound_shared_between_reaction_and_cofactor_sources_gets_union_tag(self): + result = self.run_merge_logic(self.reaction_df, self.cofactor_df) + # both the reaction-path row (EC 1.1.1.1) and the cofactor-path row + # (EC 3.3.3.3) point at the same canonical_smiles, so both entries + # should carry the combined ligand_source tag. + for ec in ["1.1.1.1", "3.3.3.3"]: + row = result.loc[result["entry"] == ec].iloc[0] + self.assertEqual(row["ligand_source"], "cofactor|reaction") + + def test_cofactor_only_entry_is_present_with_correct_ec(self): + result = self.run_merge_logic(self.reaction_df, self.cofactor_df) + self.assertIn("3.3.3.3", result["entry"].tolist()) + + def test_existing_reaction_rows_are_unaffected_by_cofactor_addition(self): + """docs/cofactor_coverage_plan.md's verification checklist: this + must be a strict addition - the unrelated reaction-only compound's + row should be identical whether or not any cofactor data is + supplied.""" + without_cofactors = self.run_merge_logic(self.reaction_df, pd.DataFrame(columns=self.cofactor_df.columns)) + with_cofactors = self.run_merge_logic(self.reaction_df, self.cofactor_df) + + row_without = without_cofactors.loc[without_cofactors["entry"] == "2.2.2.2"].iloc[0] + row_with = with_cofactors.loc[with_cofactors["entry"] == "2.2.2.2"].iloc[0] + self.assertEqual(row_without["ligand_source"], row_with["ligand_source"]) + self.assertEqual(row_without["ligand_db"], row_with["ligand_db"]) + + +if __name__ == "__main__": + unittest.main() diff --git a/nextflow/bin/tests/test_get_ec_information_glycan_fallback.py b/nextflow/bin/tests/test_get_ec_information_glycan_fallback.py index 386a1bd..bad397d 100644 --- a/nextflow/bin/tests/test_get_ec_information_glycan_fallback.py +++ b/nextflow/bin/tests/test_get_ec_information_glycan_fallback.py @@ -48,6 +48,36 @@ def test_only_missing_rows_with_a_wurcs_value_are_backfilled(self): self.assertTrue(pd.isna(df.loc[2, "smiles"])) self.assertTrue(pd.isna(df.loc[3, "smiles"])) + def test_backfill_works_when_smiles_column_starts_as_float64(self): + """Regression test: in the real pipeline, the "smiles" column is + built via df["smiles"] = df.something.apply(get_smiles_from_csdb), + which pandas can infer as float64 dtype when most/all rows resolve + to NaN (unlike a literal Python list mixing str/NaN, which pandas + infers as object dtype - the case the test above covers, and which + never reproduced this). A later df.loc[mask, "smiles"] = partial assignment into that float64 + column then raises ("Invalid value '...' for dtype + 'float64'") under this pandas version - hit for real running the + full build against real data. get_ec_information.py fixes this by + casting the "smiles" column to object dtype immediately before the + masked assignment; this test locks that in.""" + df = pd.DataFrame({ + "compound_id": ["G1", "G2"], + "wurcs": [CHITOBIOSE_WURCS, np.nan], + }) + # mimics the real pipeline: an .apply() that returns NaN for every + # row infers a float64 column, not object. + df["smiles"] = df["wurcs"].apply(lambda x: np.nan) + self.assertEqual(df["smiles"].dtype, np.float64) + + # exact fix from get_ec_information.py + df["smiles"] = df["smiles"].astype(object) + missing_smiles_mask = df["smiles"].isna() & df["wurcs"].notna() + df.loc[missing_smiles_mask, "smiles"] = df.loc[missing_smiles_mask, "wurcs"].apply(get_smiles_from_wurcs_offline) + + self.assertIsInstance(df.loc[0, "smiles"], str) + self.assertTrue(pd.isna(df.loc[1, "smiles"])) + if __name__ == "__main__": unittest.main() diff --git a/nextflow/bin/tests/test_preprocess_cofactors.py b/nextflow/bin/tests/test_preprocess_cofactors.py new file mode 100644 index 0000000..00eb16f --- /dev/null +++ b/nextflow/bin/tests/test_preprocess_cofactors.py @@ -0,0 +1,266 @@ +#!/usr/bin/env python + +""" +Tests for preprocess_cofactors.py, the new cofactor-EC association +preprocessing script for docs/cofactor_coverage_plan.md. No network +access or real reference data required - all inputs are small synthetic +fixtures built in-memory. + + python3 nextflow/bin/tests/test_preprocess_cofactors.py +""" + +import sys +import tempfile +import unittest +from pathlib import Path + +import pandas as pd +from gemmi import cif + +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) + +from preprocess_cofactors import ( + load_cofactor_db_table, + parse_uniprot_cofactor_df, + classify_and_split_ec_rows, + resolve_chebi_smiles_map, + build_cofactor_ligands_df, +) + + +class TestParseUniprotCofactorDf(unittest.TestCase): + + def test_structured_row_with_single_ec_and_chebi(self): + raw = pd.DataFrame([ + {"Entry": "P00001", "EC number": "1.1.1.1", + "Cofactor": "COFACTOR: Name=NAD(+); Xref=ChEBI:CHEBI:57540;"}, + ]) + result = parse_uniprot_cofactor_df(raw) + self.assertEqual(len(result), 1) + self.assertEqual(result.iloc[0]["ec"], "1.1.1.1") + self.assertEqual(result.iloc[0]["chebi_id"], 57540) + + def test_free_text_only_row_is_dropped(self): + # e.g. chlorophyll's real photosystem entries - Note= with no Xref= + raw = pd.DataFrame([ + {"Entry": "P56766", "EC number": "1.97.1.12", + "Cofactor": "COFACTOR: Note=P700 is a chlorophyll a/chlorophyll a' dimer."}, + ]) + result = parse_uniprot_cofactor_df(raw) + self.assertTrue(result.empty) + + def test_multiple_ecs_on_one_row_are_exploded(self): + raw = pd.DataFrame([ + {"Entry": "P00002", "EC number": "1.1.1.1; 1.1.1.2", + "Cofactor": "COFACTOR: Name=Zn(2+); Xref=ChEBI:CHEBI:29105;"}, + ]) + result = parse_uniprot_cofactor_df(raw) + self.assertCountEqual(result["ec"].tolist(), ["1.1.1.1", "1.1.1.2"]) + + def test_multiple_chebi_xrefs_on_one_row_are_exploded(self): + raw = pd.DataFrame([ + {"Entry": "P00003", "EC number": "1.1.1.1", + "Cofactor": "COFACTOR: Name=[4Fe-4S]; Xref=ChEBI:CHEBI:49883; " + "Name=[2Fe-2S]; Xref=ChEBI:CHEBI:49601;"}, + ]) + result = parse_uniprot_cofactor_df(raw) + self.assertCountEqual(result["chebi_id"].tolist(), [49883, 49601]) + + def test_row_with_no_ec_is_dropped(self): + raw = pd.DataFrame([ + {"Entry": "P00004", "EC number": float("nan"), + "Cofactor": "COFACTOR: Name=NAD(+); Xref=ChEBI:CHEBI:57540;"}, + ]) + result = parse_uniprot_cofactor_df(raw) + self.assertTrue(result.empty) + + def test_duplicate_ec_chebi_pairs_are_deduplicated(self): + raw = pd.DataFrame([ + {"Entry": "P00005", "EC number": "1.1.1.1", + "Cofactor": "COFACTOR: Name=NAD(+); Xref=ChEBI:CHEBI:57540;"}, + {"Entry": "P00006", "EC number": "1.1.1.1", + "Cofactor": "COFACTOR: Name=NAD(+); Xref=ChEBI:CHEBI:57540;"}, + ]) + result = parse_uniprot_cofactor_df(raw) + self.assertEqual(len(result), 1) + + +class TestClassifyAndSplitEcRows(unittest.TestCase): + + def setUp(self): + self.terminal_ec_list = ["1.1.1.1", "1.1.1.10", "1.1.2.1", "2.1.1.1"] + + def test_exact_ec_passes_through_unchanged(self): + df = pd.DataFrame([{"ec": "1.1.1.1", "chebi_id": 1}]) + kept, dropped = classify_and_split_ec_rows(df, self.terminal_ec_list) + self.assertEqual(kept["entry"].tolist(), ["1.1.1.1"]) + self.assertTrue(dropped.empty) + + def test_subsubclass_wildcard_broadcasts_to_matching_terminal_ecs_only(self): + df = pd.DataFrame([{"ec": "1.1.1.-", "chebi_id": 1}]) + kept, dropped = classify_and_split_ec_rows(df, self.terminal_ec_list) + self.assertCountEqual(kept["entry"].tolist(), ["1.1.1.1", "1.1.1.10"]) + self.assertTrue(dropped.empty) + + def test_class_level_wildcard_is_dropped_not_broadcast(self): + """Regression test for the bug caught during planning: naively + broadcasting a class-level wildcard (e.g. "1.-.-.-") would match + literally every terminal EC in that class, incorrectly claiming + structurally unrelated enzymes share a cofactor. This must never + reach `kept` - it belongs in `dropped` instead.""" + df = pd.DataFrame([{"ec": "1.-.-.-", "chebi_id": 1}]) + kept, dropped = classify_and_split_ec_rows(df, self.terminal_ec_list) + self.assertTrue(kept.empty) + self.assertEqual(len(dropped), 1) + self.assertIn("too coarse", dropped.iloc[0]["reason"]) + + def test_subclass_level_wildcard_is_dropped_not_broadcast(self): + df = pd.DataFrame([{"ec": "1.1.-.-", "chebi_id": 1}]) + kept, dropped = classify_and_split_ec_rows(df, self.terminal_ec_list) + self.assertTrue(kept.empty) + self.assertEqual(len(dropped), 1) + + def test_subsubclass_wildcard_with_no_match_is_logged_as_dropped(self): + df = pd.DataFrame([{"ec": "9.9.9.-", "chebi_id": 1}]) + kept, dropped = classify_and_split_ec_rows(df, self.terminal_ec_list) + self.assertTrue(kept.empty) + self.assertEqual(len(dropped), 1) + + def test_mixed_rows_only_safe_ones_survive(self): + df = pd.DataFrame([ + {"ec": "1.1.1.1", "chebi_id": 1}, # exact -> kept + {"ec": "1.1.1.-", "chebi_id": 2}, # subsubclass -> broadcast + {"ec": "1.-.-.-", "chebi_id": 3}, # class-level -> dropped + ]) + kept, dropped = classify_and_split_ec_rows(df, self.terminal_ec_list) + # exact row (1 entry) + broadcast row (2 entries: 1.1.1.1, 1.1.1.10) + self.assertEqual(len(kept), 3) + self.assertEqual(len(dropped), 1) + self.assertEqual(dropped.iloc[0]["chebi_id"], 3) + + +class TestLoadCofactorDbTable(unittest.TestCase): + + def test_merges_ec_csv_with_representative_ccd_from_details_json(self): + with tempfile.TemporaryDirectory() as tmpdir: + csv_path = Path(tmpdir) / "cofactor_ec.csv" + csv_path.write_text('"EC_NO","COFACTOR_ID","SOURCE"\n"1.1.1.1",4,cofactor_db_2010\n"1.1.1.2",4,brenda\n') + + json_path = Path(tmpdir) / "cofactors_details.json" + json_path.write_text('[{"template": "NAD2", "id": 4, "representative": "NAD", "threshold": 0.68}]') + + result = load_cofactor_db_table(str(csv_path), str(json_path)) + self.assertEqual(len(result), 2) + self.assertTrue((result["representative_ccd"] == "NAD").all()) + self.assertCountEqual(result["ec"].tolist(), ["1.1.1.1", "1.1.1.2"]) + + def test_cofactor_class_missing_from_details_json_is_dropped(self): + with tempfile.TemporaryDirectory() as tmpdir: + csv_path = Path(tmpdir) / "cofactor_ec.csv" + csv_path.write_text('"EC_NO","COFACTOR_ID","SOURCE"\n"1.1.1.1",99,cofactor_db_2010\n') + + json_path = Path(tmpdir) / "cofactors_details.json" + json_path.write_text('[{"template": "NAD2", "id": 4, "representative": "NAD", "threshold": 0.68}]') + + result = load_cofactor_db_table(str(csv_path), str(json_path)) + self.assertTrue(result.empty) + + +class TestResolveChebiSmilesMap(unittest.TestCase): + + def test_resolves_known_chebi_id(self): + structures_df = pd.DataFrame({"compound_id": [57540, 29105], "smiles": ["C1=CC...", "[Zn+2]"]}) + names_df = pd.DataFrame({"compound_id": [57540, 29105], "name": ["NAD(+)", "zinc(2+)"]}) + result = resolve_chebi_smiles_map([57540], structures_df, names_df) + self.assertEqual(len(result), 1) + self.assertEqual(result.iloc[0]["smiles"], "C1=CC...") + self.assertEqual(result.iloc[0]["name"], "NAD(+)") + + def test_unknown_chebi_id_is_silently_excluded(self): + structures_df = pd.DataFrame({"compound_id": [57540], "smiles": ["C1=CC..."]}) + names_df = pd.DataFrame({"compound_id": [57540], "name": ["NAD(+)"]}) + result = resolve_chebi_smiles_map([99999999], structures_df, names_df) + self.assertTrue(result.empty) + + +class TestBuildCofactorLigandsDf(unittest.TestCase): + """End-to-end integration across both sources with small synthetic + fixtures - no real reference data or network access.""" + + CCD_TEXT = """ +data_NAD +loop_ +_pdbx_chem_comp_descriptor.comp_id +_pdbx_chem_comp_descriptor.type +_pdbx_chem_comp_descriptor.program +_pdbx_chem_comp_descriptor.program_version +_pdbx_chem_comp_descriptor.descriptor +NAD SMILES 'OpenEye OEToolkits' 2.0.0 CC(=O)C +# +""" + + def setUp(self): + self.ccd_doc = cif.Document() + self.ccd_doc.parse_string(self.CCD_TEXT) + self.terminal_ec_list = ["1.1.1.1", "1.1.1.10", "2.2.2.2"] + + self.cofactor_db_df = pd.DataFrame([ + # subsubclass wildcard -> should broadcast to 1.1.1.1 and 1.1.1.10 + {"ec": "1.1.1.-", "cofactor_id": 4, "representative_ccd": "NAD", "source": "cofactor_db_2010"}, + ]) + self.uniprot_df = pd.DataFrame([ + {"ec": "2.2.2.2", "chebi_id": 29105}, + ]) + self.chebi_structures_df = pd.DataFrame({"compound_id": [29105], "smiles": ["[Zn+2]"]}) + self.chebi_names_df = pd.DataFrame({"compound_id": [29105], "name": ["zinc(2+)"]}) + + def test_output_shape_and_columns(self): + result, dropped = build_cofactor_ligands_df( + self.cofactor_db_df, self.uniprot_df, self.terminal_ec_list, + self.ccd_doc, self.chebi_structures_df, self.chebi_names_df, + ) + expected_cols = {"entry", "compound_id", "compound_name", "ROMol", "ligand_db", "compound_reaction", "ligand_source"} + self.assertEqual(set(result.columns), expected_cols) + self.assertTrue((result["ligand_source"] == "cofactor").all()) + self.assertTrue((result["compound_reaction"] == "").all()) + + def test_cofactor_db_row_broadcasts_to_both_terminal_ecs(self): + result, _ = build_cofactor_ligands_df( + self.cofactor_db_df, self.uniprot_df, self.terminal_ec_list, + self.ccd_doc, self.chebi_structures_df, self.chebi_names_df, + ) + nad_entries = result.loc[result["compound_id"] == "NAD", "entry"].tolist() + self.assertCountEqual(nad_entries, ["1.1.1.1", "1.1.1.10"]) + + def test_uniprot_row_resolves_via_chebi(self): + result, _ = build_cofactor_ligands_df( + self.cofactor_db_df, self.uniprot_df, self.terminal_ec_list, + self.ccd_doc, self.chebi_structures_df, self.chebi_names_df, + ) + zinc_rows = result.loc[result["compound_id"] == "CHEBI:29105"] + self.assertEqual(len(zinc_rows), 1) + self.assertEqual(zinc_rows.iloc[0]["entry"], "2.2.2.2") + self.assertEqual(zinc_rows.iloc[0]["compound_name"], "zinc(2+)") + + def test_no_wildcard_ec_ever_appears_in_entry_column(self): + """Hard invariant the downstream get_pdb_parity.py join depends on + (docs/cofactor_coverage_plan.md's verification checklist).""" + result, _ = build_cofactor_ligands_df( + self.cofactor_db_df, self.uniprot_df, self.terminal_ec_list, + self.ccd_doc, self.chebi_structures_df, self.chebi_names_df, + ) + self.assertFalse(result["entry"].str.contains("-").any()) + + def test_unresolvable_ccd_code_is_excluded_not_erroring(self): + cofactor_db_df = pd.DataFrame([ + {"ec": "1.1.1.1", "cofactor_id": 999, "representative_ccd": "NOPE", "source": "cofactor_db_2010"}, + ]) + result, _ = build_cofactor_ligands_df( + cofactor_db_df, pd.DataFrame(columns=["ec", "chebi_id"]), self.terminal_ec_list, + self.ccd_doc, self.chebi_structures_df, self.chebi_names_df, + ) + self.assertTrue(result.empty) + + +if __name__ == "__main__": + unittest.main() diff --git a/nextflow/bin/tests/test_utils_ec_broadcast.py b/nextflow/bin/tests/test_utils_ec_broadcast.py new file mode 100644 index 0000000..3320b7f --- /dev/null +++ b/nextflow/bin/tests/test_utils_ec_broadcast.py @@ -0,0 +1,101 @@ +#!/usr/bin/env python + +""" +Tests for the EC-completeness classification and subsubclass-level +broadcast helpers added to utils.py for docs/cofactor_coverage_plan.md, +plus the get_chem_comp_descriptors CCD->SMILES resolver (moved here from +process_all_pdb_contacts.py, now reused by preprocess_cofactors.py). + + python3 nextflow/bin/tests/test_utils_ec_broadcast.py +""" + +import sys +import unittest +from pathlib import Path + +from gemmi import cif + +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) + +from utils import classify_ec_completeness, broadcast_subsubclass_ec, get_chem_comp_descriptors + + +class TestClassifyEcCompleteness(unittest.TestCase): + + def test_fully_resolved_ec_is_level_4(self): + self.assertEqual(classify_ec_completeness("1.1.1.10"), 4) + + def test_subsubclass_level_wildcard_is_level_3(self): + self.assertEqual(classify_ec_completeness("1.1.1.-"), 3) + + def test_subclass_level_wildcard_is_level_2(self): + self.assertEqual(classify_ec_completeness("1.1.-.-"), 2) + + def test_class_level_wildcard_is_level_1(self): + self.assertEqual(classify_ec_completeness("1.-.-.-"), 1) + + +class TestBroadcastSubsubclassEc(unittest.TestCase): + + def setUp(self): + self.terminal_ec_list = ["1.1.1.1", "1.1.1.10", "1.1.2.1", "2.1.1.1"] + + def test_matches_only_same_subsubclass(self): + matches = broadcast_subsubclass_ec("1.1.1.-", self.terminal_ec_list) + self.assertCountEqual(matches, ["1.1.1.1", "1.1.1.10"]) + + def test_no_matches_returns_empty_list(self): + matches = broadcast_subsubclass_ec("3.3.3.-", self.terminal_ec_list) + self.assertEqual(matches, []) + + def test_exact_ec_matches_itself_only(self): + matches = broadcast_subsubclass_ec("1.1.1.1", self.terminal_ec_list) + self.assertEqual(matches, ["1.1.1.1"]) + + +class TestGetChemCompDescriptors(unittest.TestCase): + """Reproduces the OpenEye-preference and invalid-SMILES-filtering + branches with a small hand-written CCD-shaped CIF document, rather + than needing the real (multi-GB) ccd.cif.""" + + CCD_TEXT = """ +data_XXX +loop_ +_pdbx_chem_comp_descriptor.comp_id +_pdbx_chem_comp_descriptor.type +_pdbx_chem_comp_descriptor.program +_pdbx_chem_comp_descriptor.program_version +_pdbx_chem_comp_descriptor.descriptor +XXX SMILES ACD 12.0 CCO +XXX SMILES 'OpenEye OEToolkits' 2.0.0 CCO +# +data_YYY +loop_ +_pdbx_chem_comp_descriptor.comp_id +_pdbx_chem_comp_descriptor.type +_pdbx_chem_comp_descriptor.program +_pdbx_chem_comp_descriptor.program_version +_pdbx_chem_comp_descriptor.descriptor +YYY SMILES ACD 12.0 not_a_valid_smiles((( +# +""" + + def setUp(self): + self.ccd_doc = cif.Document() + self.ccd_doc.parse_string(self.CCD_TEXT) + + def test_prefers_openeye_descriptor_when_present(self): + result = get_chem_comp_descriptors(self.ccd_doc, ["XXX"]) + self.assertEqual(result["XXX"], "CCO") + + def test_invalid_smiles_resolves_to_none(self): + result = get_chem_comp_descriptors(self.ccd_doc, ["YYY"]) + self.assertIsNone(result["YYY"]) + + def test_absent_ligand_resolves_to_none(self): + result = get_chem_comp_descriptors(self.ccd_doc, ["ZZZ"]) + self.assertIsNone(result["ZZZ"]) + + +if __name__ == "__main__": + unittest.main() diff --git a/nextflow/bin/utils.py b/nextflow/bin/utils.py index eae404b..db8d34d 100644 --- a/nextflow/bin/utils.py +++ b/nextflow/bin/utils.py @@ -10,6 +10,7 @@ from Bio.ExPASy import Enzyme as EEnzyme from pdbeccdutils.helpers.mol_tools import fix_molecule from rdkit import Chem +from rdkit.Chem import PandasTools import gzip import xml.etree.ElementTree as ET import signal @@ -87,7 +88,7 @@ def get_terminal_record(entry, row, df): return row.ID def get_csdb_from_glycoct(glycoct, cache_df): - if glycoct is np.nan or glycoct == None: + if pd.isna(glycoct): return np.nan elif glycoct in cache_df.glycoct.values: csdb_linear = cache_df.loc[cache_df.glycoct == glycoct, "csdb"].values[0] @@ -110,7 +111,7 @@ def get_csdb_from_glycoct(glycoct, cache_df): return csdb_linear def get_glycoct_from_wurcs(wurcs, cache_df): - if wurcs is np.nan or wurcs == None: + if pd.isna(wurcs): return np.nan elif wurcs in cache_df.WURCS.values: glycoct = cache_df.loc[cache_df.WURCS == wurcs, "glycoct"].values[0] @@ -132,7 +133,7 @@ def get_glycoct_from_wurcs(wurcs, cache_df): return glycoct def get_smiles_from_csdb(csdb_linear, cache_df): - if csdb_linear is np.nan or csdb_linear == None: + if pd.isna(csdb_linear): return np.nan elif csdb_linear in cache_df.csdb.values: smiles = cache_df.loc[cache_df.csdb == csdb_linear, "descriptor"].values[0] @@ -169,7 +170,7 @@ def get_smiles_from_wurcs_offline(wurcs, timeout_seconds = 15): (an ANTLR grammar issue, not specific to this translator), hence the timeout - returns np.nan rather than blocking indefinitely. """ - if wurcs is np.nan or wurcs is None: + if pd.isna(wurcs): return np.nan try: iupac = translate_wurcs_to_iupac(wurcs) @@ -311,6 +312,72 @@ def return_partial_EC_list(ec, total_ec_list): else: return [ec] +def get_chem_comp_descriptors(ccd_doc, comp_id_list): + """Resolve a list of PDB chemical component codes to a single SMILES + descriptor each, using the CCD's own _pdbx_chem_comp_descriptor loop + (OpenEye descriptors preferred when present, else the first + RDKit-parseable SMILES row). Moved here from process_all_pdb_contacts.py + so it can be reused by preprocess_cofactors.py without a cross-script + import - it's a generic CCD-parsing utility, not specific to the + contacts pipeline.""" + ligand_descriptors = {} + for ligand in comp_id_list: + lig_descriptor = None + lig_block = ccd_doc.find_block(ligand) + if lig_block is not None: + lig_descriptors = pd.DataFrame(lig_block.find_mmcif_category("_pdbx_chem_comp_descriptor."), columns = ["comp_id", "type", "program", "program_version", "descriptor"]) + lig_descriptors["descriptor"] = lig_descriptors.descriptor.str.strip("\"|';").str.replace(r"\n$","", regex = True) + lig_descriptors = lig_descriptors.loc[lig_descriptors.type == "SMILES"] + PandasTools.AddMoleculeColumnToFrame(lig_descriptors, smilesCol='descriptor', molCol='pdb_ROMol') + lig_descriptors = lig_descriptors.loc[lig_descriptors.pdb_ROMol.isna() == False] + if len(lig_descriptors) == 0: + lig_descriptor = None + else: + #preference is to use openeye descriptors where available. if not, revert to the first smiles string able to be loaded into RDkit. + preferred_row = lig_descriptors.loc[lig_descriptors.program.str.startswith("OpenEye")] + if not preferred_row.empty: + lig_descriptor = preferred_row.iloc[0].descriptor + else: + # Otherwise, select the first row with a SMILES string + lig_descriptor = lig_descriptors.iloc[0].descriptor + ligand_descriptors[ligand] = lig_descriptor + else: + ligand_descriptors[ligand] = None + return ligand_descriptors + +def classify_ec_completeness(ec): + """Returns how many of an EC number's 4 dot-separated segments are + fully resolved (non "-") before the first wildcard, reading + left-to-right and stopping at the first "-". E.g. "1.1.1.10" -> 4, + "1.1.1.-" -> 3, "1.1.-.-" -> 2, "1.-.-.-" -> 1. Used by + preprocess_cofactors.py to decide whether a source's EC value is exact, + safely broadcastable (subsubclass-level, 3), or too coarse to use at + all (class/subclass-level, 1 or 2) - see docs/cofactor_coverage_plan.md + ("The broadcast problem, and its fix") for why levels 1 and 2 are + excluded: broadcasting them down to every terminal EC in that class/ + subclass would claim structurally unrelated enzymes share a cofactor + just because they share a leading EC digit or two.""" + segments = ec.split(".") + level = 0 + for segment in segments: + if segment == "-": + break + level += 1 + return level + +def broadcast_subsubclass_ec(ec, terminal_ec_list): + """Expands a subsubclass-level partial EC (e.g. "1.1.1.-") to every + matching terminal EC in terminal_ec_list. Reuses the existing + return_partial_EC_list matching logic (already used elsewhere in this + pipeline for SIFTS' partial EC annotations) rather than re-implementing + prefix matching. Callers MUST have already checked + classify_ec_completeness(ec) == 3 before calling this - it will happily + (and, per the cofactor plan, incorrectly) broadcast a class- or + subclass-level wildcard too, since return_partial_EC_list itself has no + concept of "too coarse to use".""" + matches = return_partial_EC_list(ec, terminal_ec_list) + return matches if isinstance(matches, list) else [] + def get_updated_enzyme_records(df, ec_records_df, ec_col = "protein_entity_ec"): ec_list = ec_records_df.ID.unique() ##fill the partial ec records using the original ec ids from the expasy enzyme list