From 0620549452861671e6b9df0a05f7c2adb6a5440e Mon Sep 17 00:00:00 2001 From: Adam Wright Date: Mon, 14 Sep 2026 19:42:41 -0400 Subject: [PATCH 1/2] Ship the cofactor list with the pathway artifacts, derived not typed MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A consumer deciding whether a perturbation may travel through ATP needs to know which nodes are ATP, and that knowledge has to travel with the networks. It did not, so both this repo and DeltaSignal kept their own copy, and the copies drifted. Audited against Release97, 6 of the 13 hand-written ids here are wrong: R-ALL-29438 commented "PPi" is GTP, R-ALL-29390 commented "Pi variant" is PXLP (pyridoxal 5'-phosphate), R-ALL-29360 commented "ADP variant" is NAD+, and R-ALL-217093 / R-ALL-110114 / R-ALL-29986 do not exist in the release at all. Nothing caught any of it, because the list is only ever read here. So derive it. `get_cofactor_species` returns every SimpleEntity whose ReferenceMolecule carries one of 28 curated ChEBI identifiers — 253 species at Release97, every compartment variant included. ChEBI identity is the only stable handle: a name query misses compartment variants and matches by substring, and stable ids go stale between releases. This list managed both. `export_cofactors` writes cofactors.csv into each pathway directory alongside the networks, carrying stable_id, molecule, chebi_id, name, in_network and reactome_release. Every known cofactor is listed rather than only the present ones, so a consumer can tell an empty intersection from a bundle that predates the file, and the file is written even when a pathway contains none. The networks themselves are unchanged: this adds an artifact, and the existing _COFACTOR_STIDS bridge suppression is left exactly as it is. Correcting those six entries changes which bridges are suppressed, which changes the networks, which needs its own A/B — a separate change, deliberately not bundled here. Co-Authored-By: Claude Opus 5 (1M context) --- src/logic_network_generator.py | 62 +++++++++++++++++++++++++++ src/neo4j_connector.py | 73 ++++++++++++++++++++++++++++++++ src/pathway_generator.py | 8 ++++ tests/test_provenance_exports.py | 48 +++++++++++++++++++++ 4 files changed, 191 insertions(+) diff --git a/src/logic_network_generator.py b/src/logic_network_generator.py index b49cde1..b1f6e87 100755 --- a/src/logic_network_generator.py +++ b/src/logic_network_generator.py @@ -2442,6 +2442,68 @@ def export_node_reaction_context(entity_uuid_registry: Dict[tuple, str], logger.info(f"Exported {len(rows)} node-reaction-context rows: {output_file}") +def export_cofactors(pathway_logic_network: pd.DataFrame, + reactome_id_to_uuid: Dict[str, str], + output_file: str) -> None: + """Write cofactors.csv — which nodes in THIS network are metabolic cofactors. + + The network keeps every participant Reactome records, cofactors included: + the generator represents pathways as curators intended. But a consumer + deciding whether a perturbation may travel through ATP needs to know which + nodes those are, and that knowledge must travel WITH the artifacts. A + consumer holding its own copy of the list is a copy that silently goes + stale — the generator and DeltaSignal each kept one and they had diverged, + with six of thirteen entries on this side stale or mislabelled before the + set was derived rather than typed. + + So the list ships in the bundle. Pull a pathway out of S3 and it carries + its own answer, pinned to the release it was generated from. + + Output CSV columns: + - stable_id: the Reactome stable ID of the cofactor species + - molecule: which cofactor it is (ATP, NAD+, Pi, …) + - chebi_id: the ChEBI identifier the membership was derived from + - name: the release's display name, including compartment + - in_network: 1 if this species appears as a node in this pathway + - reactome_release: the release this was derived from + + Every known cofactor is listed, not only those present, so a consumer can + tell an empty intersection from a missing file. The file is written even + when the pathway contains none. + """ + from src.neo4j_connector import get_cofactor_species, get_reactome_release + + species = get_cofactor_species() + release = get_reactome_release() + + present: set[str] = set() + if not pathway_logic_network.empty: + uuids = set(pathway_logic_network["source_id"].dropna().unique()) + uuids.update(pathway_logic_network["target_id"].dropna().unique()) + for stable_id, uuid in reactome_id_to_uuid.items(): + if uuid in uuids: + present.add(stable_id) + + present_count = sum(1 for e in species if e["stable_id"] in present) + rows = [ + { + "stable_id": entry["stable_id"], + "molecule": entry["molecule"], + "chebi_id": entry["chebi_id"], + "name": entry["name"], + "in_network": 1 if entry["stable_id"] in present else 0, + "reactome_release": release if release is not None else "", + } + for entry in species + ] + pd.DataFrame(rows, columns=["stable_id", "molecule", "chebi_id", "name", + "in_network", "reactome_release"]).to_csv( + output_file, index=False) + logger.info( + f"Exported {len(rows)} cofactor species " + f"({present_count} present in this network) to {output_file}") + + def export_node_resolution(pathway_id: str, pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFrame, diff --git a/src/neo4j_connector.py b/src/neo4j_connector.py index 637366a..69f19e4 100755 --- a/src/neo4j_connector.py +++ b/src/neo4j_connector.py @@ -720,6 +720,79 @@ def get_modifier_isoform_entity_set_ids() -> Set[str]: return _modifier_isoform_set_cache +# Metabolic cofactors, keyed by ChEBI rather than by name or stable id. +# +# ChEBI identity is the only stable handle here. A name query silently misses +# compartment variants and matches by substring ("phosphate" pulls in pyridoxal +# 5'-phosphate); stable ids go stale between releases. An audit of the +# thirteen hand-written ids this module used to carry found six that were +# stale or mislabelled against Release97 — R-ALL-29438 commented "PPi" is +# GTP, R-ALL-29390 commented "Pi variant" is PXLP, R-ALL-29360 commented +# "ADP variant" is NAD+, and three did not exist at all. +# +# Energy and phosphate carriers, redox pairs, one-carbon donors, water, +# dissolved gases and bulk ions. Deliberately NOT here: second messengers +# (Ca2+, PIP3, PI(4,5)P2, cAMP, cGMP, DAG, IP3) and modifier tags whose +# transfer is the regulatory event (ubiquitin, SUMO) — in a signalling pathway +# those ARE the signal. +_COFACTOR_CHEBI: Dict[str, List[str]] = { + "ATP": ["30616"], "ADP": ["456216"], "AMP": ["456215"], + "GTP": ["37565"], "GDP": ["58189"], "GMP": ["58115"], + "CTP": ["37563"], "CDP": ["58069"], + "UTP": ["46398"], "UDP": ["17659", "58223"], + "NAD+": ["57540"], "NADH": ["57945"], + "NADP+": ["18009", "58349"], "NADPH": ["16474", "57783"], + "FAD": ["57692"], "FADH2": ["58307"], + "CoA-SH": ["57287"], "AdoMet": ["59789"], "AdoHcy": ["57856"], + "H2O": ["15377"], "H+": ["15378"], + "Pi": ["43474"], "PPi": ["33019"], + "O2": ["15379"], "CO2": ["16526"], + "Na+": ["29101"], "K+": ["29103"], "Cl-": ["17996"], +} + +_cofactor_cache: Optional[List[Dict[str, str]]] = None + + +def get_cofactor_species() -> List[Dict[str, str]]: + """Every SimpleEntity in the connected release that IS one of the cofactors. + + Returns one dict per species with ``stable_id``, ``molecule``, ``chebi_id`` + and ``name``, covering every compartment variant the release defines. + Derived rather than hand-maintained so a new compartment appears by itself + and a renamed or retired stable id disappears by itself. + + Cached for the process. Raises if Neo4j is unreachable. + """ + global _cofactor_cache + if _cofactor_cache is not None: + return _cofactor_cache + query = """ + MATCH (se:SimpleEntity)-[:referenceEntity]->(rm:ReferenceMolecule) + WHERE rm.identifier IN $ids + RETURN DISTINCT se.stId AS stable_id, se.displayName AS name, + rm.identifier AS chebi_id + ORDER BY stable_id + """ + by_chebi = {c: mol for mol, ids in _COFACTOR_CHEBI.items() for c in ids} + try: + rows = get_graph().run(query, ids=list(by_chebi)).data() + except Exception: + logger.error("Error in get_cofactor_species", **_traceback_kwargs()) + raise + out = [ + { + "stable_id": r["stable_id"], + "molecule": by_chebi[r["chebi_id"]], + "chebi_id": r["chebi_id"], + "name": r["name"] or "", + } + for r in rows + if r.get("stable_id") + ] + _cofactor_cache = sorted(out, key=lambda d: (d["molecule"], d["stable_id"])) + return _cofactor_cache + + def get_reference_entity_id(entity_id: str) -> Union[str, None]: if entity_id in _reference_entity_cache: return _reference_entity_cache[entity_id] diff --git a/src/pathway_generator.py b/src/pathway_generator.py index 10e89f1..e594c6c 100755 --- a/src/pathway_generator.py +++ b/src/pathway_generator.py @@ -11,6 +11,7 @@ from src.decomposed_uid_mapping import decomposed_uid_mapping_column_types from src.logic_network_generator import ( create_pathway_logic_network, + export_cofactors, export_entity_reaction_proxy_mapping, export_node_reaction_context, export_node_resolution, @@ -411,6 +412,13 @@ def generate_pathway_file( str(pathway_output_dir / "node_resolution.csv"), str(pathway_output_dir / "node_exclusions.csv"), ) + # Ships WITH the networks so an artifact bundle pulled from S3 + # carries its own answer to "which of these nodes is ATP". + export_cofactors( + result.logic_network, + result.uuid_mapping, + str(pathway_output_dir / "cofactors.csv"), + ) except Exception as e: logger.error(f"Failed to write node provenance files: {e}", exc_info=True) # Don't raise - supplementary diff --git a/tests/test_provenance_exports.py b/tests/test_provenance_exports.py index 8148d8c..cbc699e 100644 --- a/tests/test_provenance_exports.py +++ b/tests/test_provenance_exports.py @@ -187,3 +187,51 @@ def test_glyph_id_and_diagram_are_written_together(tmp_path, monkeypatch): f"glyph_id={row['glyph_id']!r} and diagram_stid={row['diagram_stid']!r} " "must be present together or absent together" ) + + +def test_export_cofactors_lists_all_and_flags_present(tmp_path, monkeypatch): + """Every known cofactor is listed; only those in the network are flagged. + + Listing all of them is what lets a consumer tell "this pathway has no + cofactors" from "this bundle predates the file". + """ + monkeypatch.setattr(neo4j_connector, "get_cofactor_species", lambda: [ + {"stable_id": "R-ALL-113592", "molecule": "ATP", "chebi_id": "30616", + "name": "ATP [cytosol]"}, + {"stable_id": "R-ALL-29356", "molecule": "H2O", "chebi_id": "15377", + "name": "H2O [cytosol]"}, + ]) + monkeypatch.setattr(neo4j_connector, "get_reactome_release", lambda: 97) + + edges = pd.DataFrame([{"source_id": "u-atp", "target_id": "u-rxn"}]) + out = tmp_path / "cofactors.csv" + m.export_cofactors(edges, {"R-ALL-113592": "u-atp"}, str(out)) + + df = pd.read_csv(out) + assert list(df.columns) == ["stable_id", "molecule", "chebi_id", "name", + "in_network", "reactome_release"] + assert len(df) == 2, "every known cofactor is listed, not only the present ones" + assert set(df.loc[df.in_network == 1, "stable_id"]) == {"R-ALL-113592"} + assert set(df["reactome_release"]) == {97}, "the release must travel with the list" + + +def test_export_cofactors_writes_a_file_even_when_none_are_present(tmp_path, monkeypatch): + """A pathway with no cofactors still gets the file, all flags zero. + + A missing file and an empty intersection mean different things and a + consumer must be able to distinguish them. + """ + monkeypatch.setattr(neo4j_connector, "get_cofactor_species", lambda: [ + {"stable_id": "R-ALL-113592", "molecule": "ATP", "chebi_id": "30616", + "name": "ATP [cytosol]"}, + ]) + monkeypatch.setattr(neo4j_connector, "get_reactome_release", lambda: 97) + + edges = pd.DataFrame([{"source_id": "u-x", "target_id": "u-y"}]) + out = tmp_path / "cofactors.csv" + m.export_cofactors(edges, {"R-HSA-9999": "u-x"}, str(out)) + + df = pd.read_csv(out) + assert out.exists() + assert len(df) == 1 + assert int(df.in_network.sum()) == 0 From 5275551b190ba6080b2de553d17388336dcb9c29 Mon Sep 17 00:00:00 2001 From: Adam Wright Date: Tue, 15 Sep 2026 11:02:43 -0400 Subject: [PATCH 2/2] Fix export_cofactors reading the uuid mapping backwards MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Found by adversarial review. `reactome_id_to_uuid` is built as `reactome_id_to_uuid[entity_uuid] = entity_stId` — it is uuid -> stId despite its name, which is why `_uuid_to_stable_id_map` exists and why every other consumer calls `.get(uuid)`. `export_cofactors` was the only new export not using that helper: it unpacked `for stable_id, uuid in ...items()`, binding the names to the wrong halves, then tested a stId against a set of uuids. That is never true, so `present` stayed empty and EVERY row of EVERY cofactors.csv would have shipped with in_network=0 — a file asserting that no pathway contains any cofactor. The same scan was wrong a second way: one stable id routinely maps to several uuids (the positional-decomposition silo — GPVI carries four separate GTP nodes), and a dict keyed by stable id holds only one of them, so a split entity was undercounted even in the right direction. Both existing tests passed the mapping in the direction production never uses, so they locked the bug in rather than catching it. The two new tests exercise both directions and a split entity, and both fail against the previous code — verified by reverting the fix. mypy and ruff clean. pytest: 955 passed; the 8 failures are identical before and after this change (6 are a local scipy stub for this box's Python 3.14, 1 needs Neo4j, 1 pre-exists on main). Co-Authored-By: Claude Opus 5 (1M context) --- src/logic_network_generator.py | 17 +++++++---- tests/test_provenance_exports.py | 50 ++++++++++++++++++++++++++++++++ 2 files changed, 62 insertions(+), 5 deletions(-) diff --git a/src/logic_network_generator.py b/src/logic_network_generator.py index b1f6e87..e135ed7 100755 --- a/src/logic_network_generator.py +++ b/src/logic_network_generator.py @@ -2476,13 +2476,20 @@ def export_cofactors(pathway_logic_network: pd.DataFrame, species = get_cofactor_species() release = get_reactome_release() + # `reactome_id_to_uuid` can be stored in either direction, and one stable + # id routinely maps to SEVERAL uuids (that is the positional-decomposition + # silo — GPVI carries four separate GTP nodes). Both facts break a naive + # `for stable_id, uuid in mapping.items()` scan: the wrong direction marks + # every row absent, and the right one still undercounts a split entity. + # `_uuid_to_stable_id_map` already solves both and is what the other + # exporters use. present: set[str] = set() if not pathway_logic_network.empty: - uuids = set(pathway_logic_network["source_id"].dropna().unique()) - uuids.update(pathway_logic_network["target_id"].dropna().unique()) - for stable_id, uuid in reactome_id_to_uuid.items(): - if uuid in uuids: - present.add(stable_id) + for node_id in _uuid_to_stable_id_map( + pathway_logic_network, reactome_id_to_uuid).values(): + # A set_variant node is "{parent}::variant::{members}"; a cofactor + # is always a plain stId, so the split is a cheap exact match. + present.add(node_id) present_count = sum(1 for e in species if e["stable_id"] in present) rows = [ diff --git a/tests/test_provenance_exports.py b/tests/test_provenance_exports.py index cbc699e..17f533b 100644 --- a/tests/test_provenance_exports.py +++ b/tests/test_provenance_exports.py @@ -235,3 +235,53 @@ def test_export_cofactors_writes_a_file_even_when_none_are_present(tmp_path, mon assert out.exists() assert len(df) == 1 assert int(df.in_network.sum()) == 0 + + +def test_export_cofactors_handles_both_mapping_directions(tmp_path, monkeypatch): + """`reactome_id_to_uuid` is stored either direction depending on caller. + + The first version of this exporter assumed stable_id -> uuid. Given the + other direction it marked EVERY row absent and shipped a file saying no + pathway contains any cofactor. The original test constructed the mapping in + the assumed direction, so it passed either way and could not catch this. + """ + monkeypatch.setattr(neo4j_connector, "get_cofactor_species", lambda: [ + {"stable_id": "R-ALL-113592", "molecule": "ATP", "chebi_id": "30616", + "name": "ATP [cytosol]"}, + ]) + monkeypatch.setattr(neo4j_connector, "get_reactome_release", lambda: 97) + + uuid = "aaaaaaaa-0000-0000-0000-000000000001" + edges = pd.DataFrame([{"source_id": uuid, "target_id": "u-rxn"}]) + + for direction, mapping in ( + ("stable_id -> uuid", {"R-ALL-113592": uuid}), + ("uuid -> stable_id", {uuid: "R-ALL-113592"}), + ): + out = tmp_path / f"cofactors_{direction.split()[0]}.csv" + m.export_cofactors(edges, mapping, str(out)) + df = pd.read_csv(out) + assert int(df.in_network.sum()) == 1, f"ATP missed with {direction}" + + +def test_export_cofactors_finds_an_entity_split_across_uuids(tmp_path, monkeypatch): + """One stable id routinely maps to several uuids (the silo). + + A dict keyed by stable id holds only one of them, so scanning the mapping + by key undercounts a split entity. GPVI carries four separate GTP nodes. + """ + monkeypatch.setattr(neo4j_connector, "get_cofactor_species", lambda: [ + {"stable_id": "R-ALL-29438", "molecule": "GTP", "chebi_id": "37565", + "name": "GTP [cytosol]"}, + ]) + monkeypatch.setattr(neo4j_connector, "get_reactome_release", lambda: 97) + + u1 = "aaaaaaaa-0000-0000-0000-00000000000a" + u2 = "aaaaaaaa-0000-0000-0000-00000000000b" + # Only the SECOND occurrence appears in the network. + edges = pd.DataFrame([{"source_id": u2, "target_id": "u-rxn"}]) + out = tmp_path / "cofactors.csv" + m.export_cofactors(edges, {u1: "R-ALL-29438", u2: "R-ALL-29438"}, str(out)) + + df = pd.read_csv(out) + assert int(df.in_network.sum()) == 1, "split entity missed when only one uuid is used"