diff --git a/src/logic_network_generator.py b/src/logic_network_generator.py index 25ea14b..ba51f0e 100755 --- a/src/logic_network_generator.py +++ b/src/logic_network_generator.py @@ -3409,6 +3409,421 @@ def parts(node_id: str): logger.info(f"Exported {len(rows)} drug-derived entities to {output_file}") +POOL_MAX_STEPS = 6 # a state -> intermediates -> state path longer than this is dropped +POOL_PATH_BUDGET = 20000 # a pool with more paths (node + reaction sequences) than this is dropped and counted +POOL_STATS_KEYS = ("pools", "states", "intermediates", "paths", "multi_step_pools", "long_paths_dropped", + "shared_nodes_dropped", "ties", "carriers", "autocat_source", "autocat_product", + "autocat_other", "carrier_loops", "no_path_pools", "off_path_intermediates", + "off_path_states", "ambiguous_copies", "budget_dropped", "merged_proteins", "copy_rows", + "carrier_conflicts") + + +def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFrame, + uuid_mapping: Dict[str, str], steps, profiles: Dict[str, Dict[str, Any]], + donor_reactions: Optional[Set[str]] = None, + set_members: Optional[Dict[str, Set[str]]] = None, + stats: Optional[Dict[str, int]] = None): + """Modification-cycle pools at NODE level (deltasignal specs/039 amendment 2). + + ``steps`` are the pathway's R-steps at stId level, ``(reaction, R, input + form, output form, enzyme_driven)``: a reaction whose exactly one input and + exactly one output contain the reference entity R, both at stoichiometry 1, + and which differ (see :func:`r_steps`). ``profiles`` is + :func:`get_form_profiles` over every node's entity (and a variant's chosen + members); ``set_members`` says which nodes a set-valued participant became. + + A step exists at node level only as a reaction node whose R-containing + input node and output node map to the step's forms (the node's entity, its + variant parent, or a member of the step's set). Per R, the R-graph over + form nodes is condensed. In a strongly connected set with >= 2 forms, the + CORE is the proteins every form carries; only a form made of exactly the + core (small molecules allowed) can be a state, and the set is a pool only + if its core-only forms carry >= 2 modification signatures (amendment 3; + one signature is a carrier loop: an enzyme's E -> E:S -> E, and a set + whose forms differ only by what else is bound is not a pool). In each + signature the core-only forms with the fewest non-small slots are states, + the rest intermediates. Transitions are + simple directed paths state -> intermediates* -> another state of at most + ``POOL_MAX_STEPS`` steps, longer ones dropped and counted. A path's + identity is its node sequence plus the reaction stId at each step; the + reaction NODES of that stId between the step's two nodes (variant and set + copies) are the step's parallel copies and are listed one row each, so a + consumer reads a step as the mean of its copies and a path as the product + of its steps. Different reactions between the same two nodes (GEF + exchange beside intrinsic exchange) stay different paths, since each is + weighed on its own as enzyme-driven or not. A pool with more than + ``POOL_PATH_BUDGET`` paths is dropped and counted (``budget_dropped``). The base state follows amendment 1 + (residues, then the donor-consuming direction, then components; undecided + pools fall back to the smaller stId and are counted). Pools identical for + several R (the members of one set travelling together) are one pool; a node + claimed by two DIFFERENT pools is removed from both, everything is + recomputed, and the node is counted. Carriers are the enzyme's free-form + nodes (a joining input of an enzyme-driven step, containing a protein and + not R, that a step of the pool releases) with the reaction nodes releasing + them; a carrier that any pool manages (state, intermediate or step + node) is excluded and counted (``carrier_conflicts``). + + Returns (forms, transitions, carriers): + forms: [(pool_id, node_uuid, stable_id, role, is_base)], role state | intermediate + transitions: [(pool_id, path_id, step, source_uuid, target_uuid, reaction_uuid, reaction_stid, enzyme_driven)] + one row per (path, step, copy): rows sharing (pool_id, path_id, step) have the same + source, target and reaction_stid and differ in reaction_uuid. ``enzyme_driven`` is the + STEP's flag; a path is enzyme-driven if any step is. + carriers: [(pool_id, carrier_uuid, release_reaction_uuid)] + ``stats`` receives the counts in ``POOL_STATS_KEYS``; the autocatalysis + counts are distinct (reaction node, catalyst) pairs where a form of the pool + catalyses a step of its own pool, split by whether the catalyst is the + step's source form, its product form, or another form. + """ + import networkx as nx # type: ignore[import-untyped] + stats = stats if stats is not None else {} + for k in POOL_STATS_KEYS: + stats[k] = 0 + donor_reactions = donor_reactions or set() + set_members = set_members or {} + steps = list(steps) + if pathway_logic_network.empty or not steps: + return [], [], [] + ent = _uuid_to_stable_id_map(pathway_logic_network, uuid_mapping) + vr = {str(u): str(r) for u, r in zip(reaction_id_map["uid"], reaction_id_map["reactome_id"])} + ins: Dict[str, Dict[str, float]] = {} + outs: Dict[str, Dict[str, float]] = {} + cats: Dict[str, Set[str]] = {} + has_st = "stoichiometry" in pathway_logic_network.columns + for _, e in pathway_logic_network.iterrows(): + s, t, et = str(e["source_id"]), str(e["target_id"]), e.get("edge_type") + st = 1.0 + if has_st: + try: + st = float(e["stoichiometry"]) + except (TypeError, ValueError): + st = 1.0 + if st != st: # NaN + st = 1.0 + if et == "input" and t in vr: + ins.setdefault(t, {})[s] = st + elif et == "output" and s in vr: + outs.setdefault(s, {})[t] = st + elif et == "catalyst" and t in vr: + cats.setdefault(t, set()).add(s) + copies: Dict[str, List[str]] = {} + for u, r in vr.items(): + copies.setdefault(r, []).append(u) + + def parts(u: str) -> Tuple[str, List[str]]: + s = ent.get(u, "") + if "::variant::" not in s: + return s, [] + parent, members = s.split("::variant::", 1) + return parent, [m for m in members.split("_") if m] + + def pe_of(u: str) -> str: + return parts(u)[0] + + prot_cache: Dict[str, frozenset] = {} + + def proteins(u: str) -> frozenset: + if u not in prot_cache: + parent, members = parts(u) + if not members: + out = set(profiles.get(parent, {}).get("proteins", ())) + else: + out = set(profiles.get(parent, {}).get("fixed", ())) + for m in members: + out |= set(profiles.get(m, {}).get("proteins", ())) + prot_cache[u] = frozenset(out) + return prot_cache[u] + + def prof(u: str, key: str, default=0): + return profiles.get(pe_of(u), {}).get(key, default) + + def sig(u: str, r): + return profiles.get(pe_of(u), {}).get("sig", {}).get(r, frozenset()) + + def maps(u: str, pe: str) -> bool: + p = pe_of(u) + return p == pe or p in set_members.get(pe, ()) + + # uuid-level steps: (R, in node, out node, reaction node) -> (reaction stId, enzyme_driven) + ustep: Dict[Tuple[Any, str, str, str], Tuple[str, bool]] = {} + for rx, r, a, b, enz in steps: + for rx_u in sorted(copies.get(rx, [])): + cin = [u for u, st in ins.get(rx_u, {}).items() if st == 1 and maps(u, a) and r in proteins(u)] + cout = [u for u, st in outs.get(rx_u, {}).items() if st == 1 and maps(u, b) and r in proteins(u)] + if len(cin) != 1 or len(cout) != 1: + if len(cin) > 1 or len(cout) > 1: + stats["ambiguous_copies"] += 1 + continue + if cin[0] != cout[0]: + ustep[(r, cin[0], cout[0], rx_u)] = (rx, bool(enz)) + if not ustep: + return [], [], [] + + def build(removed: Set[str]): + local = {k: 0 for k in POOL_STATS_KEYS} + by_r: Dict[Any, List[Tuple[str, str, str, str, bool]]] = {} + for (r, a_u, b_u, rx_u), (rs, enz) in ustep.items(): + if a_u in removed or b_u in removed or rx_u in removed: + continue + by_r.setdefault(r, []).append((a_u, b_u, rx_u, rs, enz)) + pools: Dict[tuple, Dict[str, Any]] = {} + for r in sorted(by_r): + g = nx.DiGraph() + g.add_edges_from((a_u, b_u) for a_u, b_u, _, _, _ in by_r[r]) + for scc in nx.strongly_connected_components(g): + if len(scc) < 2: + continue + # amendment 3: the core is the proteins every form of the set carries; + # only a form made of exactly the core (small molecules allowed) can + # be a state, and the set is a pool only if such forms carry >= 2 + # signatures. An enzyme's loop has one core-only form: not a pool. + core = frozenset.intersection(*(proteins(u) for u in scc)) + by_sig: Dict[frozenset, List[str]] = {} + for u in scc: + if proteins(u) == core: + by_sig.setdefault(sig(u, r), []).append(u) + if len(by_sig) < 2: + local["carrier_loops"] += 1 + continue + states: Set[str] = set() + for members in by_sig.values(): + fewest = min(prof(u, "slots") for u in members) + states |= {u for u in members if prof(u, "slots") == fewest} + # a step is (source node, target node, reaction stId); the reaction + # nodes of that stId between those nodes are its parallel copies + copies_of: Dict[Tuple[str, str, str], List[str]] = {} + enz_of: Dict[Tuple[str, str, str], bool] = {} + for a_u, b_u, rx_u, rs, enz in sorted(by_r[r]): + if a_u in scc and b_u in scc: + copies_of.setdefault((a_u, b_u, rs), []).append(rx_u) + enz_of[(a_u, b_u, rs)] = enz_of.get((a_u, b_u, rs), False) or enz + adj: Dict[str, List[Tuple[str, str]]] = {} + for a_u, b_u, rs in sorted(copies_of): + adj.setdefault(a_u, []).append((b_u, rs)) + paths: List[Tuple[Tuple[str, str, str, bool, Tuple[str, ...]], ...]] = [] + over_budget = False + + def walk(u, so_far, seen): + nonlocal over_budget + for v, rs in adj.get(u, ()): + if over_budget: + return + if v in seen: + continue + new = so_far + ((u, v, rs, enz_of[(u, v, rs)], tuple(copies_of[(u, v, rs)])),) + if v in states: + if len(new) > POOL_MAX_STEPS: + local["long_paths_dropped"] += 1 # a whole path, not a prefix + else: + paths.append(new) + if len(paths) + local["long_paths_dropped"] > POOL_PATH_BUDGET: + over_budget = True + return + continue + walk(v, new, seen | {v}) + for s0 in sorted(states): + walk(s0, (), {s0}) + if over_budget: + local["budget_dropped"] += 1 + continue + if not paths: + local["no_path_pools"] += 1 + continue + on_path: Set[str] = {n for p in paths for st in p for n in (st[0], st[1])} + kept_states = states & on_path + if len({sig(u, r) for u in kept_states}) < 2: + local["no_path_pools"] += 1 + continue + inter = (scc - states) & on_path + local["off_path_states"] += len(states - on_path) + local["off_path_intermediates"] += len(scc - states - on_path) + key = (frozenset(kept_states), frozenset(inter), + frozenset(tuple((a_u, b_u, rs, cp) for a_u, b_u, rs, _, cp in p) for p in paths)) + if key in pools: + pools[key]["r"].add(r) + local["merged_proteins"] += 1 + else: + pools[key] = {"r": {r}, "states": kept_states, "inter": inter, "paths": paths, + "nodes": kept_states | inter | {c for p in paths for st in p for c in st[4]}} + return list(pools.values()), local + + removed: Set[str] = set() + dropped_shared = 0 + while True: + pools, local = build(removed) + claimed: Dict[str, int] = {} + for p in pools: + for n in p["nodes"]: + claimed[n] = claimed.get(n, 0) + 1 + shared = {n for n, c in claimed.items() if c > 1} + if not shared: + break + removed |= shared + dropped_shared += len(shared) + stats.update(local) + managed: Set[str] = set().union(*(p["nodes"] for p in pools)) if pools else set() + stats["shared_nodes_dropped"] = dropped_shared + + def path_key(p): + return tuple((ent.get(a, ""), a, ent.get(b, ""), b, rs) for a, b, rs, _, _ in p) + + def modified_side(x, y, xy, yx): + """The state a pair marks as modified, or None when undecidable.""" + mx, my = prof(x, "mods"), prof(y, "mods") + if mx != my: + return y if my > mx else x + fd = any(rs in donor_reactions for p in xy for _, _, rs, _, _ in p) + rd = any(rs in donor_reactions for p in yx for _, _, rs, _, _ in p) + if fd != rd: + return y if fd else x + cx, cy = prof(x, "comps"), prof(y, "comps") + if cx != cy: + return y if cy > cx else x + return None + + forms: List[Tuple] = [] + transitions: List[Tuple] = [] + carriers: List[Tuple] = [] + pools.sort(key=lambda p: min(ent.get(u, "") + u for u in p["states"])) + for i, p in enumerate(pools, start=1): + pid = f"pool{i}" + states, inter, paths, rset = p["states"], p["inter"], sorted(p["paths"], key=path_key), p["r"] + between: Dict[Tuple[str, str], List[tuple]] = {} + for path in paths: + between.setdefault((path[0][0], path[-1][1]), []).append(path) + marked: Dict[str, int] = {} + for x in sorted(states): + for y in sorted(states): + if x >= y or not (between.get((x, y)) or between.get((y, x))): + continue + side = modified_side(x, y, between.get((x, y), []), between.get((y, x), [])) + if side is not None: + marked[side] = marked.get(side, 0) + 1 + base = min(states, key=lambda u: (marked.get(u, 0), prof(u, "mods"), prof(u, "comps"), ent.get(u, ""), u)) + if not marked: + stats["ties"] += 1 + for u in sorted(states | inter, key=lambda u: (ent.get(u, ""), u)): + forms.append((pid, u, ent.get(u, ""), "state" if u in states else "intermediate", u == base)) + for path_no, path in enumerate(paths, start=1): + for step_no, (a_u, b_u, rs, enz, cp) in enumerate(path, start=1): + for rx_u in sorted(cp): + transitions.append((pid, f"{pid}_p{path_no}", step_no, a_u, b_u, rx_u, rs, enz)) + rx_on_paths = sorted({c for path in paths for st in path for c in st[4]}) + found: Set[Tuple[str, str]] = set() + autocat: Set[Tuple[str, str, str]] = set() + for path in paths: + for a_u, b_u, _, enz, cp in path: + for rx_u in cp: + for c in cats.get(rx_u, ()): + if c in p["nodes"]: # a form of the pool catalyses a step of its own pool + autocat.add((rx_u, c, "autocat_source" if c == a_u else "autocat_product" if c == b_u + else "autocat_other")) + if not enz: + continue + for j in ins.get(rx_u, {}): + if j == a_u or j in p["nodes"] or pe_of(j) in _UBIQUITIN_STIDS: + continue + pj = proteins(j) + if not pj or pj & rset: + continue + for rel in rx_on_paths: + if j in outs.get(rel, {}): + found.add((j, rel)) + # a carrier that a pool manages (a state, an intermediate or a step node of + # any pool) is not a carrier: its value is that pool's to write + conflicts = {j for j, _ in found if j in managed} + stats["carrier_conflicts"] += len(conflicts) + found = {(j, rel) for j, rel in found if j not in conflicts} + carriers.extend((pid, j, rel) for j, rel in sorted(found, key=lambda t: (ent.get(t[0], ""), t[0], t[1]))) + for _, _, kind in autocat: + stats[kind] += 1 + stats["pools"] += 1 + stats["states"] += len(states) + stats["intermediates"] += len(inter) + stats["paths"] += len(paths) + stats["copy_rows"] += sum(len(cp) for path in paths for *_, cp in path) + stats["multi_step_pools"] += int(any(len(path) > 1 for path in paths)) + stats["carriers"] += len({j for j, _ in found}) + return forms, transitions, carriers + + +def r_steps(participants, profiles: Dict[str, Dict[str, Any]], ubiquitin_stids=()) -> List[Tuple[str, Any, str, str, bool]]: + """The R-steps of a pathway at stId level (deltasignal specs/039 amendment + 2), from :func:`get_reaction_participants` rows: for each reaction and each + reference entity R, the reaction is a step if exactly one input and exactly + one output contain R, both at stoichiometry 1, and they differ. A step is + enzyme-driven if a catalyst or a joining input (another input; ubiquitin + is a co-substrate) contains a protein and does not contain R: a catalyst + that is itself a form of R (RAS:GTP's intrinsic hydrolysis, curated with + the set p21 RAS:GTP as its own catalyst) is self-catalysis, which counts as + non-enzyme. Returns sorted ``(reaction, R, input, output, enzyme_driven)``.""" + by_rx: Dict[str, Dict[str, List[Tuple[str, int]]]] = {} + for rx, role, pe, st in participants: + by_rx.setdefault(rx, {}).setdefault(role, []).append((pe, st)) + out = [] + for rx in sorted(by_rx): + roles = by_rx[rx] + rs: Set[Any] = set() + for pe, _ in roles.get("input", []): + rs |= set(profiles.get(pe, {}).get("proteins", ())) + for r in sorted(rs): + a = [(pe, st) for pe, st in roles.get("input", []) if r in profiles.get(pe, {}).get("proteins", ())] + b = [(pe, st) for pe, st in roles.get("output", []) if r in profiles.get(pe, {}).get("proteins", ())] + if len(a) != 1 or len(b) != 1 or a[0][1] != 1 or b[0][1] != 1 or a[0][0] == b[0][0]: + continue + helpers = [pe for pe, _ in roles.get("catalyst", [])] + helpers += [pe for pe, _ in roles.get("input", []) if pe != a[0][0] and pe not in ubiquitin_stids] + enz = any(profiles.get(pe, {}).get("proteins") and r not in profiles[pe]["proteins"] for pe in helpers) + out.append((rx, r, a[0][0], b[0][0], enz)) + return out + + +def export_pools(pathway_id: str, pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFrame, + uuid_mapping: Dict[str, str], pools_file: str, transitions_file: str, + carriers_file: str) -> None: + """Write pools.csv, pool_transitions.csv and pool_carriers.csv (deltasignal + specs/039 amendment 2): each protein's modification-cycle pool, for a + consumer that solves it at steady state (pi = pi P) instead of multiplying + around the loop. All three files are written even when empty, so a consumer + can tell "no pools" from "no table".""" + from src.neo4j_connector import (get_reaction_participants, get_form_profiles, + get_donor_reactions, get_set_descendants) + participants = get_reaction_participants(pathway_id) + ent = _uuid_to_stable_id_map(pathway_logic_network, uuid_mapping) + wanted = {pe for _, _, pe, _ in participants} + for s in ent.values(): + if "::variant::" in s: + parent, members = s.split("::variant::", 1) + wanted.add(parent) + wanted.update(m for m in members.split("_") if m) + else: + wanted.add(s) + profiles = get_form_profiles(wanted) + steps = r_steps(participants, profiles, _UBIQUITIN_STIDS) + donors = get_donor_reactions({s[0] for s in steps}, _UBIQUITIN_STIDS) + set_members = get_set_descendants({x for s in steps for x in (s[2], s[3])}) + stats: Dict[str, int] = {} + forms, transitions, carriers = find_pools(pathway_logic_network, reaction_id_map, uuid_mapping, steps, + profiles, donors, set_members, stats) + pd.DataFrame(forms, columns=["pool_id", "node_uuid", "stable_id", "role", "is_base"]).to_csv(pools_file, index=False) + pd.DataFrame(transitions, columns=["pool_id", "path_id", "step", "source_uuid", "target_uuid", "reaction_uuid", + "reaction_stid", "enzyme_driven"]).to_csv(transitions_file, index=False) + pd.DataFrame(carriers, columns=["pool_id", "carrier_uuid", "release_reaction_uuid"]).to_csv(carriers_file, index=False) + logger.info( + f"Exported {stats['pools']} pools ({stats['states']} states, {stats['intermediates']} intermediates, " + f"{stats['paths']} paths in {stats['copy_rows']} step-copy rows, {stats['multi_step_pools']} multi-step; " + f"{len(steps)} R-steps at stId level; " + f"{stats['long_paths_dropped']} paths over {POOL_MAX_STEPS} steps dropped, " + f"{stats['shared_nodes_dropped']} shared nodes dropped, {stats['ties']} oriented by tie-break, " + f"{stats['carriers']} carriers ({stats['carrier_conflicts']} excluded as pool nodes), " + f"{stats['carrier_loops']} carrier loops, " + f"{stats['merged_proteins']} co-travelling proteins merged, {stats['no_path_pools']} candidates with no path, " + f"{stats['off_path_states']}/{stats['off_path_intermediates']} states/intermediates off every path, " + f"{stats['ambiguous_copies']} ambiguous reaction copies, {stats['budget_dropped']} over the path budget; " + f"autocatalysis by source form {stats['autocat_source']}, by product form {stats['autocat_product']}, " + f"by another form {stats['autocat_other']}) to {pools_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 cdc3b45..2873847 100755 --- a/src/neo4j_connector.py +++ b/src/neo4j_connector.py @@ -932,6 +932,170 @@ def derived(s: str, depth: int = 0) -> bool: for s in ids if derived(s)} +# Group donors: a reaction consuming one ADDS a group to its substrate, so its +# product is the modified form of an interconversion pool (specs/039). ChEBI +# ids: ATP, GTP, S-adenosyl-L-methionine, acetyl-CoA, NAD+. +_DONOR_CHEBI = ("30616", "37565", "15414", "59789", "15351", "57288", "15846", "57540") + + +def get_donor_reactions(reaction_ids, ubiquitin_stids=()) -> Set[str]: + """Reactions (stIds) that consume a group donor or ubiquitin.""" + ids = sorted({r for r in reaction_ids if r}) + if not ids: + return set() + query = """ + UNWIND $ids AS rid MATCH (r:ReactionLikeEvent {stId: rid})-[:input]->(i:PhysicalEntity) + OPTIONAL MATCH (i)-[:referenceEntity]->(re:ReferenceMolecule) + WITH rid, i, re WHERE re.identifier IN $chebi OR i.stId IN $ub + RETURN DISTINCT rid + """ + try: + rows = get_graph().run(query, ids=ids, chebi=list(_DONOR_CHEBI), ub=list(ubiquitin_stids)).data() + except Exception: + logger.error("Error in get_donor_reactions", **_traceback_kwargs()) + raise + return {r["rid"] for r in rows} + + +def get_reaction_participants(pathway_id: str) -> List[Tuple[str, str, str, int]]: + """(reaction stId, role, entity stId, stoichiometry) for every input, output + and catalyst of every reaction-like event of a pathway (deltasignal + specs/039 amendment 2). Roles are ``input``, ``output`` and ``catalyst``; + a catalyst's stoichiometry is 1. One round trip per role.""" + io_query = """ + MATCH (p:Pathway {stId: $pid})-[:hasEvent*]->(rx:ReactionLikeEvent) + WITH DISTINCT rx + MATCH (rx)-[e:input|output]->(pe:PhysicalEntity) + RETURN rx.stId AS rx, type(e) AS role, pe.stId AS pe, coalesce(e.stoichiometry, 1) AS st + """ + cat_query = """ + MATCH (p:Pathway {stId: $pid})-[:hasEvent*]->(rx:ReactionLikeEvent) + WITH DISTINCT rx + MATCH (rx)-[:catalystActivity]->(:CatalystActivity)-[:physicalEntity]->(pe:PhysicalEntity) + RETURN DISTINCT rx.stId AS rx, 'catalyst' AS role, pe.stId AS pe, 1 AS st + """ + try: + rows = get_graph().run(io_query, pid=pathway_id).data() + rows += get_graph().run(cat_query, pid=pathway_id).data() + except Exception: + logger.error("Error in get_reaction_participants", **_traceback_kwargs()) + raise + return sorted({(r["rx"], r["role"], r["pe"], int(r["st"] or 1)) for r in rows}) + + +_form_profile_cache: Dict[str, Dict[str, Any]] = {} + + +def get_form_profiles(stable_ids) -> Dict[str, Dict[str, Any]]: + """stId -> what a form is, for pool detection (deltasignal specs/039 + amendment 2). Per entity: + + - ``proteins``: reference entities (dbIds) of its EWAS leaves, through + components, members and candidates; + - ``fixed``: the subset reached through components only (no set on the + way), which every variant of the entity carries; + - ``sig``: per reference entity R, its modification signature: the + modified residues on R's leaves plus the small molecules (ChEBI id, else + name) that are direct components of the complex R's leaf sits in, looking + through sets (GTP attaches to RAS in p21 RAS:GTP = [GTP, {H,K,N}RAS], and + not to the GAP in p21 RAS:GTP:RAS GAPs, which nests it); + - ``mods``: modified residues summed over all leaves (amendment 1's base + rule), ``comps``: components at any depth, ``slots``: non-small-molecule + positions (a set inside a complex is one position; a bare protein or set + is one). + + Cached per process: the 92 catalog pathways share most entities.""" + ids = sorted({s for s in stable_ids if s}) + missing = [s for s in ids if s not in _form_profile_cache] + leaves_query = """ + UNWIND $ids AS s MATCH (e:PhysicalEntity {stId: s}) + OPTIONAL MATCH p = (e)-[:hasComponent|hasMember|hasCandidate*0..8]->(l:EntityWithAccessionedSequence)-[:referenceEntity]->(re) + RETURN s, re.dbId AS r, + CASE WHEN p IS NULL THEN [] ELSE + [n IN nodes(p) WHERE n:PhysicalEntity | + [n.stId, CASE WHEN n:Complex THEN 'C' WHEN n:EntitySet THEN 'S' ELSE 'L' END]] END AS chain, + CASE WHEN l IS NULL THEN [] ELSE [(l)-[:hasModifiedResidue]->(m) | m.displayName] END AS res + """ + counts_query = """ + UNWIND $ids AS s MATCH (e:PhysicalEntity {stId: s}) + RETURN s, + size([(e)-[:hasComponent*1..6]->(c) | c]) AS comps, + size([(e)-[:hasComponent*0..6]->(x) WHERE NOT x:Complex AND NOT x:SimpleEntity | x]) AS slots + """ + smalls_query = """ + UNWIND $ids AS s MATCH (c:Complex {stId: s})-[:hasComponent]->(sm:SimpleEntity) + OPTIONAL MATCH (sm)-[:referenceEntity]->(rm) + RETURN s, coalesce(rm.identifier, head(sm.name), sm.displayName) AS mol + """ + try: + for i in range(0, len(missing), 200): + chunk = missing[i:i + 200] + prof: Dict[str, Dict[str, Any]] = { + s: {"proteins": set(), "fixed": set(), "sig": {}, "mods": 0, "comps": 0, "slots": 0, + "_res": {}, "_cx": {}} for s in chunk} + leaf_rows = get_graph().run(leaves_query, ids=chunk).data() + for r in get_graph().run(counts_query, ids=chunk).data(): + prof[r["s"]]["comps"] = int(r["comps"]) + prof[r["s"]]["slots"] = int(r["slots"]) + complexes = set() + for r in leaf_rows: + if r["r"] is None: + continue + pr, rid, chain, res = prof[r["s"]], int(r["r"]), r["chain"], list(r["res"]) + pr["proteins"].add(rid) + if all(k != "S" for _, k in chain[:-1]): + pr["fixed"].add(rid) + pr["mods"] += len(res) + pr["_res"].setdefault(rid, set()).update(res) + # the complex R's leaf sits in: the nearest Complex above it with + # only set steps between + cx = None + for st_id, kind in reversed(chain[:-1]): + if kind == "C": + cx = st_id + break + if kind != "S": + break + if cx is not None: + pr["_cx"].setdefault(rid, set()).add(cx) + complexes.add(cx) + smalls: Dict[str, Set[str]] = {} + for j in range(0, len(complexes), 200): + for r in get_graph().run(smalls_query, ids=sorted(complexes)[j:j + 200]).data(): + smalls.setdefault(r["s"], set()).add(str(r["mol"])) + for s, pr in prof.items(): + for rid in pr["proteins"]: + sig = {("res", x) for x in pr["_res"].get(rid, ())} + for cx in pr["_cx"].get(rid, ()): + sig |= {("mol", x) for x in smalls.get(cx, ())} + pr["sig"][rid] = frozenset(sig) + del pr["_res"], pr["_cx"] + _form_profile_cache[s] = pr + except Exception: + logger.error("Error in get_form_profiles", **_traceback_kwargs()) + raise + return {s: _form_profile_cache[s] for s in ids if s in _form_profile_cache} + + +def get_set_descendants(stable_ids) -> Dict[str, Set[str]]: + """stId -> the stIds an EntitySet reaches through members and candidates + (deltasignal specs/039 amendment 2): a set-valued participant appears in + the network as its member nodes, and this says which nodes those are.""" + ids = sorted({s for s in stable_ids if s}) + if not ids: + return {} + query = """ + UNWIND $ids AS s MATCH (e:EntitySet {stId: s}) + RETURN s, [(e)-[:hasMember|hasCandidate*1..6]->(m) | m.stId] AS members + """ + try: + rows = get_graph().run(query, ids=ids).data() + except Exception: + logger.error("Error in get_set_descendants", **_traceback_kwargs()) + raise + return {r["s"]: {str(m) for m in r["members"] if m} for r in rows} + + def get_reactome_release() -> Optional[int]: """The Reactome release number of the connected graph, or None. diff --git a/src/pathway_generator.py b/src/pathway_generator.py index 8b9d873..897b6bf 100755 --- a/src/pathway_generator.py +++ b/src/pathway_generator.py @@ -12,6 +12,7 @@ create_pathway_logic_network, export_cofactors, export_drugs, + export_pools, export_containment, export_entity_reaction_proxy_mapping, export_node_reaction_context, @@ -510,6 +511,20 @@ def generate_pathway_file( except Exception as e: logger.error(f"Failed to write drugs.csv: {e}", exc_info=True) + # Interconversion pools (deltasignal specs/039), in their own block. + try: + export_pools( + pathway_id, + result.logic_network, + result.reaction_id_map, + result.uuid_mapping, + str(pathway_output_dir / "pools.csv"), + str(pathway_output_dir / "pool_transitions.csv"), + str(pathway_output_dir / "pool_carriers.csv"), + ) + except Exception as e: + logger.error(f"Failed to write pools.csv: {e}", exc_info=True) + logger.info(f"Output directory: {pathway_output_dir}") except (ConnectionError, ValueError) as e: diff --git a/tests/test_pools.py b/tests/test_pools.py new file mode 100644 index 0000000..bd1914f --- /dev/null +++ b/tests/test_pools.py @@ -0,0 +1,387 @@ +"""Modification-cycle pools (deltasignal specs/039 amendment 2): a protein's +forms joined by R-steps, found at NODE level; states are the least-bound form +of each modification signature, the enzyme complexes on the way are +intermediates, and a transition is a path state -> intermediates* -> state. +Shapes: RAS (GEF / intrinsic / GAP bind-hydrolyse-release), a six-form +kinase/phosphatase ring, an enzyme's own binding loop.""" +import pandas as pd + +import src.logic_network_generator as m + +# _uuid_to_stable_id_map detects the mapping direction by uuid SHAPE +NAMES = ["g", "t", "tx", "gap", "gef", "vgef", "vhyd", "vbind", "vrel", "a2", + "s", "se", "sxe", "sx", "sxp", "sp", "e", "p", "v1", "v2", "v3", "v4", "v5", "v6", + "x", "y", "xy", "va", "vb", "i1", "i2", "i3", "i4", "i5", "i6", "v7", "v8", "vback", + "xsy", "xs", "ys", "vc", "vd", "ve"] +U = {k: f"00000000-0000-0000-0000-{i:012d}" for i, k in enumerate(NAMES, 1)} +RAS, GAP_P, GEF_P, SUB, KIN, PHOS = 1, 2, 3, 10, 11, 12 # reference entities + + +def e(s, t, et, st=1): + return {"source_id": U[s], "target_id": U[t], "edge_type": et, "stoichiometry": st} + + +def prof(proteins, sig=None, slots=1, mods=0, comps=0, fixed=None): + return {"proteins": set(proteins), "fixed": set(fixed if fixed is not None else proteins), + "sig": {r: frozenset(sig or ()) for r in proteins}, "mods": mods, "comps": comps, "slots": slots} + + +def network(edges, reactions, umap): + return (pd.DataFrame(edges), + pd.DataFrame({"uid": [U[k] for k in reactions], "reactome_id": [f"R-HSA-{k}" for k in reactions]}), + {U[k]: s for k, s in umap.items()}) + + +# --- RAS: GDP <-> GTP by GEF (enzyme), intrinsic exchange / hydrolysis (self), GAP bind + release ---------- + +def ras_network(hydrolysis_catalyst="t"): + edges = [e("g", "vgef", "input"), e("vgef", "t", "output"), e("gef", "vgef", "catalyst"), + e("t", "vhyd", "input"), e("vhyd", "g", "output"), e(hydrolysis_catalyst, "vhyd", "catalyst"), + e("t", "vbind", "input"), e("gap", "vbind", "input"), e("vbind", "tx", "output"), + e("tx", "vrel", "input"), e("vrel", "g", "output"), e("vrel", "gap", "output")] + net, rmap, umap = network(edges, ["vgef", "vhyd", "vbind", "vrel"], + {"g": "R-HSA-G", "t": "R-HSA-T", "tx": "R-HSA-TX", "gap": "R-HSA-GAP", "gef": "R-HSA-GEF"}) + profiles = {"R-HSA-G": prof([RAS], [("mol", "GDP")], slots=1, comps=2), + "R-HSA-T": prof([RAS], [("mol", "GTP")], slots=1, comps=2), + "R-HSA-TX": {"proteins": {RAS, GAP_P}, "fixed": {RAS, GAP_P}, + "sig": {RAS: frozenset([("mol", "GTP")]), GAP_P: frozenset()}, "mods": 0, "comps": 3, "slots": 2}, + "R-HSA-GAP": prof([GAP_P]), "R-HSA-GEF": prof([GEF_P])} + participants = [("R-HSA-vgef", "input", "R-HSA-G", 1), ("R-HSA-vgef", "output", "R-HSA-T", 1), + ("R-HSA-vgef", "catalyst", "R-HSA-GEF", 1), + ("R-HSA-vhyd", "input", "R-HSA-T", 1), ("R-HSA-vhyd", "output", "R-HSA-G", 1), + ("R-HSA-vhyd", "catalyst", "R-HSA-T", 1), + ("R-HSA-vbind", "input", "R-HSA-T", 1), ("R-HSA-vbind", "input", "R-HSA-GAP", 1), + ("R-HSA-vbind", "output", "R-HSA-TX", 1), + ("R-HSA-vrel", "input", "R-HSA-TX", 1), ("R-HSA-vrel", "output", "R-HSA-G", 1), + ("R-HSA-vrel", "output", "R-HSA-GAP", 1)] + return net, rmap, umap, profiles, participants + + +def test_r_steps_flag_self_catalysis_as_non_enzyme(): + _, _, _, profiles, participants = ras_network() + steps = m.r_steps(participants, profiles) + by = {(rx, r): (a, b, enz) for rx, r, a, b, enz in steps} + assert by[("R-HSA-vgef", RAS)] == ("R-HSA-G", "R-HSA-T", True) # GEF catalyst: enzyme + assert by[("R-HSA-vhyd", RAS)] == ("R-HSA-T", "R-HSA-G", False) # catalysed by RAS:GTP itself + assert by[("R-HSA-vbind", RAS)] == ("R-HSA-T", "R-HSA-TX", True) # joining input: the GAP + assert by[("R-HSA-vrel", RAS)] == ("R-HSA-TX", "R-HSA-G", False) # release: nothing joins + assert by[("R-HSA-vbind", GAP_P)] == ("R-HSA-GAP", "R-HSA-TX", True) # for the GAP, RAS:GTP joins + assert ("R-HSA-vrel", GAP_P) in by + + +def test_ras_pool_has_gdp_base_gtp_state_and_the_gap_complex_as_intermediate(): + net, rmap, umap, profiles, participants = ras_network() + st = {} + forms, trans, carriers = m.find_pools(net, rmap, umap, m.r_steps(participants, profiles), profiles, + {"R-HSA-vgef"}, {}, st) + assert [(u, role, base) for _, u, _, role, base in forms] == [ + (U["g"], "state", True), (U["t"], "state", False), (U["tx"], "intermediate", False)] + paths = {} + for _, pid, step, a, b, rx, rs, enz in trans: + paths.setdefault(pid, []).append((step, a, b, rx, rs, enz)) + shapes = sorted(tuple((a, b, rs, enz) for _, a, b, _, rs, enz in sorted(p)) for p in paths.values()) + assert shapes == sorted([ + ((U["g"], U["t"], "R-HSA-vgef", True),), + ((U["t"], U["g"], "R-HSA-vhyd", False),), + ((U["t"], U["tx"], "R-HSA-vbind", True), (U["tx"], U["g"], "R-HSA-vrel", False)), + ]) + assert carriers == [("pool1", U["gap"], U["vrel"])] + assert st["pools"] == 1 and st["states"] == 2 and st["intermediates"] == 1 and st["paths"] == 3 + assert st["multi_step_pools"] == 1 and st["carriers"] == 1 and st["ties"] == 0 + assert st["carrier_loops"] == 1 # the GAP's own GAP -> RAS:GTP:GAP -> GAP loop + assert st["autocat_source"] == 1 and st["autocat_product"] == 0 and st["autocat_other"] == 0 + + +def test_a_step_whose_node_loop_is_broken_is_not_a_transition(): + net, rmap, umap, profiles, participants = ras_network() + # the release reaction returns a DIFFERENT copy of RAS:GDP + net.loc[(net["source_id"] == U["vrel"]) & (net["target_id"] == U["g"]), "target_id"] = U["a2"] + umap[U["a2"]] = "R-HSA-G" + st = {} + forms, trans, _ = m.find_pools(net, rmap, umap, m.r_steps(participants, profiles), profiles, set(), {}, st) + assert {u for _, u, _, _, _ in forms} == {U["g"], U["t"]} and st["paths"] == 2 and st["intermediates"] == 0 + + +def test_a_form_node_that_is_a_member_of_the_steps_set_maps_to_it(): + net, rmap, umap, profiles, participants = ras_network() + participants = [(rx, role, "R-HSA-GSET" if pe == "R-HSA-G" else pe, s) for rx, role, pe, s in participants] + profiles["R-HSA-GSET"] = profiles["R-HSA-G"] + steps = m.r_steps(participants, profiles) + assert m.find_pools(net, rmap, umap, steps, profiles, set(), {}, {})[0] == [] + forms, _, _ = m.find_pools(net, rmap, umap, steps, profiles, set(), {"R-HSA-GSET": {"R-HSA-G"}}, {}) + assert len(forms) == 3 + + +# --- six-form ring S -> S:E -> S*:E -> S* -> S*:P -> S:P -> S ----------------------------------------------- + +def ring_network(): + edges = [e("s", "v1", "input"), e("e", "v1", "input"), e("v1", "se", "output"), + e("se", "v2", "input"), e("v2", "sxe", "output"), + e("sxe", "v3", "input"), e("v3", "sx", "output"), e("v3", "e", "output"), + e("sx", "v4", "input"), e("p", "v4", "input"), e("v4", "sxp", "output"), + e("sxp", "v5", "input"), e("v5", "sp", "output"), + e("sp", "v6", "input"), e("v6", "s", "output"), e("v6", "p", "output")] + net, rmap, umap = network(edges, ["v1", "v2", "v3", "v4", "v5", "v6"], + {"s": "R-HSA-S", "se": "R-HSA-SE", "sxe": "R-HSA-SXE", "sx": "R-HSA-SX", + "sxp": "R-HSA-SXP", "sp": "R-HSA-SP", "e": "R-HSA-E", "p": "R-HSA-P"}) + ph = [("res", "phospho-S at 1")] + profiles = {"R-HSA-S": prof([SUB]), "R-HSA-SX": prof([SUB], ph, mods=1), + "R-HSA-E": prof([KIN]), "R-HSA-P": prof([PHOS]), + "R-HSA-SE": {"proteins": {SUB, KIN}, "fixed": {SUB, KIN}, "sig": {SUB: frozenset(), KIN: frozenset()}, + "mods": 0, "comps": 2, "slots": 2}, + "R-HSA-SXE": {"proteins": {SUB, KIN}, "fixed": {SUB, KIN}, "sig": {SUB: frozenset(ph), KIN: frozenset()}, + "mods": 1, "comps": 2, "slots": 2}, + "R-HSA-SXP": {"proteins": {SUB, PHOS}, "fixed": {SUB, PHOS}, "sig": {SUB: frozenset(ph), PHOS: frozenset()}, + "mods": 1, "comps": 2, "slots": 2}, + "R-HSA-SP": {"proteins": {SUB, PHOS}, "fixed": {SUB, PHOS}, "sig": {SUB: frozenset(), PHOS: frozenset()}, + "mods": 0, "comps": 2, "slots": 2}} + participants = [("R-HSA-v1", "input", "R-HSA-S", 1), ("R-HSA-v1", "input", "R-HSA-E", 1), ("R-HSA-v1", "output", "R-HSA-SE", 1), + ("R-HSA-v2", "input", "R-HSA-SE", 1), ("R-HSA-v2", "output", "R-HSA-SXE", 1), + ("R-HSA-v3", "input", "R-HSA-SXE", 1), ("R-HSA-v3", "output", "R-HSA-SX", 1), ("R-HSA-v3", "output", "R-HSA-E", 1), + ("R-HSA-v4", "input", "R-HSA-SX", 1), ("R-HSA-v4", "input", "R-HSA-P", 1), ("R-HSA-v4", "output", "R-HSA-SXP", 1), + ("R-HSA-v5", "input", "R-HSA-SXP", 1), ("R-HSA-v5", "output", "R-HSA-SP", 1), + ("R-HSA-v6", "input", "R-HSA-SP", 1), ("R-HSA-v6", "output", "R-HSA-S", 1), ("R-HSA-v6", "output", "R-HSA-P", 1)] + return net, rmap, umap, profiles, participants + + +def test_six_form_ring_gives_two_states_four_intermediates_two_paths_two_carriers(): + net, rmap, umap, profiles, participants = ring_network() + st = {} + forms, trans, carriers = m.find_pools(net, rmap, umap, m.r_steps(participants, profiles), profiles, set(), {}, st) + roles = {u: (role, base) for _, u, _, role, base in forms} + assert roles == {U["s"]: ("state", True), U["sx"]: ("state", False), + U["se"]: ("intermediate", False), U["sxe"]: ("intermediate", False), + U["sxp"]: ("intermediate", False), U["sp"]: ("intermediate", False)} + paths = {} + for _, pid, step, a, b, rx, rs, enz in trans: + paths.setdefault(pid, []).append((step, a, b, rx, enz)) + assert len(paths) == 2 and st["paths"] == 2 and st["multi_step_pools"] == 1 + for p in paths.values(): + assert [s for s, *_ in sorted(p)] == [1, 2, 3] + assert any(enz for *_, enz in p) # the binding step joins the enzyme + assert [enz for *_, enz in sorted(p)] == [True, False, False] + assert {(p[0][1], p[-1][2]) for p in (sorted(v) for v in paths.values())} == {(U["s"], U["sx"]), (U["sx"], U["s"])} + assert sorted(carriers) == sorted([("pool1", U["e"], U["v3"]), ("pool1", U["p"], U["v6"])]) + assert st["carriers"] == 2 and st["carrier_loops"] == 2 and st["ties"] == 0 + + +def test_an_enzymes_own_binding_loop_is_not_a_pool(): + # E -> E:S -> E: one signature for E, so a carrier loop, not a pool + edges = [e("e", "v1", "input"), e("v1", "se", "output"), e("se", "v3", "input"), e("v3", "e", "output")] + net, rmap, umap = network(edges, ["v1", "v3"], {"e": "R-HSA-E", "se": "R-HSA-SE"}) + profiles = {"R-HSA-E": prof([KIN]), "R-HSA-SE": prof([KIN, SUB], slots=2, comps=2)} + steps = [("R-HSA-v1", KIN, "R-HSA-E", "R-HSA-SE", True), ("R-HSA-v3", KIN, "R-HSA-SE", "R-HSA-E", False)] + st = {} + assert m.find_pools(net, rmap, umap, steps, profiles, set(), {}, st) == ([], [], []) + assert st["carrier_loops"] == 1 and st["pools"] == 0 + + +# --- shared nodes, path cap, orientation -------------------------------------------------------------------- + +def test_a_node_claimed_by_two_pools_is_dropped_and_counted(): + # X + Y -> X:Y -> X*:Y* -> X* + Y*; X* -> X and Y* -> Y directly. X's pool and + # Y's pool both claim the two complexes and the three reaction nodes. + edges = [e("x", "va", "input"), e("y", "va", "input"), e("va", "xy", "output"), + e("xy", "vb", "input"), e("vb", "xsy", "output"), + e("xsy", "vc", "input"), e("vc", "xs", "output"), e("vc", "ys", "output"), + e("xs", "vd", "input"), e("vd", "x", "output"), e("ys", "ve", "input"), e("ve", "y", "output")] + net, rmap, umap = network(edges, ["va", "vb", "vc", "vd", "ve"], + {"x": "R-HSA-X", "y": "R-HSA-Y", "xy": "R-HSA-XY", "xsy": "R-HSA-XSY", + "xs": "R-HSA-XS", "ys": "R-HSA-YS"}) + p1, p2 = [("res", "p1")], [("res", "p2")] + profiles = {"R-HSA-X": prof([1]), "R-HSA-Y": prof([2]), "R-HSA-XS": prof([1], p1, mods=1), "R-HSA-YS": prof([2], p2, mods=1), + "R-HSA-XY": {"proteins": {1, 2}, "fixed": {1, 2}, "mods": 0, "comps": 2, "slots": 2, + "sig": {1: frozenset(), 2: frozenset()}}, + "R-HSA-XSY": {"proteins": {1, 2}, "fixed": {1, 2}, "mods": 2, "comps": 2, "slots": 2, + "sig": {1: frozenset(p1), 2: frozenset(p2)}}} + steps = [("R-HSA-va", 1, "R-HSA-X", "R-HSA-XY", True), ("R-HSA-vb", 1, "R-HSA-XY", "R-HSA-XSY", False), + ("R-HSA-vc", 1, "R-HSA-XSY", "R-HSA-XS", False), ("R-HSA-vd", 1, "R-HSA-XS", "R-HSA-X", True), + ("R-HSA-va", 2, "R-HSA-Y", "R-HSA-XY", True), ("R-HSA-vb", 2, "R-HSA-XY", "R-HSA-XSY", False), + ("R-HSA-vc", 2, "R-HSA-XSY", "R-HSA-YS", False), ("R-HSA-ve", 2, "R-HSA-YS", "R-HSA-Y", True)] + # each alone is a pool + for keep in (1, 2): + st = {} + m.find_pools(net, rmap, umap, [t for t in steps if t[1] == keep], profiles, set(), {}, st) + assert st["pools"] == 1 and st["intermediates"] == 2 + st = {} + assert m.find_pools(net, rmap, umap, steps, profiles, set(), {}, st) == ([], [], []) + assert st["shared_nodes_dropped"] == 5 and st["pools"] == 0 # X:Y, X*:Y* and the three reaction nodes + + +def test_forms_differing_only_by_a_bound_partner_are_not_a_pool(): + # Activated FGFR4 <-> FGFR4:PLCG1 <-> FGFR4:p-PLCG1 (R-HSA-5654743): the PLCG1 + # phosphorylation gives FGFR4 no new signature; the only core-only form is + # FGFR4 itself, so this is a carrier loop, not a pool (amendment 3) + edges = [e("x", "va", "input"), e("y", "va", "input"), e("va", "xy", "output"), + e("xy", "vb", "input"), e("vb", "xsy", "output"), + e("xsy", "vc", "input"), e("vc", "x", "output"), e("vc", "ys", "output")] + net, rmap, umap = network(edges, ["va", "vb", "vc"], + {"x": "R-HSA-X", "y": "R-HSA-Y", "xy": "R-HSA-XY", "xsy": "R-HSA-XSY", "ys": "R-HSA-YS"}) + p2 = [("res", "p2")] + profiles = {"R-HSA-X": prof([1]), "R-HSA-Y": prof([2]), "R-HSA-YS": prof([2], p2, mods=1), + "R-HSA-XY": {"proteins": {1, 2}, "fixed": {1, 2}, "mods": 0, "comps": 2, "slots": 2, + "sig": {1: frozenset(), 2: frozenset()}}, + "R-HSA-XSY": {"proteins": {1, 2}, "fixed": {1, 2}, "mods": 1, "comps": 2, "slots": 2, + "sig": {1: frozenset(), 2: frozenset(p2)}}} + steps = [("R-HSA-va", 1, "R-HSA-X", "R-HSA-XY", True), ("R-HSA-vb", 1, "R-HSA-XY", "R-HSA-XSY", False), + ("R-HSA-vc", 1, "R-HSA-XSY", "R-HSA-X", False)] + st = {} + assert m.find_pools(net, rmap, umap, steps, profiles, set(), {}, st) == ([], [], []) + assert st["carrier_loops"] == 1 + # and the same shape where the RECEPTOR is what gets phosphorylated IS a pool + profiles["R-HSA-XSY"]["sig"] = {1: frozenset(p2), 2: frozenset()} + profiles["R-HSA-XS"] = prof([1], p2, mods=1) + umap[U["ys"]] = "R-HSA-XS" + steps = [("R-HSA-va", 1, "R-HSA-X", "R-HSA-XY", True), ("R-HSA-vb", 1, "R-HSA-XY", "R-HSA-XSY", False), + ("R-HSA-vc", 1, "R-HSA-XSY", "R-HSA-XS", False)] + edges += [e("ys", "vd", "input"), e("vd", "x", "output")] + net, rmap, _ = network(edges, ["va", "vb", "vc", "vd"], {}) + steps.append(("R-HSA-vd", 1, "R-HSA-XS", "R-HSA-X", True)) + st = {} + forms, _, _ = m.find_pools(net, rmap, umap, steps, profiles, set(), {}, st) + assert st["pools"] == 1 and {u: r for _, u, _, r, _ in forms} == { + U["x"]: "state", U["ys"]: "state", U["xy"]: "intermediate", U["xsy"]: "intermediate"} + + +def test_a_carrier_that_another_pool_manages_is_excluded_and_counted(): + # the ring's kinase E is also a form of its own pool (E <-> E* by an outside + # kinase and phosphatase): E is not a carrier of the substrate's pool + net, rmap, umap, profiles, participants = ring_network() + edges = [e("e", "v7", "input"), e("v7", "i1", "output"), e("i1", "v8", "input"), e("v8", "e", "output")] + net = pd.concat([net, pd.DataFrame(edges)], ignore_index=True) + rmap = pd.concat([rmap, pd.DataFrame({"uid": [U["v7"], U["v8"]], "reactome_id": ["R-HSA-v7", "R-HSA-v8"]})], + ignore_index=True) + umap[U["i1"]] = "R-HSA-ES" + profiles["R-HSA-ES"] = prof([KIN], [("res", "pE")], mods=1) + participants += [("R-HSA-v7", "input", "R-HSA-E", 1), ("R-HSA-v7", "output", "R-HSA-ES", 1), + ("R-HSA-v7", "catalyst", "R-HSA-P", 1), + ("R-HSA-v8", "input", "R-HSA-ES", 1), ("R-HSA-v8", "output", "R-HSA-E", 1), + ("R-HSA-v8", "catalyst", "R-HSA-P", 1)] + st = {} + forms, _, carriers = m.find_pools(net, rmap, umap, m.r_steps(participants, profiles), profiles, set(), {}, st) + assert st["pools"] == 2 and st["shared_nodes_dropped"] == 0 + assert carriers == [(next(p for p, u, *_ in forms if u == U["s"]), U["p"], U["v6"])] + assert st["carrier_conflicts"] == 1 and st["carriers"] == 1 + + +def test_co_travelling_proteins_of_one_set_are_one_pool(): + # p21 RAS:GDP <-> RAS:GTP: the H/K/NRAS graphs are identical, so one pool, counted as merged + net, rmap, umap, profiles, participants = ras_network() + for s in ("R-HSA-G", "R-HSA-T"): + profiles[s] = prof([RAS, 4], profiles[s]["sig"][RAS], comps=2) + profiles["R-HSA-TX"]["proteins"].add(4) + profiles["R-HSA-TX"]["sig"][4] = profiles["R-HSA-TX"]["sig"][RAS] + st = {} + forms, _, _ = m.find_pools(net, rmap, umap, m.r_steps(participants, profiles), profiles, set(), {}, st) + assert st["pools"] == 1 and st["merged_proteins"] == 1 and st["shared_nodes_dropped"] == 0 and len(forms) == 3 + + +def chain_network(n_intermediates): + inter = ["i1", "i2", "i3", "i4", "i5", "i6"][:n_intermediates] + seq = ["s"] + inter + ["sx"] + rxs = ["v1", "v2", "v3", "v4", "v5", "v6", "v7"][:len(seq) - 1] + edges, umap, steps = [], {"s": "R-HSA-S", "sx": "R-HSA-SX"}, [] + profiles = {"R-HSA-S": prof([SUB]), "R-HSA-SX": prof([SUB], [("res", "p")], mods=1)} + for k, rx in enumerate(rxs): + a, b = seq[k], seq[k + 1] + edges += [e(a, rx, "input"), e(rx, b, "output")] + for x in (a, b): + if x in inter: + umap[x] = f"R-HSA-{x.upper()}" + profiles[f"R-HSA-{x.upper()}"] = prof([SUB], slots=2, comps=2) + steps.append((f"R-HSA-{rx}", SUB, umap[a], umap[b], True)) + edges += [e("sx", "vback", "input"), e("vback", "s", "output")] + steps.append(("R-HSA-vback", SUB, "R-HSA-SX", "R-HSA-S", True)) + net, rmap, umap = network(edges, rxs + ["vback"], umap) + return net, rmap, umap, steps, profiles + + +def test_paths_over_six_steps_are_dropped_and_counted(): + net, rmap, umap, steps, profiles = chain_network(5) # 6 steps: kept + st = {} + m.find_pools(net, rmap, umap, steps, profiles, set(), {}, st) + assert st["paths"] == 2 and st["long_paths_dropped"] == 0 and st["intermediates"] == 5 + net, rmap, umap, steps, profiles = chain_network(6) # 7 steps: dropped + st = {} + forms, _, _ = m.find_pools(net, rmap, umap, steps, profiles, set(), {}, st) + assert st["paths"] == 1 and st["long_paths_dropped"] == 1 + assert st["off_path_intermediates"] == 6 and [r for _, _, _, r, _ in forms] == ["state", "state"] + # a second reaction for one of its steps is a second dropped PATH, not a prefix + net = pd.concat([net, pd.DataFrame([e("i3", "v8", "input"), e("v8", "i4", "output")])], ignore_index=True) + rmap = pd.concat([rmap, pd.DataFrame({"uid": [U["v8"]], "reactome_id": ["R-HSA-v8"]})], ignore_index=True) + steps.append(("R-HSA-v8", SUB, "R-HSA-I3", "R-HSA-I4", True)) + st = {} + m.find_pools(net, rmap, umap, steps, profiles, set(), {}, st) + assert st["paths"] == 1 and st["long_paths_dropped"] == 2 + + +def test_orientation_by_residues_then_donor_then_components_and_ties_are_counted(): + net, rmap, umap, profiles, participants = ras_network() + steps = m.r_steps(participants, profiles) + + def base_of(profiles, donors, st=None): + forms, _, _ = m.find_pools(net, rmap, umap, steps, profiles, donors, {}, st if st is not None else {}) + return next(u for _, u, _, role, base in forms if base) + # no residue difference: the donor is consumed on the way to GTP, so GDP is base + assert base_of(profiles, {"R-HSA-vgef"}) == U["g"] + assert base_of(profiles, {"R-HSA-vhyd"}) == U["t"] + # residues decide first, whatever the donor says + profiles["R-HSA-G"]["mods"] = 2 + assert base_of(profiles, {"R-HSA-vgef"}) == U["t"] + profiles["R-HSA-G"]["mods"] = 0 + # no residue or donor difference: the form with more components is modified + profiles["R-HSA-T"]["comps"] = 3 + assert base_of(profiles, set()) == U["g"] + profiles["R-HSA-T"]["comps"] = 2 + # nothing decides: smaller stId, counted as a tie + st = {} + assert base_of(profiles, set(), st) == U["g"] and st["ties"] == 1 + + +def test_no_steps_or_empty_network_give_no_pools(): + net, rmap, umap, profiles, participants = ras_network() + assert m.find_pools(net, rmap, umap, [], profiles, set()) == ([], [], []) + empty = pd.DataFrame(columns=["source_id", "target_id", "edge_type"]) + assert m.find_pools(empty, rmap, umap, m.r_steps(participants, profiles), profiles, set()) == ([], [], []) + + +def test_a_stoichiometry_two_input_is_not_a_step(): + _, _, _, profiles, participants = ras_network() + participants = [(rx, role, pe, 2 if (rx, role) == ("R-HSA-vgef", "input") else s) for rx, role, pe, s in participants] + assert ("R-HSA-vgef", RAS) not in {(rx, r) for rx, r, *_ in m.r_steps(participants, profiles)} + + +def test_copies_of_one_reaction_collapse_into_one_path_with_one_row_per_copy(): + # RAS with the GAP binding curated for two GAP variants (two reaction NODES of + # one stId) and both release copies, plus an intrinsic exchange beside the GEF: + # copies of one stId are one step (k rows); different reactions stay apart. + net, rmap, umap, profiles, participants = ras_network() + extra = [e("t", "v7", "input"), e("a2", "v7", "input"), e("v7", "tx", "output"), + e("tx", "v8", "input"), e("v8", "g", "output"), e("v8", "a2", "output"), + e("g", "vback", "input"), e("vback", "t", "output")] + net = pd.concat([net, pd.DataFrame(extra)], ignore_index=True) + rmap = pd.concat([rmap, pd.DataFrame({"uid": [U["v7"], U["v8"], U["vback"]], + "reactome_id": ["R-HSA-vbind", "R-HSA-vrel", "R-HSA-vint"]})], ignore_index=True) + umap[U["a2"]] = "R-HSA-GAP2" + profiles["R-HSA-GAP2"] = prof([GAP_P]) + participants += [("R-HSA-vint", "input", "R-HSA-G", 1), ("R-HSA-vint", "output", "R-HSA-T", 1)] + st = {} + _, trans, carriers = m.find_pools(net, rmap, umap, m.r_steps(participants, profiles), profiles, set(), {}, st) + rows = {} + for _, pid, step, a, b, rx, rs, enz in trans: + rows.setdefault((pid, step), []).append((a, b, rx, rs, enz)) + by_shape = {} + for (pid, step), rs_ in rows.items(): + assert len({(a, b, rs, enz) for a, b, _, rs, enz in rs_}) == 1 # one source, target, stId, flag per step + a, b, _, rs, enz = rs_[0] + by_shape.setdefault(pid, []).append((step, a, b, rs, enz, sorted(rx for _, _, rx, _, _ in rs_))) + shapes = sorted(tuple(x[1:] for x in sorted(v)) for v in by_shape.values()) + assert shapes == sorted([ + ((U["g"], U["t"], "R-HSA-vgef", True, [U["vgef"]]),), + ((U["g"], U["t"], "R-HSA-vint", False, [U["vback"]]),), # not merged into the GEF step + ((U["t"], U["g"], "R-HSA-vhyd", False, [U["vhyd"]]),), + ((U["t"], U["tx"], "R-HSA-vbind", True, sorted([U["vbind"], U["v7"]])), + (U["tx"], U["g"], "R-HSA-vrel", False, sorted([U["vrel"], U["v8"]]))), + ]) + assert st["paths"] == 4 and st["copy_rows"] == 7 and st["multi_step_pools"] == 1 + assert sorted(carriers) == sorted([("pool1", U["gap"], U["vrel"]), ("pool1", U["a2"], U["v8"])])