Ship the cofactor list with the pathway artifacts, derived not typed

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,

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]

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

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"