From cf7424d832b196a589c78882e66a6604b6d0732d Mon Sep 17 00:00:00 2001 From: m-crown <30960941+m-crown@users.noreply.github.com> Date: Sun, 9 Aug 2026 20:37:26 +0200 Subject: [PATCH 1/2] Add cofactor-EC coverage as a second, independent cognate-ligand source Combines CoFactor DB 2010 + BRENDA (vendored from RelLig's cofactor_ec.csv and cofactors_details.json) with UniProt's structured COFACTOR annotations into cofactor_ligands_df.pkl (preprocess_cofactors.py, mirrors preprocess_rhea.py's pattern), then concats it into get_ec_information.py's cognate_ligands_df alongside the existing reaction-derived sources, tagged via a new ligand_source column. Partial ECs are handled per the plan doc's broadcast rule: subsubclass-level wildcards (N.N.N.-) are expanded to matching terminal ECs, class/subclass-level wildcards are dropped rather than broadcast (verified unsafe - would match every terminal EC in that class). See docs/cofactor_coverage_plan.md for the full design and the coverage numbers behind it. Co-Authored-By: Claude Sonnet 5 --- docs/cofactor_coverage_plan.md | 302 ++++++++++++++++++ docs/installation.md | 18 ++ nextflow/bin/get_ec_information.py | 30 +- nextflow/bin/preprocess_cofactors.py | 261 +++++++++++++++ nextflow/bin/process_all_pdb_contacts.py | 28 +- nextflow/bin/reference_data_manifest.yaml | 64 +++- .../test_get_ec_information_cofactor_merge.py | 119 +++++++ .../bin/tests/test_preprocess_cofactors.py | 266 +++++++++++++++ nextflow/bin/tests/test_utils_ec_broadcast.py | 101 ++++++ nextflow/bin/utils.py | 67 ++++ 10 files changed, 1206 insertions(+), 50 deletions(-) create mode 100644 docs/cofactor_coverage_plan.md create mode 100644 nextflow/bin/preprocess_cofactors.py create mode 100644 nextflow/bin/tests/test_get_ec_information_cofactor_merge.py create mode 100644 nextflow/bin/tests/test_preprocess_cofactors.py create mode 100644 nextflow/bin/tests/test_utils_ec_broadcast.py diff --git a/docs/cofactor_coverage_plan.md b/docs/cofactor_coverage_plan.md new file mode 100644 index 0000000..cfe7611 --- /dev/null +++ b/docs/cofactor_coverage_plan.md @@ -0,0 +1,302 @@ +# 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. + +## Suggested order of work + +1. Vendor and parse `cofactor_ec.csv` / `cofactors_details.json`; resolve + representative CCD codes to SMILES. Smallest, self-contained first + step. +2. Add the UniProt bulk cofactor TSV pull to the reference data manifest; + parse structured rows; implement the exact/broadcast/drop split by EC + completeness. +3. Add `ligand_source` column (backfilling existing reaction rows to + `"reaction"`); union and concat into `cognate_ligands_df` construction. +4. Re-run the unmatched-ligand analysis to quantify real improvement + against the ~4,843-EC combined coverage. +5. Only then, as a separate follow-on: LLM-assisted extraction for the 742 + free-text-only UniProt rows (chlorophyll's actual fix). 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/nextflow/bin/get_ec_information.py b/nextflow/bin/get_ec_information.py index 3568409..6c10a45 100644 --- a/nextflow/bin/get_ec_information.py +++ b/nextflow/bin/get_ec_information.py @@ -365,6 +365,7 @@ def main(): 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('--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: @@ -499,6 +500,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 @@ -672,23 +684,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) 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/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/reference_data_manifest.yaml b/nextflow/bin/reference_data_manifest.yaml index 935f365..f5516d6 100644 --- a/nextflow/bin/reference_data_manifest.yaml +++ b/nextflow/bin/reference_data_manifest.yaml @@ -326,29 +326,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: > + 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: For the cofactor-coverage plan (docs/cofactor_coverage_plan.md), not the current pipeline. + 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_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..df81e65 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 @@ -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 From 5fdc7a2c2cccec3e02a427169bc24d0248619794 Mon Sep 17 00:00:00 2001 From: m-crown <30960941+m-crown@users.noreply.github.com> Date: Sun, 9 Aug 2026 23:56:11 +0200 Subject: [PATCH 2/2] Fix three real-data bugs in preprocess_rhea.py blocking Rhea preprocessing Found while running preprocess_rhea.py end-to-end for the first time against real reference data (unrelated to the cofactor-coverage work, but this script is a prerequisite for it). None of these are cofactor- specific: - rdfreader's .metadata property fixed-column-parses a mol block's whole header, including the program/timestamp line - Rhea's CDK-written molfiles don't conform to that fixed-width format (e.g. " CDK 2/12/10,15:27" instead of MDL's fixed MMDDYY field), which broke metadata parsing on ~97% of real Rhea rd/ files. Read the molecule name directly from the mol block's raw first line instead, since that's all this script actually needed from .metadata. - Assigning int -1 into a pandas StringDtype column raises under modern pandas (silently allowed under older/object-dtype behavior). Assign the string "-1" instead - the column gets cast to int on the next line regardless. - ChEBI's names.tsv.gz moved to lowercase columns upstream (already documented in download_reference_data.py's derive_chebi_results docstring, but never applied here) - rename after reading. Co-Authored-By: Claude Sonnet 5 --- nextflow/bin/preprocess_rhea.py | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) 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)