diff --git a/src/logic_network_generator.py b/src/logic_network_generator.py index b49cde1..e135ed7 100755 --- a/src/logic_network_generator.py +++ b/src/logic_network_generator.py @@ -2442,6 +2442,75 @@ 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() + + # `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: + 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 = [ + { + "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..17f533b 100644 --- a/tests/test_provenance_exports.py +++ b/tests/test_provenance_exports.py @@ -187,3 +187,101 @@ 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 + + +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"