Add cofactor-EC coverage as a second, independent cognate-ligand source

README.md

@@ -215,15 +215,32 @@ 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/)
 - [RHEA](https://www.rhea-db.org/)
 - [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
 

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

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:<n>` (≈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.