Tracing a TGF-beta fibrosis response in one AnnNet object¶
TGF-beta stimulation of human kidney PDGFRb+ mesenchymal cells produces a timed fibrosis response. This notebook asks which responsive signaling, regulatory, complex and metabolic contexts line up around that matrix-remodeling program.
The graph combines OmniPath signaling, DoRothEA regulation, OmniPath complexes and Human-GEM metabolism in one AnnNet object. The same object is then used for biological joins, receptor-to-effector paths, Cytoscape export and late-response forecasting.
Two aspects carry the experiment: mechanism and time. A node-layer such as
(prot:COL1A1, signaling, 24h) means COL1A1 is represented as a responsive
protein in the signaling layer at 24h.
In [1]:
Copied!
from __future__ import annotations
import json
import time
import warnings
from collections import Counter
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import polars as pl
import pyreadr
import annnet as an
DATA, OUT = Path("data"), Path("outputs")
OMICS = DATA if (DATA / "responsive.parquet").exists() else Path("f/UC2b/data")
TABLES, FIGS = OUT / "tables", OUT / "figures"
TABLES.mkdir(parents=True, exist_ok=True)
FIGS.mkdir(parents=True, exist_ok=True)
SNAPSHOT = DATA / "uc2.annnet"
SEED = 7
ADJP = 0.05
BASELINE = "0h"
TIMES = ("1h", "12h", "24h", "48h", "72h", "96h")
MECHANISMS = ("signaling", "regulatory")
MIN_FRACTION, MIN_HITS = 0.5, 2
RECEPTORS = ("TGFB1", "TGFBR1", "TGFBR2")
EFFECTORS = (
"COL1A1", "COL1A2", "COL3A1", "FN1", "ACTA2", "SERPINE1", "CCN2", "TAGLN",
"TIMP1", "LOX", "THBS1", "TNC", "ITGB1", "MMP2", "VIM", "SMAD7", "JUNB", "ID1",
)
FAMILIES = ("signaling", "regulatory", "coupling_translation", "coupling_time")
an.info()
from __future__ import annotations
import json
import time
import warnings
from collections import Counter
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import polars as pl
import pyreadr
import annnet as an
DATA, OUT = Path("data"), Path("outputs")
OMICS = DATA if (DATA / "responsive.parquet").exists() else Path("f/UC2b/data")
TABLES, FIGS = OUT / "tables", OUT / "figures"
TABLES.mkdir(parents=True, exist_ok=True)
FIGS.mkdir(parents=True, exist_ok=True)
SNAPSHOT = DATA / "uc2.annnet"
SEED = 7
ADJP = 0.05
BASELINE = "0h"
TIMES = ("1h", "12h", "24h", "48h", "72h", "96h")
MECHANISMS = ("signaling", "regulatory")
MIN_FRACTION, MIN_HITS = 0.5, 2
RECEPTORS = ("TGFB1", "TGFBR1", "TGFBR2")
EFFECTORS = (
"COL1A1", "COL1A2", "COL3A1", "FN1", "ACTA2", "SERPINE1", "CCN2", "TAGLN",
"TIMP1", "LOX", "THBS1", "TNC", "ITGB1", "MMP2", "VIM", "SMAD7", "JUNB", "ID1",
)
FAMILIES = ("signaling", "regulatory", "coupling_translation", "coupling_time")
an.info()
Out[1]:
annnetv0.2.0graph: networkxplot: graphvizI/O: 12/12
Show environment details
| Version | v0.2.0 |
|---|---|
| License | BSD-3-Clause |
| Authors | Youssef Zerta ✉, Daniele Bottazzi ✉, Denes Turei ✉ |
| Repository | https://github.com/saezlab/annnet |
| Documentation | https://saezlab.github.io/annnet/reference/ |
| Default adapter backend | networkx |
| Default plot backend | graphviz |
| Graph backends | ✓networkx✓igraph✗graph-tool✓pyg |
| Plot backends | ✓graphviz✓pydot✓matplotlib |
| Tabular data backends | ✓polars✓pandas✓pyarrow |
| I/O modules | ✓annnet✓json/ndjson✓dataframes✓csv✓excel✓graphml/gexf✓sif✓cx2✓parquet✓zarr✓sbml✓scverse |
| Installable bundles | backendsplotiostorageall |
The responsive gate¶
One gate defines the universe: a gene symbol enters a time layer only where the study calls it differential at that timepoint. No fold-change floor, so the early layers stay populated.
In [2]:
Copied!
# Network inputs live in data/; cached response tables live in OMICS.
measurements = pl.read_parquet(OMICS / "measured.parquet")
responsive = pl.read_parquet(OMICS / "responsive.parquet").select(
"symbol", "time", "best_logFC"
)
measured_symbols = set(measurements["symbol"])
responsive_at = {
t: set(group["symbol"]) for t, group in responsive.to_pandas().groupby("time", observed=True)
}
logfc_at = {(r.symbol, r.time): r.best_logFC for r in responsive.to_pandas().itertuples()}
print(f"measured universe: {len(measured_symbols):,} symbols")
print("responsive per timepoint:", {t: len(responsive_at[t]) for t in TIMES})
# Network inputs live in data/; cached response tables live in OMICS.
measurements = pl.read_parquet(OMICS / "measured.parquet")
responsive = pl.read_parquet(OMICS / "responsive.parquet").select(
"symbol", "time", "best_logFC"
)
measured_symbols = set(measurements["symbol"])
responsive_at = {
t: set(group["symbol"]) for t, group in responsive.to_pandas().groupby("time", observed=True)
}
logfc_at = {(r.symbol, r.time): r.best_logFC for r in responsive.to_pandas().itertuples()}
print(f"measured universe: {len(measured_symbols):,} symbols")
print("responsive per timepoint:", {t: len(responsive_at[t]) for t in TIMES})
measured universe: 14,966 symbols
responsive per timepoint: {'1h': 529, '12h': 2762, '24h': 3297, '48h': 4464, '72h': 5078, '96h': 7129}
The aspect grid¶
The two aspects and their elementary layers are declared in the constructor.
The layer space is their product, so every vertex placement below names a
coordinate. G.layers.set_aspects(...) does the same thing after the fact.
In [3]:
Copied!
G = an.AnnNet(
directed=True,
aspects={
"mechanism": ["signaling", "regulatory", "metabolic"],
"time": [BASELINE, *TIMES],
},
)
G.history.enable(True)
G.history.snapshot("init")
print("aspects:", G.layers.list_aspects())
print("elementary layers:", G.layers.list_layers())
G = an.AnnNet(
directed=True,
aspects={
"mechanism": ["signaling", "regulatory", "metabolic"],
"time": [BASELINE, *TIMES],
},
)
G.history.enable(True)
G.history.snapshot("init")
print("aspects:", G.layers.list_aspects())
print("elementary layers:", G.layers.list_layers())
aspects: ('mechanism', 'time')
elementary layers: {'mechanism': ['metabolic', 'regulatory', 'signaling'], 'time': ['0h', '12h', '1h', '24h', '48h', '72h', '96h']}
Prior knowledge¶
The signaling prior is read from the cached CX2 backbone. DoRothEA confidence A/B supplies the regulatory prior.
In [4]:
Copied!
cx2_backbone = json.loads((OUT / "signaling.cx2").read_text())
node_symbol = {
node["id"]: node["v"].get("gene_symbol")
for section in cx2_backbone
for node in section.get("nodes", [])
}
signaling_pairs = []
for section in cx2_backbone:
for edge in section.get("edges", []):
attrs = edge.get("v", {})
if attrs.get("edge_kind") != "signaling":
continue
source, target = node_symbol.get(edge["s"]), node_symbol.get(edge["t"])
if source and target:
signaling_pairs.append((source, target, float(attrs.get("weight", 0.0))))
def edge_sign(stimulation, inhibition):
return 1.0 if stimulation and not inhibition else -1.0 if inhibition and not stimulation else 0.0
dorothea = pl.read_csv(DATA / "dorothea.tsv", separator=" ", infer_schema_length=5000)
regulatory_pairs = [
(r["source_genesymbol"], r["target_genesymbol"], edge_sign(r["is_stimulation"], r["is_inhibition"]))
for r in dorothea.iter_rows(named=True)
]
tf_symbols = set(dorothea["source_genesymbol"])
print(f"prior knowledge: {len(signaling_pairs):,} signaling, {len(regulatory_pairs):,} regulatory")
cx2_backbone = json.loads((OUT / "signaling.cx2").read_text())
node_symbol = {
node["id"]: node["v"].get("gene_symbol")
for section in cx2_backbone
for node in section.get("nodes", [])
}
signaling_pairs = []
for section in cx2_backbone:
for edge in section.get("edges", []):
attrs = edge.get("v", {})
if attrs.get("edge_kind") != "signaling":
continue
source, target = node_symbol.get(edge["s"]), node_symbol.get(edge["t"])
if source and target:
signaling_pairs.append((source, target, float(attrs.get("weight", 0.0))))
def edge_sign(stimulation, inhibition):
return 1.0 if stimulation and not inhibition else -1.0 if inhibition and not stimulation else 0.0
dorothea = pl.read_csv(DATA / "dorothea.tsv", separator=" ", infer_schema_length=5000)
regulatory_pairs = [
(r["source_genesymbol"], r["target_genesymbol"], edge_sign(r["is_stimulation"], r["is_inhibition"]))
for r in dorothea.iter_rows(named=True)
]
tf_symbols = set(dorothea["source_genesymbol"])
print(f"prior knowledge: {len(signaling_pairs):,} signaling, {len(regulatory_pairs):,} regulatory")
prior knowledge: 31,898 signaling, 15,267 regulatory
Responsive layers — signaling and regulatory¶
For each mechanism, walk the six timepoints and keep only the edges whose
endpoints are both responsive then — a pure 0-hop restriction on the gate. The
same protein therefore appears once per timepoint it responds in, as a distinct
supra-node (vertex, mechanism, time).
In [5]:
Copied!
def add_time_layers(G, mechanism, pairs, prefix, kind):
"""Place responsive-gated supra-nodes and edges for one mechanism, per timepoint.
Returns the {(symbol, time)} placed, which the coupling step consumes.
"""
placed = set()
for t in TIMES:
active = responsive_at[t]
kept = [(a, b, w) for a, b, w in pairs if a in active and b in active]
symbols = {v for a, b, _ in kept for v in (a, b)}
G.add_nodes(
[{"node_id": f"{prefix}{v}", "gene_symbol": v, "kind": kind} for v in sorted(symbols)],
layer=(mechanism, t),
)
G.add_edges(
[
{
"source": (f"{prefix}{a}", (mechanism, t)),
"target": (f"{prefix}{b}", (mechanism, t)),
"weight": w,
"edge_kind": mechanism,
}
for a, b, w in kept
],
default_edge_directed=True,
)
placed |= {(v, t) for v in symbols}
return placed
signaling_placed = add_time_layers(G, "signaling", signaling_pairs, "prot:", "protein")
regulatory_placed = add_time_layers(G, "regulatory", regulatory_pairs, "gene:", "gene")
G.history.snapshot("after_response")
# `layer_node_set` / `layer_edge_set` read a node-layer population and its
# intra-layer edge set straight off the graph — no bookkeeping on the side.
topology = pd.DataFrame(
{
f"{mechanism}_{part}": [
len(getattr(G.layers, f"layer_{"node" if part == "Node" else "edge"}_set")((mechanism, t))) for t in TIMES
]
for mechanism in MECHANISMS
for part in ("Node", "Edge")
},
index=list(TIMES),
).rename(columns=lambda c: c.replace("_Node", " |N|").replace("_Edge", " |E|"))
topology.to_csv(TABLES / "per_time_topology.csv")
print(topology.to_string())
def add_time_layers(G, mechanism, pairs, prefix, kind):
"""Place responsive-gated supra-nodes and edges for one mechanism, per timepoint.
Returns the {(symbol, time)} placed, which the coupling step consumes.
"""
placed = set()
for t in TIMES:
active = responsive_at[t]
kept = [(a, b, w) for a, b, w in pairs if a in active and b in active]
symbols = {v for a, b, _ in kept for v in (a, b)}
G.add_nodes(
[{"node_id": f"{prefix}{v}", "gene_symbol": v, "kind": kind} for v in sorted(symbols)],
layer=(mechanism, t),
)
G.add_edges(
[
{
"source": (f"{prefix}{a}", (mechanism, t)),
"target": (f"{prefix}{b}", (mechanism, t)),
"weight": w,
"edge_kind": mechanism,
}
for a, b, w in kept
],
default_edge_directed=True,
)
placed |= {(v, t) for v in symbols}
return placed
signaling_placed = add_time_layers(G, "signaling", signaling_pairs, "prot:", "protein")
regulatory_placed = add_time_layers(G, "regulatory", regulatory_pairs, "gene:", "gene")
G.history.snapshot("after_response")
# `layer_node_set` / `layer_edge_set` read a node-layer population and its
# intra-layer edge set straight off the graph — no bookkeeping on the side.
topology = pd.DataFrame(
{
f"{mechanism}_{part}": [
len(getattr(G.layers, f"layer_{"node" if part == "Node" else "edge"}_set")((mechanism, t))) for t in TIMES
]
for mechanism in MECHANISMS
for part in ("Node", "Edge")
},
index=list(TIMES),
).rename(columns=lambda c: c.replace("_Node", " |N|").replace("_Edge", " |E|"))
topology.to_csv(TABLES / "per_time_topology.csv")
print(topology.to_string())
signaling |N| signaling |E| regulatory |N| regulatory |E| 1h 105 132 116 158 12h 813 1714 513 838 24h 1125 2643 750 1346 48h 1645 4858 1063 2044 72h 1950 5961 1221 2184 96h 2974 10884 1840 3592
Interlayer coupling¶
Two families, and the distinction between them is the whole point of having two aspects.
Time coupling is ordinal: it links the same entity at consecutive
timepoints, (v, mech, tᵢ) → (v, mech, tᵢ₊₁), wherever it is responsive in
both. Time runs one way, so these edges are directed.
Translation coupling is categorical: it links a responsive gene to its
protein inside one timepoint, (v, regulatory, t) — (v, signaling, t). This is
an identity correspondence between two readings of the same entity, not a
reaction, so it is undirected — which is also what makes the causal chain
receptor → TF protein → TF-as-regulator → target genes traversable in the
propagation section.
In [6]:
Copied!
def at(placed, t):
"""Symbols placed in layer `t`."""
return {v for v, tt in placed if tt == t}
time_edges = [
{
"edge_id": f"time:{mechanism}:{v}:{a}->{b}",
"edge_kind": "coupling_time",
"weight": 1.0,
"source": (f"{prefix}{v}", (mechanism, a)),
"target": (f"{prefix}{v}", (mechanism, b)),
}
for mechanism, prefix, placed in [
("signaling", "prot:", signaling_placed),
("regulatory", "gene:", regulatory_placed),
]
for a, b in zip(TIMES, TIMES[1:])
for v in sorted(at(placed, a) & at(placed, b))
]
translation_edges = [
{
"edge_id": f"trans:{v}:{t}",
"edge_kind": "coupling_translation",
"weight": 1.0,
"source": (f"gene:{v}", ("regulatory", t)),
"target": (f"prot:{v}", ("signaling", t)),
}
for t in TIMES
for v in sorted(at(regulatory_placed, t) & at(signaling_placed, t))
]
G.add_edges(time_edges, default_edge_directed=True)
G.add_edges(translation_edges, default_edge_directed=False)
G.history.snapshot("after_coupling")
print(f"time-coupling edges (directed) : {len(time_edges):,}")
print(f"translation-coupling edges (undirected): {len(translation_edges):,}")
def at(placed, t):
"""Symbols placed in layer `t`."""
return {v for v, tt in placed if tt == t}
time_edges = [
{
"edge_id": f"time:{mechanism}:{v}:{a}->{b}",
"edge_kind": "coupling_time",
"weight": 1.0,
"source": (f"{prefix}{v}", (mechanism, a)),
"target": (f"{prefix}{v}", (mechanism, b)),
}
for mechanism, prefix, placed in [
("signaling", "prot:", signaling_placed),
("regulatory", "gene:", regulatory_placed),
]
for a, b in zip(TIMES, TIMES[1:])
for v in sorted(at(placed, a) & at(placed, b))
]
translation_edges = [
{
"edge_id": f"trans:{v}:{t}",
"edge_kind": "coupling_translation",
"weight": 1.0,
"source": (f"gene:{v}", ("regulatory", t)),
"target": (f"prot:{v}", ("signaling", t)),
}
for t in TIMES
for v in sorted(at(regulatory_placed, t) & at(signaling_placed, t))
]
G.add_edges(time_edges, default_edge_directed=True)
G.add_edges(translation_edges, default_edge_directed=False)
G.history.snapshot("after_coupling")
print(f"time-coupling edges (directed) : {len(time_edges):,}")
print(f"translation-coupling edges (undirected): {len(translation_edges):,}")
time-coupling edges (directed) : 7,644 translation-coupling edges (undirected): 3,314
Baseline scaffold — two kinds of hyperedge¶
Complexes and metabolism are time-invariant, so they live in the 0h layer.
Each OmniPath complex becomes one undirected hyperedge over its subunits —
a single incidence-matrix column, not a clique of pairwise edges. One
from_sbml call reads Human-GEM: reactions become signed directed
hyperedges carrying stoichiometry, and compartments become slices.
from_sbml adds metabolite and boundary vertices without a kind, so they are
typed here for the PyG export further down.
In [7]:
Copied!
complexes = pl.read_csv(DATA / "omnipath_complexes.tsv", separator="\t", infer_schema_length=5000)
complex_specs, members = [], set()
for row in complexes.iter_rows(named=True):
subunits = (row["components_genesymbols"] or "").split("_")
if subunits == [""] or not set(subunits) <= measured_symbols:
continue
members.update(f"prot:{s}" for s in subunits)
complex_specs.append(
{
"edge_id": f"cpx:{row['name']}",
"edge_kind": "complex",
"weight": 1.0,
"members": [(f"prot:{s}", ("signaling", BASELINE)) for s in subunits],
}
)
G.add_nodes(
[
{"node_id": v, "kind": "protein", "gene_symbol": v.removeprefix("prot:")}
for v in sorted(members)
],
layer=("signaling", BASELINE),
)
G.add_edges(complex_specs, layer=("signaling", BASELINE))
G.history.snapshot("after_complex")
an.from_sbml(
str(DATA / "Human-GEM.xml"),
graph=G,
slice="metabolic",
layer=("metabolic", BASELINE),
preserve_stoichiometry=True,
)
untyped = (
G.views.nodes()
.filter(pl.col("kind").is_null())
.get_column("node_id")
.to_list()
)
G.attrs.set_node_attrs_bulk(
{v: {"kind": "boundary" if v.startswith("__") else "metabolite"} for v in untyped}
)
G.history.snapshot("after_metabolic")
# Complex names are not unique in OmniPath, so same-named entries share an edge_id.
n_complex = sum(str(e).startswith("cpx:") for e in G.hyperedge_definitions)
print(f"complex hyperedges : {n_complex:,} (from {len(complex_specs):,} specs)")
print(f"reaction hyperedges : {len(G.hyperedge_definitions) - n_complex:,}")
print(f"|V| = {G.global_count('nodes'):,} |E| = {G.global_count('edges'):,}")
complexes = pl.read_csv(DATA / "omnipath_complexes.tsv", separator="\t", infer_schema_length=5000)
complex_specs, members = [], set()
for row in complexes.iter_rows(named=True):
subunits = (row["components_genesymbols"] or "").split("_")
if subunits == [""] or not set(subunits) <= measured_symbols:
continue
members.update(f"prot:{s}" for s in subunits)
complex_specs.append(
{
"edge_id": f"cpx:{row['name']}",
"edge_kind": "complex",
"weight": 1.0,
"members": [(f"prot:{s}", ("signaling", BASELINE)) for s in subunits],
}
)
G.add_nodes(
[
{"node_id": v, "kind": "protein", "gene_symbol": v.removeprefix("prot:")}
for v in sorted(members)
],
layer=("signaling", BASELINE),
)
G.add_edges(complex_specs, layer=("signaling", BASELINE))
G.history.snapshot("after_complex")
an.from_sbml(
str(DATA / "Human-GEM.xml"),
graph=G,
slice="metabolic",
layer=("metabolic", BASELINE),
preserve_stoichiometry=True,
)
untyped = (
G.views.nodes()
.filter(pl.col("kind").is_null())
.get_column("node_id")
.to_list()
)
G.attrs.set_node_attrs_bulk(
{v: {"kind": "boundary" if v.startswith("__") else "metabolite"} for v in untyped}
)
G.history.snapshot("after_metabolic")
# Complex names are not unique in OmniPath, so same-named entries share an edge_id.
n_complex = sum(str(e).startswith("cpx:") for e in G.hyperedge_definitions)
print(f"complex hyperedges : {n_complex:,} (from {len(complex_specs):,} specs)")
print(f"reaction hyperedges : {len(G.hyperedge_definitions) - n_complex:,}")
print(f"|V| = {G.global_count('nodes'):,} |E| = {G.global_count('edges'):,}")
complex hyperedges : 3,446 (from 7,991 specs) reaction hyperedges : 12,971 |V| = 19,580 |E| = 63,729
Organelle slices¶
A slice is a membership overlay on the shared graph: it duplicates no topology, where the usual alternative keeps one graph object per compartment. Compartments are time-invariant, so a protein joins its organelle regardless of when it responds. The cross-compartment query below reads the membership back out of the graph.
In [8]:
Copied!
ORGANELLES = {
"Mitochondria": "mitochondria",
"Nucleus": "nucleus",
"Nucleoplasm": "nucleus",
"Nucleoli": "nucleus",
"Nuclear membrane": "nucleus",
"Endoplasmic reticulum": "er",
"Golgi apparatus": "golgi",
"Lysosome": "lysosome",
"Cytosol": "cytosol",
"Cytoplasm": "cytosol",
"Plasma membrane": "plasma_membrane",
"Cell Junctions": "plasma_membrane",
"Peroxisome": "peroxisome",
"Vesicles": "vesicles",
}
hpa = pl.read_csv(DATA / "proteinatlas.tsv", separator="\t", infer_schema_length=10000)
location = dict(zip(hpa["Gene"].to_list(), hpa["Subcellular main location"].to_list()))
def organelle_of(symbol):
"""Organelle slice name for a gene symbol; cytosol by default."""
text = location.get(symbol) or ""
return next((name for key, name in ORGANELLES.items() if key.lower() in text.lower()), "cytosol")
proteins = G.views.nodes().filter(pl.col("node_id").str.starts_with("prot:"))
organelle = {
r["node_id"]: organelle_of(r["gene_symbol"]) for r in proteins.iter_rows(named=True)
}
for name in sorted(set(organelle.values())):
vids = {v for v, o in organelle.items() if o == name}
G.slices.add(name, role="organelle")
for v in sorted(vids):
G.slices.add_node_to_slice(name, v)
G.slices.add_edges(
name,
[e for e, (s, t, _) in G.edge_definitions.items() if s[0] in vids and t[0] in vids],
)
G.history.snapshot("after_slices")
print(f"{len(set(organelle.values()))} organelle slices over {len(organelle):,} proteins")
print(Counter(organelle.values()).most_common())
ORGANELLES = {
"Mitochondria": "mitochondria",
"Nucleus": "nucleus",
"Nucleoplasm": "nucleus",
"Nucleoli": "nucleus",
"Nuclear membrane": "nucleus",
"Endoplasmic reticulum": "er",
"Golgi apparatus": "golgi",
"Lysosome": "lysosome",
"Cytosol": "cytosol",
"Cytoplasm": "cytosol",
"Plasma membrane": "plasma_membrane",
"Cell Junctions": "plasma_membrane",
"Peroxisome": "peroxisome",
"Vesicles": "vesicles",
}
hpa = pl.read_csv(DATA / "proteinatlas.tsv", separator="\t", infer_schema_length=10000)
location = dict(zip(hpa["Gene"].to_list(), hpa["Subcellular main location"].to_list()))
def organelle_of(symbol):
"""Organelle slice name for a gene symbol; cytosol by default."""
text = location.get(symbol) or ""
return next((name for key, name in ORGANELLES.items() if key.lower() in text.lower()), "cytosol")
proteins = G.views.nodes().filter(pl.col("node_id").str.starts_with("prot:"))
organelle = {
r["node_id"]: organelle_of(r["gene_symbol"]) for r in proteins.iter_rows(named=True)
}
for name in sorted(set(organelle.values())):
vids = {v for v, o in organelle.items() if o == name}
G.slices.add(name, role="organelle")
for v in sorted(vids):
G.slices.add_node_to_slice(name, v)
G.slices.add_edges(
name,
[e for e, (s, t, _) in G.edge_definitions.items() if s[0] in vids and t[0] in vids],
)
G.history.snapshot("after_slices")
print(f"{len(set(organelle.values()))} organelle slices over {len(organelle):,} proteins")
print(Counter(organelle.values()).most_common())
9 organelle slices over 8,999 proteins
[('cytosol', 3465), ('nucleus', 3249), ('mitochondria', 609), ('plasma_membrane', 513), ('vesicles', 489), ('golgi', 369), ('er', 277), ('lysosome', 17), ('peroxisome', 11)]
Part 2 — Biological joins¶
The first readouts look for structure around the fibrotic program before any shortest-path or learning model is used.
In [9]:
Copied!
edge_kind = G.attrs.get_attr_from_edges("edge_kind")
E = pd.DataFrame(
[(e, s[0], *s[1], t[0], *t[1], edge_kind[e]) for e, (s, t, _) in G.edge_definitions.items()],
columns=["edge_id", "src", "src_mech", "src_time", "tgt", "tgt_mech", "tgt_time", "edge_kind"],
)
intra = E[(E.src_mech == E.tgt_mech) & (E.src_time == E.tgt_time)]
subunits = {
e: {m[0].removeprefix("prot:") for m in spec["members"]}
for e, spec in G.hyperedge_definitions.items()
if str(e).startswith("cpx:")
}
print(E.edge_kind.value_counts().to_string())
print(f"\nintra-layer {len(intra):,} (both endpoints in one (mechanism, time) layer)")
print(f"interlayer {len(E) - len(intra):,} (the two coupling families)")
print(f"complexes {len(subunits):,} hyperedges")
edge_kind = G.attrs.get_attr_from_edges("edge_kind")
E = pd.DataFrame(
[(e, s[0], *s[1], t[0], *t[1], edge_kind[e]) for e, (s, t, _) in G.edge_definitions.items()],
columns=["edge_id", "src", "src_mech", "src_time", "tgt", "tgt_mech", "tgt_time", "edge_kind"],
)
intra = E[(E.src_mech == E.tgt_mech) & (E.src_time == E.tgt_time)]
subunits = {
e: {m[0].removeprefix("prot:") for m in spec["members"]}
for e, spec in G.hyperedge_definitions.items()
if str(e).startswith("cpx:")
}
print(E.edge_kind.value_counts().to_string())
print(f"\nintra-layer {len(intra):,} (both endpoints in one (mechanism, time) layer)")
print(f"interlayer {len(E) - len(intra):,} (the two coupling families)")
print(f"complexes {len(subunits):,} hyperedges")
edge_kind signaling 26192 regulatory 10162 coupling_time 7644 coupling_translation 3314 intra-layer 36,354 (both endpoints in one (mechanism, time) layer) interlayer 10,958 (the two coupling families) complexes 3,446 hyperedges
Q1 — Which TFs organize responsive protein complexes?¶
For each timepoint, a TF is paired with an OmniPath complex when its responsive DoRothEA targets cover enough measured subunits of that complex. This highlights candidate regulatory programs whose targets appear as coordinated protein machinery, not just individual genes.
In [10]:
Copied!
regulatory = intra[intra.edge_kind == "regulatory"].assign(
src_symbol=lambda d: d.src.str.removeprefix("gene:"),
tgt_symbol=lambda d: d.tgt.str.removeprefix("gene:"),
)
big_complexes = {e: g for e, g in subunits.items() if len(g) >= 3}
def coregulated(targets_of, complex_sets=big_complexes):
"""(TF, complex) pairs whose targets cover enough of the complex."""
return {
(tf, complex_id)
for complex_id, genes in complex_sets.items()
for tf, targets in targets_of.items()
if len(genes & targets) >= MIN_HITS and len(genes & targets) / len(genes) >= MIN_FRACTION
}
targets_at = {
t: regulatory[regulatory.src_time == t].groupby("src_symbol")["tgt_symbol"].apply(set)
for t in TIMES
}
pairs_at = {t: coregulated(targets_at[t]) for t in TIMES}
print("responsive TF-complex programs:", {t: len(p) for t, p in pairs_at.items()})
late_programs = pd.DataFrame(sorted(pairs_at[TIMES[-1]]), columns=["tf", "complex_id"])
late_programs["fibrosis_subunits"] = late_programs.complex_id.map(
lambda complex_id: ",".join(sorted(big_complexes[complex_id] & set(EFFECTORS)))
)
late_programs = late_programs.sort_values(
["fibrosis_subunits", "tf", "complex_id"], ascending=[False, True, True]
)
late_programs.to_csv(TABLES / f"tf_complex_programs_{TIMES[-1]}.csv", index=False)
print(f"\nlate programs with fibrosis-effector subunits first:")
print(late_programs.head(10).to_string(index=False))
regulatory = intra[intra.edge_kind == "regulatory"].assign(
src_symbol=lambda d: d.src.str.removeprefix("gene:"),
tgt_symbol=lambda d: d.tgt.str.removeprefix("gene:"),
)
big_complexes = {e: g for e, g in subunits.items() if len(g) >= 3}
def coregulated(targets_of, complex_sets=big_complexes):
"""(TF, complex) pairs whose targets cover enough of the complex."""
return {
(tf, complex_id)
for complex_id, genes in complex_sets.items()
for tf, targets in targets_of.items()
if len(genes & targets) >= MIN_HITS and len(genes & targets) / len(genes) >= MIN_FRACTION
}
targets_at = {
t: regulatory[regulatory.src_time == t].groupby("src_symbol")["tgt_symbol"].apply(set)
for t in TIMES
}
pairs_at = {t: coregulated(targets_at[t]) for t in TIMES}
print("responsive TF-complex programs:", {t: len(p) for t, p in pairs_at.items()})
late_programs = pd.DataFrame(sorted(pairs_at[TIMES[-1]]), columns=["tf", "complex_id"])
late_programs["fibrosis_subunits"] = late_programs.complex_id.map(
lambda complex_id: ",".join(sorted(big_complexes[complex_id] & set(EFFECTORS)))
)
late_programs = late_programs.sort_values(
["fibrosis_subunits", "tf", "complex_id"], ascending=[False, True, True]
)
late_programs.to_csv(TABLES / f"tf_complex_programs_{TIMES[-1]}.csv", index=False)
print(f"\nlate programs with fibrosis-effector subunits first:")
print(late_programs.head(10).to_string(index=False))
responsive TF-complex programs: {'1h': 3, '12h': 53, '24h': 63, '48h': 130, '72h': 132, '96h': 278}
late programs with fibrosis-effector subunits first:
tf complex_id fibrosis_subunits
JUND cpx:ITGAV-ITGB3-THBS1 complex THBS1
SP1 cpx:Ubiquitin E3 ligase (SMAD7-SMURF1) - TGF(beta) SMAD7
SP1 cpx:Ubiquitin E3 ligase (SMURF2, SMAD7) - TGF(beta) SMAD7
E2F1 cpx:ITGA3-ITGB1-THBS1 complex ITGB1,THBS1
E2F1 cpx:ITGA4-ITGB1-THBS1 complex ITGB1,THBS1
MYC cpx:ITGA3-ITGB1-THBS1 complex ITGB1,THBS1
MYC cpx:ITGA4-ITGB1-THBS1 complex ITGB1,THBS1
MYC cpx:ITGA3-ITGB1-BSG complex ITGB1
MYC cpx:ITGA6-ITGB1-CD151 complex ITGB1
MYC cpx:ITGA6-ITGB1-CYR61 complex ITGB1
Q2 — Where does responsive signaling cross compartments?¶
Organelle slices turn each responsive signaling edge into a compartment-to- compartment observation. The table asks whether the signaling layer is mostly local or whether the response increasingly crosses organelle boundaries, which is relevant for receptor, cytosolic kinase and nuclear TF handoffs.
In [11]:
Copied!
organelle_of_vertex = {
v: name
for name in G.slices.list(include_default=False)
if G.attrs.get_slice_attr(name, "role") == "organelle"
for v in G.slices.nodes(name)
}
print(f"{len(set(organelle_of_vertex.values()))} organelle slices covering "
f"{len(organelle_of_vertex):,} proteins\n")
signaling = intra[intra.edge_kind == "signaling"].assign(
src_organelle=lambda d: d.src.map(organelle_of_vertex),
tgt_organelle=lambda d: d.tgt.map(organelle_of_vertex),
)
cross = signaling[signaling.src_organelle != signaling.tgt_organelle]
per_time = pd.DataFrame(
{
"signaling": signaling.groupby("src_time").size().reindex(TIMES, fill_value=0),
"cross_compartment": cross.groupby("src_time").size().reindex(TIMES, fill_value=0),
}
)
per_time["fraction"] = (per_time.cross_compartment / per_time.signaling).round(3)
per_time.to_csv(TABLES / "cross_compartment.csv")
print(per_time.to_string())
organelle_of_vertex = {
v: name
for name in G.slices.list(include_default=False)
if G.attrs.get_slice_attr(name, "role") == "organelle"
for v in G.slices.nodes(name)
}
print(f"{len(set(organelle_of_vertex.values()))} organelle slices covering "
f"{len(organelle_of_vertex):,} proteins\n")
signaling = intra[intra.edge_kind == "signaling"].assign(
src_organelle=lambda d: d.src.map(organelle_of_vertex),
tgt_organelle=lambda d: d.tgt.map(organelle_of_vertex),
)
cross = signaling[signaling.src_organelle != signaling.tgt_organelle]
per_time = pd.DataFrame(
{
"signaling": signaling.groupby("src_time").size().reindex(TIMES, fill_value=0),
"cross_compartment": cross.groupby("src_time").size().reindex(TIMES, fill_value=0),
}
)
per_time["fraction"] = (per_time.cross_compartment / per_time.signaling).round(3)
per_time.to_csv(TABLES / "cross_compartment.csv")
print(per_time.to_string())
9 organelle slices covering 8,999 proteins
signaling cross_compartment fraction
src_time
1h 132 81 0.614
12h 1714 1166 0.680
24h 2643 1770 0.670
48h 4858 3139 0.646
72h 5961 3793 0.636
96h 10884 6537 0.601
Q3 — Which complexes sit near the responsive signaling backbone?¶
Each baseline complex is scored by the signaling degree of its measured subunits across the time-resolved signaling layers. Complexes containing fibrosis effectors are kept visible in the same table, so the query can point to machinery that is both topologically exposed and biologically on-theme.
In [12]:
Copied!
degree = (
pd.concat([signaling.src.value_counts(), signaling.tgt.value_counts()], axis=1)
.fillna(0)
.sum(axis=1)
)
hubs = (
pd.DataFrame(
[
{
"complex_id": complex_id,
"n_subunits": len(genes),
"signaling_degree": float(sum(degree.get(f"prot:{g}", 0) for g in genes)),
"fibrosis_subunits": ",".join(sorted(genes & set(EFFECTORS))),
}
for complex_id, genes in subunits.items()
]
)
.sort_values(["fibrosis_subunits", "signaling_degree"], ascending=[False, False])
.head(10)
)
hubs.to_csv(TABLES / "hub_complexes.csv", index=False)
print(hubs.head(5).to_string(index=False))
degree = (
pd.concat([signaling.src.value_counts(), signaling.tgt.value_counts()], axis=1)
.fillna(0)
.sum(axis=1)
)
hubs = (
pd.DataFrame(
[
{
"complex_id": complex_id,
"n_subunits": len(genes),
"signaling_degree": float(sum(degree.get(f"prot:{g}", 0) for g in genes)),
"fibrosis_subunits": ",".join(sorted(genes & set(EFFECTORS))),
}
for complex_id, genes in subunits.items()
]
)
.sort_values(["fibrosis_subunits", "signaling_degree"], ascending=[False, False])
.head(10)
)
hubs.to_csv(TABLES / "hub_complexes.csv", index=False)
print(hubs.head(5).to_string(index=False))
complex_id n_subunits signaling_degree fibrosis_subunits
cpx:pERK-vimentin-KPNA2 complex 3 498.0 VIM
cpx:Tenascin-C complex 1 42.0 TNC
cpx:ITGAV-ITGB3-THBS1 complex 3 295.0 THBS1
cpx:Thrombospondin 1 complex 1 49.0 THBS1
cpx:Ubiquitin E3 ligase (SMAD7-SMURF1) - TGF(beta) 4 367.0 SMAD7
Part 4 — Receptor-to-effector paths¶
Starting from TGF-beta receptors at 1h, which fibrosis effectors become reachable, through which intermediates, and at what timepoint?
In [13]:
Copied!
t0 = time.perf_counter()
with warnings.catch_warnings():
warnings.simplefilter("ignore", RuntimeWarning)
flow = G.ig.backend()
keys = [str(n) for n in flow.vs["name"]]
flow = flow.permute_vertices(np.argsort(keys, kind="stable").tolist())
names = flow.vs["name"]
supra = {n: i for i, n in enumerate(names) if isinstance(n, tuple)}
seeds = [supra[(f"prot:{r}", ("signaling", TIMES[0]))] for r in RECEPTORS
if (f"prot:{r}", ("signaling", TIMES[0])) in supra]
print(f"AnnNet -> igraph: {time.perf_counter() - t0:.1f}s "
f"({flow.vcount():,} vertices, {flow.ecount():,} arcs)")
print(f"seeds: {[names[i] for i in seeds]}")
print(Counter(flow.es["edge_kind"]).most_common())
t0 = time.perf_counter()
with warnings.catch_warnings():
warnings.simplefilter("ignore", RuntimeWarning)
flow = G.ig.backend()
keys = [str(n) for n in flow.vs["name"]]
flow = flow.permute_vertices(np.argsort(keys, kind="stable").tolist())
names = flow.vs["name"]
supra = {n: i for i, n in enumerate(names) if isinstance(n, tuple)}
seeds = [supra[(f"prot:{r}", ("signaling", TIMES[0]))] for r in RECEPTORS
if (f"prot:{r}", ("signaling", TIMES[0])) in supra]
print(f"AnnNet -> igraph: {time.perf_counter() - t0:.1f}s "
f"({flow.vcount():,} vertices, {flow.ecount():,} arcs)")
print(f"seeds: {[names[i] for i in seeds]}")
print(Counter(flow.es["edge_kind"]).most_common())
AnnNet -> igraph: 2.8s (46,109 vertices, 50,626 arcs)
seeds: [('prot:TGFBR1', ('signaling', '1h')), ('prot:TGFBR2', ('signaling', '1h'))]
[('signaling', 26192), ('regulatory', 10162), ('coupling_time', 7644), ('coupling_translation', 6628)]
Coupling ablation¶
The ablation asks which edge families are needed for receptor-origin paths to reach regulatory effector genes. The middle rows are diagnostic: signaling can propagate among proteins, time coupling can carry that signal forward, and translation coupling is what lets a reached TF protein become a regulator of responsive target genes.
In [14]:
Copied!
def reachable(families):
kept = [e.index for e in flow.es if e["edge_kind"] in families]
sub = flow.subgraph_edges(kept, delete_vertices=False)
return np.asarray(sub.distances(source=seeds, mode="out")).min(axis=0)
def earliest(distance, symbol, mechanism, prefix):
hits = [
(t, distance[supra[(f"{prefix}{symbol}", (mechanism, t))]])
for t in TIMES
if (f"{prefix}{symbol}", (mechanism, t)) in supra
]
hits = [(t, d) for t, d in hits if np.isfinite(d)]
return min(hits, key=lambda x: (x[1], TIMES.index(x[0]))) if hits else None
configs = [
("signaling", ("signaling",)),
("+ regulatory", ("signaling", "regulatory")),
("+ time coupling", ("signaling", "regulatory", "coupling_time")),
("+ translation", FAMILIES),
]
rows, distances = [], {}
for tag, families in configs:
d = reachable(families)
distances[tag] = d
reached = [s for s in EFFECTORS if earliest(d, s, "regulatory", "gene:")]
rows.append({"configuration": tag, "supra_nodes": int(np.isfinite(d).sum()),
"effectors": len(reached)})
ablation = pd.DataFrame(rows)
ablation.to_csv(TABLES / "signal_flow_ablation.csv", index=False)
print(ablation.to_string(index=False))
def reachable(families):
kept = [e.index for e in flow.es if e["edge_kind"] in families]
sub = flow.subgraph_edges(kept, delete_vertices=False)
return np.asarray(sub.distances(source=seeds, mode="out")).min(axis=0)
def earliest(distance, symbol, mechanism, prefix):
hits = [
(t, distance[supra[(f"{prefix}{symbol}", (mechanism, t))]])
for t in TIMES
if (f"{prefix}{symbol}", (mechanism, t)) in supra
]
hits = [(t, d) for t, d in hits if np.isfinite(d)]
return min(hits, key=lambda x: (x[1], TIMES.index(x[0]))) if hits else None
configs = [
("signaling", ("signaling",)),
("+ regulatory", ("signaling", "regulatory")),
("+ time coupling", ("signaling", "regulatory", "coupling_time")),
("+ translation", FAMILIES),
]
rows, distances = [], {}
for tag, families in configs:
d = reachable(families)
distances[tag] = d
reached = [s for s in EFFECTORS if earliest(d, s, "regulatory", "gene:")]
rows.append({"configuration": tag, "supra_nodes": int(np.isfinite(d).sum()),
"effectors": len(reached)})
ablation = pd.DataFrame(rows)
ablation.to_csv(TABLES / "signal_flow_ablation.csv", index=False)
print(ablation.to_string(index=False))
configuration supra_nodes effectors
signaling 3 0
+ regulatory 3 0
+ time coupling 6164 0
+ translation 12264 18
Earliest arrival, and the path that delivers it¶
With all four families in play, distances gives each effector the timepoint at
which it first becomes reachable and the number of hops it takes;
get_shortest_paths gives the actual chain of supra-nodes.
In [15]:
Copied!
full = distances["+ translation"]
arrival = pd.DataFrame(
[
{"effector": s, "first_reachable": layer, "hops": int(hops)}
for s in EFFECTORS
for layer, hops in [earliest(full, s, "regulatory", "gene:") or (None, np.nan)]
if layer is not None
]
).sort_values(["first_reachable", "hops", "effector"])
arrival.to_csv(TABLES / "effector_arrival.csv", index=False)
print(arrival.to_string(index=False))
paths = {}
for r in arrival.itertuples():
target = supra[(f"gene:{r.effector}", ("regulatory", r.first_reachable))]
source = min(seeds, key=lambda s: flow.distances(source=[s], target=[target], mode="out")[0][0])
paths[r.effector] = flow.get_shortest_paths(source, to=target, mode="out")[0]
for s in ("COL1A1", "SERPINE1", "ACTA2"):
chain = " -> ".join(f"{names[v][0]}@{names[v][1][1]}" for v in paths[s])
print(f"\n{s}: {chain}")
full = distances["+ translation"]
arrival = pd.DataFrame(
[
{"effector": s, "first_reachable": layer, "hops": int(hops)}
for s in EFFECTORS
for layer, hops in [earliest(full, s, "regulatory", "gene:") or (None, np.nan)]
if layer is not None
]
).sort_values(["first_reachable", "hops", "effector"])
arrival.to_csv(TABLES / "effector_arrival.csv", index=False)
print(arrival.to_string(index=False))
paths = {}
for r in arrival.itertuples():
target = supra[(f"gene:{r.effector}", ("regulatory", r.first_reachable))]
source = min(seeds, key=lambda s: flow.distances(source=[s], target=[target], mode="out")[0][0])
paths[r.effector] = flow.get_shortest_paths(source, to=target, mode="out")[0]
for s in ("COL1A1", "SERPINE1", "ACTA2"):
chain = " -> ".join(f"{names[v][0]}@{names[v][1][1]}" for v in paths[s])
print(f"\n{s}: {chain}")
effector first_reachable hops
COL1A1 12h 4
SERPINE1 12h 4
SMAD7 12h 4
TNC 12h 4
FN1 12h 5
ITGB1 12h 5
MMP2 12h 5
TAGLN 12h 5
VIM 12h 5
ID1 12h 6
CCN2 12h 7
LOX 12h 7
JUNB 24h 6
THBS1 24h 6
ACTA2 24h 7
COL1A2 24h 8
COL3A1 72h 7
TIMP1 72h 7
COL1A1: prot:TGFBR2@1h -> prot:TGFBR2@12h -> prot:SMAD3@12h -> gene:SMAD3@12h -> gene:COL1A1@12h
SERPINE1: prot:TGFBR2@1h -> prot:TGFBR2@12h -> prot:SMAD3@12h -> gene:SMAD3@12h -> gene:SERPINE1@12h
ACTA2: prot:TGFBR2@1h -> prot:TGFBR2@12h -> prot:PARD6A@12h -> prot:PARD6A@24h -> prot:PRKCA@24h -> prot:TP53@24h -> gene:TP53@24h -> gene:ACTA2@24h
Effector arrival¶
In [16]:
Copied!
COLOURS = {
"signaling": "#4c72b0", "regulatory": "#55a868",
"coupling_time": "#c44e52", "coupling_translation": "#dd8452",
}
BAND = {"regulatory": 0, "signaling": 1}
nodes = sorted({v for path in paths.values() for v in path})
arcs = sorted({edge for path in paths.values() for edge in zip(path, path[1:])})
endpoints = {path[-1] for path in paths.values()}
arc_kind = {(a, b): flow.es[flow.get_eid(a, b)]["edge_kind"] for a, b in arcs}
slots = {}
for v in nodes:
mechanism, layer = names[v][1]
slot = (mechanism, layer)
rank = slots.setdefault(slot, {})
rank[v] = len(rank)
pos = {}
for v in nodes:
mechanism, layer = names[v][1]
k = slots[(mechanism, layer)][v]
n = len(slots[(mechanism, layer)])
pos[v] = (TIMES.index(layer), BAND[mechanism] + (k - (n - 1) / 2) * 0.055)
fig, ax = plt.subplots(figsize=(12, 5.5))
for mechanism, y in BAND.items():
ax.axhspan(y - 0.35, y + 0.35, color=COLOURS[mechanism], alpha=0.06)
for x in np.arange(len(TIMES) - 1) + 0.5:
ax.axvline(x, color="#dddddd", lw=0.7, ls=":")
for a, b in arcs:
(x0, y0), (x1, y1) = pos[a], pos[b]
ax.annotate("", (x1, y1), (x0, y0), arrowprops={
"arrowstyle": "-|>", "lw": 0.9, "color": COLOURS.get(arc_kind[(a, b)], "#999999"),
"alpha": 0.8, "connectionstyle": "arc3,rad=0.12",
})
for v in nodes:
symbol = names[v][0].split(":", 1)[1]
style = ("#22223b", "white", "bold") if v in seeds else ("#dd8452", "black", "bold") if v in endpoints else ("white", "black", "normal")
ax.text(*pos[v], symbol, ha="center", va="center", fontsize=6.8,
color=style[1], fontweight=style[2],
bbox=dict(boxstyle="round,pad=0.2", fc=style[0], ec="#888888", lw=0.7))
ax.set(xticks=range(len(TIMES)), xticklabels=TIMES, yticks=[0, 1],
yticklabels=["regulatory\n(genes)", "signaling\n(proteins)"],
xlim=(-0.6, len(TIMES) - 0.4), ylim=(-0.45, 1.45), xlabel="time after TGF-beta")
ax.set_title("TGF-beta receptor-to-effector shortest paths", fontsize=11)
ax.legend([plt.Line2D([], [], color=c, lw=2) for c in COLOURS.values()],
[k.replace("coupling_", "") for k in COLOURS], ncol=4, frameon=False,
loc="upper center", bbox_to_anchor=(0.5, -0.11), fontsize=9)
for spine in ("top", "right", "left"):
ax.spines[spine].set_visible(False)
fig.tight_layout()
fig.savefig(FIGS / "signal_flow_paths.png", dpi=150)
plt.show()
COLOURS = {
"signaling": "#4c72b0", "regulatory": "#55a868",
"coupling_time": "#c44e52", "coupling_translation": "#dd8452",
}
BAND = {"regulatory": 0, "signaling": 1}
nodes = sorted({v for path in paths.values() for v in path})
arcs = sorted({edge for path in paths.values() for edge in zip(path, path[1:])})
endpoints = {path[-1] for path in paths.values()}
arc_kind = {(a, b): flow.es[flow.get_eid(a, b)]["edge_kind"] for a, b in arcs}
slots = {}
for v in nodes:
mechanism, layer = names[v][1]
slot = (mechanism, layer)
rank = slots.setdefault(slot, {})
rank[v] = len(rank)
pos = {}
for v in nodes:
mechanism, layer = names[v][1]
k = slots[(mechanism, layer)][v]
n = len(slots[(mechanism, layer)])
pos[v] = (TIMES.index(layer), BAND[mechanism] + (k - (n - 1) / 2) * 0.055)
fig, ax = plt.subplots(figsize=(12, 5.5))
for mechanism, y in BAND.items():
ax.axhspan(y - 0.35, y + 0.35, color=COLOURS[mechanism], alpha=0.06)
for x in np.arange(len(TIMES) - 1) + 0.5:
ax.axvline(x, color="#dddddd", lw=0.7, ls=":")
for a, b in arcs:
(x0, y0), (x1, y1) = pos[a], pos[b]
ax.annotate("", (x1, y1), (x0, y0), arrowprops={
"arrowstyle": "-|>", "lw": 0.9, "color": COLOURS.get(arc_kind[(a, b)], "#999999"),
"alpha": 0.8, "connectionstyle": "arc3,rad=0.12",
})
for v in nodes:
symbol = names[v][0].split(":", 1)[1]
style = ("#22223b", "white", "bold") if v in seeds else ("#dd8452", "black", "bold") if v in endpoints else ("white", "black", "normal")
ax.text(*pos[v], symbol, ha="center", va="center", fontsize=6.8,
color=style[1], fontweight=style[2],
bbox=dict(boxstyle="round,pad=0.2", fc=style[0], ec="#888888", lw=0.7))
ax.set(xticks=range(len(TIMES)), xticklabels=TIMES, yticks=[0, 1],
yticklabels=["regulatory\n(genes)", "signaling\n(proteins)"],
xlim=(-0.6, len(TIMES) - 0.4), ylim=(-0.45, 1.45), xlabel="time after TGF-beta")
ax.set_title("TGF-beta receptor-to-effector shortest paths", fontsize=11)
ax.legend([plt.Line2D([], [], color=c, lw=2) for c in COLOURS.values()],
[k.replace("coupling_", "") for k in COLOURS], ncol=4, frameon=False,
loc="upper center", bbox_to_anchor=(0.5, -0.11), fontsize=9)
for spine in ("top", "right", "left"):
ax.spines[spine].set_visible(False)
fig.tight_layout()
fig.savefig(FIGS / "signal_flow_paths.png", dpi=150)
plt.show()
The result goes back onto the graph¶
igraph computed it; AnnNet stores it. The arrival layer and hop count are written as ordinary vertex attributes, so they survive the snapshot and are available to every later section — the round-trip at the end checks exactly that.
In [17]:
Copied!
flow_attrs = {}
for symbol in {n[0] for n in supra if n[0].startswith(("prot:", "gene:"))}:
mechanism = "signaling" if symbol.startswith("prot:") else "regulatory"
hit = earliest(full, symbol.split(":", 1)[1], mechanism, symbol.split(":", 1)[0] + ":")
if hit:
flow_attrs[symbol] = {"tgfb_first_layer": hit[0], "tgfb_hops": int(hit[1])}
G.attrs.set_node_attrs_bulk(flow_attrs)
G.history.snapshot("after_signal_flow")
print(f"{len(flow_attrs):,} entities carry a TGF-β arrival time")
print(
G.views.nodes()
.filter(pl.col("tgfb_hops").is_not_null())
.group_by("tgfb_first_layer")
.len()
.sort("tgfb_first_layer")
.to_pandas()
.to_string(index=False)
)
flow_attrs = {}
for symbol in {n[0] for n in supra if n[0].startswith(("prot:", "gene:"))}:
mechanism = "signaling" if symbol.startswith("prot:") else "regulatory"
hit = earliest(full, symbol.split(":", 1)[1], mechanism, symbol.split(":", 1)[0] + ":")
if hit:
flow_attrs[symbol] = {"tgfb_first_layer": hit[0], "tgfb_hops": int(hit[1])}
G.attrs.set_node_attrs_bulk(flow_attrs)
G.history.snapshot("after_signal_flow")
print(f"{len(flow_attrs):,} entities carry a TGF-β arrival time")
print(
G.views.nodes()
.filter(pl.col("tgfb_hops").is_not_null())
.group_by("tgfb_first_layer")
.len()
.sort("tgfb_first_layer")
.to_pandas()
.to_string(index=False)
)
4,944 entities carry a TGF-β arrival time
tgfb_first_layer len
12h 1017
1h 4
24h 825
48h 978
72h 594
96h 1526
Part 5 — Cytoscape handoff¶
The Cytoscape view uses the receptor-to-effector path union from Part 4.
Cytoscape path subnetwork¶
In [18]:
Copied!
cx2 = an.to_cx2(G, path=str(OUT / "uc2.cx2"), export_name="UC2 response graph", hyperedges="skip")
print(f"whole graph CX2 : {len(cx2[4]['nodes']):,} nodes")
path_demo = an.AnnNet(directed=True)
path_label, path_vertices = {}, []
for v in sorted(nodes):
vid, (mechanism, layer) = names[v]
symbol = vid.split(":", 1)[1]
label = f"{symbol}|{mechanism}|{layer}"
role = "receptor" if v in seeds else "effector" if v in endpoints else "intermediate"
path_label[v] = label
path_vertices.append({"node_id": label, "gene_symbol": symbol, "mechanism": mechanism,
"time": layer, "role": role})
path_demo.add_nodes(path_vertices)
path_demo.add_edges(
[
{"edge_id": f"path:{i}", "source": path_label[a], "target": path_label[b],
"edge_kind": arc_kind[(a, b)], "weight": 1.0}
for i, (a, b) in enumerate(sorted(arcs, key=lambda e: (path_label[e[0]], path_label[e[1]])))
],
default_edge_directed=True,
)
path_cx2 = an.to_cx2(path_demo, path=str(OUT / "tgfb_effector_paths.cx2"),
export_name="TGF-beta receptor-to-effector paths", hyperedges="skip")
print(f"path subnetwork : {len(path_cx2[4]['nodes']):,} nodes, {len(path_cx2[5]['edges']):,} edges")
an.show_cx2(path_demo, hyperedges="skip", inline=True,
export_name="TGF-beta receptor-to-effector paths")
cx2 = an.to_cx2(G, path=str(OUT / "uc2.cx2"), export_name="UC2 response graph", hyperedges="skip")
print(f"whole graph CX2 : {len(cx2[4]['nodes']):,} nodes")
path_demo = an.AnnNet(directed=True)
path_label, path_vertices = {}, []
for v in sorted(nodes):
vid, (mechanism, layer) = names[v]
symbol = vid.split(":", 1)[1]
label = f"{symbol}|{mechanism}|{layer}"
role = "receptor" if v in seeds else "effector" if v in endpoints else "intermediate"
path_label[v] = label
path_vertices.append({"node_id": label, "gene_symbol": symbol, "mechanism": mechanism,
"time": layer, "role": role})
path_demo.add_nodes(path_vertices)
path_demo.add_edges(
[
{"edge_id": f"path:{i}", "source": path_label[a], "target": path_label[b],
"edge_kind": arc_kind[(a, b)], "weight": 1.0}
for i, (a, b) in enumerate(sorted(arcs, key=lambda e: (path_label[e[0]], path_label[e[1]])))
],
default_edge_directed=True,
)
path_cx2 = an.to_cx2(path_demo, path=str(OUT / "tgfb_effector_paths.cx2"),
export_name="TGF-beta receptor-to-effector paths", hyperedges="skip")
print(f"path subnetwork : {len(path_cx2[4]['nodes']):,} nodes, {len(path_cx2[5]['edges']):,} edges")
an.show_cx2(path_demo, hyperedges="skip", inline=True,
export_name="TGF-beta receptor-to-effector paths")
whole graph CX2 : 19,580 nodes path subnetwork : 46 nodes, 45 edges
Out[18]:
TGF-beta receptor-to-effector paths —
46 nodes, 45 edges
Part 6 — Late-response forecasting¶
Using only layers through 72h, the model ranks regulatory genes that newly join the response at 96h.
In [19]:
Copied!
import torch
from sklearn.metrics import average_precision_score, roc_auc_score
from torch_geometric.nn import SAGEConv
import torch
from sklearn.metrics import average_precision_score, roc_auc_score
from torch_geometric.nn import SAGEConv
/home/l1boll/miniconda3/lib/python3.13/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html from .autonotebook import tqdm as notebook_tqdm
In [20]:
Copied!
INPUT_TIMES = TIMES[:-1]
TARGET_TIME = TIMES[-1]
HIDDEN, LR, EPOCHS, TEST_FRAC = 64, 1e-2, 60, 0.30
window = G.layers.layer_union(
[(m, t) for m in MECHANISMS for t in INPUT_TIMES],
include_inter=True, # translation coupling: two different vertex ids, so `inter`
include_coupling=True, # time coupling: one vertex id across layers, so `coupling`
)
inside = E[E.edge_id.isin(window["edges"]) & E.edge_kind.isin(FAMILIES)]
inside = inside[inside.src_time.isin(INPUT_TIMES) & inside.tgt_time.isin(INPUT_TIMES)]
print(f"layer_union: {len(window['nodes']):,} nodes, {len(window['edges']):,} edges;"
f" {len(window['edges']) - len(inside):,} dropped as crossing the 72h boundary")
node_index, family_edges = {}, {k: [] for k in FAMILIES}
for r in inside.itertuples():
src = node_index.setdefault((r.src, r.src_mech, r.src_time), len(node_index))
tgt = node_index.setdefault((r.tgt, r.tgt_mech, r.tgt_time), len(node_index))
family_edges[r.edge_kind].append((src, tgt))
position = {t: i / (len(INPUT_TIMES) - 1) for i, t in enumerate(INPUT_TIMES)}
features = torch.tensor(
[
[
float(mechanism == "signaling"),
float(mechanism == "regulatory"),
position[t],
float(vid.split(":", 1)[1] in tf_symbols),
float(logfc_at.get((vid.split(":", 1)[1], t), 0.0)),
]
for vid, mechanism, t in node_index
],
dtype=torch.float32,
)
# Each candidate gene is read out at its latest regulatory supra-node <= 72h.
latest_node = {}
for (vid, mechanism, t), i in node_index.items():
symbol = vid.removeprefix("gene:")
if mechanism == "regulatory" and (
symbol not in latest_node or position[t] > position[latest_node[symbol][1]]
):
latest_node[symbol] = (i, t)
candidates = sorted(latest_node)
y = np.array(
[
int(s in responsive_at[TARGET_TIME] and s not in responsive_at[INPUT_TIMES[-1]])
for s in candidates
]
)
readout = torch.tensor([latest_node[s][0] for s in candidates], dtype=torch.long)
print(f"supra-nodes {len(node_index):,} | candidate genes {len(candidates):,}")
print(f"{TARGET_TIME}-specific positives: {y.sum()} ({y.mean():.1%} base rate)")
print("edge families:", {k: len(v) for k, v in family_edges.items()})
INPUT_TIMES = TIMES[:-1]
TARGET_TIME = TIMES[-1]
HIDDEN, LR, EPOCHS, TEST_FRAC = 64, 1e-2, 60, 0.30
window = G.layers.layer_union(
[(m, t) for m in MECHANISMS for t in INPUT_TIMES],
include_inter=True, # translation coupling: two different vertex ids, so `inter`
include_coupling=True, # time coupling: one vertex id across layers, so `coupling`
)
inside = E[E.edge_id.isin(window["edges"]) & E.edge_kind.isin(FAMILIES)]
inside = inside[inside.src_time.isin(INPUT_TIMES) & inside.tgt_time.isin(INPUT_TIMES)]
print(f"layer_union: {len(window['nodes']):,} nodes, {len(window['edges']):,} edges;"
f" {len(window['edges']) - len(inside):,} dropped as crossing the 72h boundary")
node_index, family_edges = {}, {k: [] for k in FAMILIES}
for r in inside.itertuples():
src = node_index.setdefault((r.src, r.src_mech, r.src_time), len(node_index))
tgt = node_index.setdefault((r.tgt, r.tgt_mech, r.tgt_time), len(node_index))
family_edges[r.edge_kind].append((src, tgt))
position = {t: i / (len(INPUT_TIMES) - 1) for i, t in enumerate(INPUT_TIMES)}
features = torch.tensor(
[
[
float(mechanism == "signaling"),
float(mechanism == "regulatory"),
position[t],
float(vid.split(":", 1)[1] in tf_symbols),
float(logfc_at.get((vid.split(":", 1)[1], t), 0.0)),
]
for vid, mechanism, t in node_index
],
dtype=torch.float32,
)
# Each candidate gene is read out at its latest regulatory supra-node <= 72h.
latest_node = {}
for (vid, mechanism, t), i in node_index.items():
symbol = vid.removeprefix("gene:")
if mechanism == "regulatory" and (
symbol not in latest_node or position[t] > position[latest_node[symbol][1]]
):
latest_node[symbol] = (i, t)
candidates = sorted(latest_node)
y = np.array(
[
int(s in responsive_at[TARGET_TIME] and s not in responsive_at[INPUT_TIMES[-1]])
for s in candidates
]
)
readout = torch.tensor([latest_node[s][0] for s in candidates], dtype=torch.long)
print(f"supra-nodes {len(node_index):,} | candidate genes {len(candidates):,}")
print(f"{TARGET_TIME}-specific positives: {y.sum()} ({y.mean():.1%} base rate)")
print("edge families:", {k: len(v) for k, v in family_edges.items()})
layer_union: 4,013 nodes, 31,655 edges; 2,814 dropped as crossing the 72h boundary
supra-nodes 9,301 | candidate genes 1,543
96h-specific positives: 108 (7.0% base rate)
edge families: {'signaling': 15308, 'regulatory': 6570, 'coupling_translation': 2133, 'coupling_time': 4830}
Model and priority table¶
In [21]:
Copied!
UNDIRECTED = {"signaling", "regulatory", "coupling_translation"}
def edge_index(families):
forward = [e for k in families for e in family_edges[k]]
index = torch.tensor(forward, dtype=torch.long).t().contiguous()
reverse = [e for k in families if k in UNDIRECTED for e in family_edges[k]]
return torch.cat([index, torch.tensor(reverse, dtype=torch.long).t().flip(0)], dim=1)
class Forecaster(torch.nn.Module):
def __init__(self, in_dim, hidden):
super().__init__()
self.conv1 = SAGEConv(in_dim, hidden)
self.conv2 = SAGEConv(hidden, hidden)
self.head = torch.nn.Linear(hidden, 1)
def forward(self, x, index):
x = self.conv1(x, index).relu()
return self.head(self.conv2(x, index).relu()).squeeze(-1)
rng = np.random.default_rng(SEED)
test_idx = rng.choice(len(candidates), int(TEST_FRAC * len(candidates)), replace=False)
train_idx = np.setdiff1d(np.arange(len(candidates)), test_idx)
target = torch.tensor(y, dtype=torch.float32)
index = edge_index(FAMILIES)
model = Forecaster(features.size(1), HIDDEN)
optimizer = torch.optim.Adam(model.parameters(), lr=LR)
for _ in range(EPOCHS):
optimizer.zero_grad()
loss = torch.nn.functional.binary_cross_entropy_with_logits(model(features, index)[readout[train_idx]], target[train_idx])
loss.backward()
optimizer.step()
with torch.no_grad():
scores = torch.sigmoid(model(features, index)[readout]).numpy()
print(f"AUROC={roc_auc_score(y[test_idx], scores[test_idx]):.3f} "
f"AP={average_precision_score(y[test_idx], scores[test_idx]):.3f}")
priority_rows = []
for gene, score, label in zip(candidates, scores, y.astype(bool)):
hit = flow_attrs.get(f"gene:{gene}", {})
priority_rows.append({"gene": gene, "score": score, "is_new_96h_responder": label,
"latest_input_layer": latest_node[gene][1],
"tgfb_first_layer": hit.get("tgfb_first_layer"),
"tgfb_hops": hit.get("tgfb_hops")})
priority = pd.DataFrame(priority_rows).iloc[test_idx].sort_values("score", ascending=False)
priority.to_csv(TABLES / "late_response_priorities.csv", index=False)
print(priority.head(15).to_string(index=False))
UNDIRECTED = {"signaling", "regulatory", "coupling_translation"}
def edge_index(families):
forward = [e for k in families for e in family_edges[k]]
index = torch.tensor(forward, dtype=torch.long).t().contiguous()
reverse = [e for k in families if k in UNDIRECTED for e in family_edges[k]]
return torch.cat([index, torch.tensor(reverse, dtype=torch.long).t().flip(0)], dim=1)
class Forecaster(torch.nn.Module):
def __init__(self, in_dim, hidden):
super().__init__()
self.conv1 = SAGEConv(in_dim, hidden)
self.conv2 = SAGEConv(hidden, hidden)
self.head = torch.nn.Linear(hidden, 1)
def forward(self, x, index):
x = self.conv1(x, index).relu()
return self.head(self.conv2(x, index).relu()).squeeze(-1)
rng = np.random.default_rng(SEED)
test_idx = rng.choice(len(candidates), int(TEST_FRAC * len(candidates)), replace=False)
train_idx = np.setdiff1d(np.arange(len(candidates)), test_idx)
target = torch.tensor(y, dtype=torch.float32)
index = edge_index(FAMILIES)
model = Forecaster(features.size(1), HIDDEN)
optimizer = torch.optim.Adam(model.parameters(), lr=LR)
for _ in range(EPOCHS):
optimizer.zero_grad()
loss = torch.nn.functional.binary_cross_entropy_with_logits(model(features, index)[readout[train_idx]], target[train_idx])
loss.backward()
optimizer.step()
with torch.no_grad():
scores = torch.sigmoid(model(features, index)[readout]).numpy()
print(f"AUROC={roc_auc_score(y[test_idx], scores[test_idx]):.3f} "
f"AP={average_precision_score(y[test_idx], scores[test_idx]):.3f}")
priority_rows = []
for gene, score, label in zip(candidates, scores, y.astype(bool)):
hit = flow_attrs.get(f"gene:{gene}", {})
priority_rows.append({"gene": gene, "score": score, "is_new_96h_responder": label,
"latest_input_layer": latest_node[gene][1],
"tgfb_first_layer": hit.get("tgfb_first_layer"),
"tgfb_hops": hit.get("tgfb_hops")})
priority = pd.DataFrame(priority_rows).iloc[test_idx].sort_values("score", ascending=False)
priority.to_csv(TABLES / "late_response_priorities.csv", index=False)
print(priority.head(15).to_string(index=False))
AUROC=0.956 AP=0.475 gene score is_new_96h_responder latest_input_layer tgfb_first_layer tgfb_hops BACH1 0.474527 True 24h 24h 7.0 NCOA3 0.467743 True 48h 48h 8.0 MXI1 0.467102 False 24h None NaN HMGCR 0.460829 False 24h 24h 7.0 DFFA 0.459496 False 24h 24h 6.0 CD63 0.459167 False 24h 24h 6.0 OGG1 0.454197 False 24h 24h 7.0 RIMS2 0.452184 True 48h 48h 7.0 PSMB10 0.451155 False 48h 48h 9.0 GULP1 0.450894 False 48h 48h 7.0 ST7 0.450444 False 24h 24h 6.0 PLAC8 0.444492 True 48h 24h 7.0 PTEN 0.440350 False 24h 24h 6.0 RIGI 0.434931 False 48h 48h 7.0 EXOC6B 0.434765 True 48h 24h 6.0
Round-trip¶
The final .annnet file preserves the graph, slices, hyperedges, arrival
attributes and build history.
In [22]:
Copied!
G.history.snapshot("final")
G.write(str(SNAPSHOT), overwrite=True)
reloaded = an.AnnNet.read(str(SNAPSHOT))
print(f"{SNAPSHOT.name} ({SNAPSHOT.stat().st_size / 1e6:.1f} MB)")
print(f"|V|={reloaded.global_count('nodes'):,} |E|={reloaded.global_count('edges'):,}"
f" aspects={reloaded.layers.list_aspects()} slices={len(reloaded.slices.list()):,}")
print(f"hyperedges : {len(reloaded.hyperedge_definitions):,}")
back = reloaded.views.nodes().filter(pl.col("tgfb_hops").is_not_null())
print(f"arrival times survived: {len(back):,} / {len(flow_attrs):,}")
print(f"attribute columns: {sorted(reloaded.views.nodes().columns)}")
G.history.snapshot("final")
G.write(str(SNAPSHOT), overwrite=True)
reloaded = an.AnnNet.read(str(SNAPSHOT))
print(f"{SNAPSHOT.name} ({SNAPSHOT.stat().st_size / 1e6:.1f} MB)")
print(f"|V|={reloaded.global_count('nodes'):,} |E|={reloaded.global_count('edges'):,}"
f" aspects={reloaded.layers.list_aspects()} slices={len(reloaded.slices.list()):,}")
print(f"hyperedges : {len(reloaded.hyperedge_definitions):,}")
back = reloaded.views.nodes().filter(pl.col("tgfb_hops").is_not_null())
print(f"arrival times survived: {len(back):,} / {len(flow_attrs):,}")
print(f"attribute columns: {sorted(reloaded.views.nodes().columns)}")
uc2.annnet (3.7 MB)
|V|=19,580 |E|=63,729 aspects=('mechanism', 'time') slices=20
hyperedges : 16,417
arrival times survived: 4,944 / 4,944
attribute columns: ['boundary_condition', 'compartment', 'constant', 'gene_symbol', 'has_only_substance_units', 'initial_concentration', 'kind', 'meta_id', 'name', 'node_id', 'sbo_term', 'tgfb_first_layer', 'tgfb_hops']
Scope¶
This is a software case study on a real biological substrate. Shortest paths are routes to inspect, not mechanisms by themselves. GNN scores are late-response priorities, not causal claims.
The 0-hop responsive restriction keeps time informative: a node-layer exists only where the entity is differential, and an intra-layer edge only where both endpoints are differential. Translation coupling is undirected because a gene and its protein are two readings of the same entity; time coupling is directed because time has an order.
In [ ]:
Copied!