From d6f163ed8130fd743c9241eb1357b2b2e67e5a36 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 20 Sep 2026 01:22:56 +0200 Subject: [PATCH] feat: summarise enzyme usage by system Capacity usage per enzyme and condition, the annotated systems drawn one panel per pathway, and the numbers written beside the figure. The transformation from per-enzyme usage to capacity usage per system is pure arithmetic on a table, so it is checked against the published output: fed the committed enzymeUsages.txt it reproduces the committed capUsage.txt exactly, to the row order. --- README.md | 10 ++ src/overflow/analyze_usage.py | 103 ++++++++++++++++++++ src/overflow/plots.py | 108 +++++++++++++++++++++ src/overflow/usage.py | 119 +++++++++++++++++++++++ tests/test_usage.py | 174 ++++++++++++++++++++++++++++++++++ 5 files changed, 514 insertions(+) create mode 100644 src/overflow/analyze_usage.py create mode 100644 src/overflow/plots.py create mode 100644 src/overflow/usage.py create mode 100644 tests/test_usage.py diff --git a/README.md b/README.md index 1b933ae..49a72d8 100644 --- a/README.md +++ b/README.md @@ -100,6 +100,16 @@ ribosomal subunits, at 10.5 amino acids per second per ribosome. The core is the 48 subunits whose average abundance across all conditions reaches 1e-5 mmol/gDW. +## Summarising enzyme usage + +```bash +python -m overflow.analyze_usage +``` + +Writes, per enzyme and condition, how much of it the model uses and what +fraction of what was available that is, plus the capacity usage of the +annotated systems and the two figures over them. + ## Tests ```bash diff --git a/src/overflow/analyze_usage.py b/src/overflow/analyze_usage.py new file mode 100644 index 0000000..0481717 --- /dev/null +++ b/src/overflow/analyze_usage.py @@ -0,0 +1,103 @@ +"""Summarising enzyme usage across the condition models. + +Run after the models are built: ``python -m overflow.analyze_usage``. +By default it reads the models with the ribosome in them. +""" +from __future__ import annotations + +import argparse +from pathlib import Path +from typing import Optional + +import pandas as pd +from geckopy import load_ec_model + +from overflow.adapter import build_adapter +from overflow.config import ( + BIO_RXN, + CONDITION_ORDER, + MODELS_DIR, + POOL_RXN, + RESULTS_DIR, + load_conditions, +) +from overflow.plots import ( + SELECTED_SYSTEMS, + SUPPLEMENT_SYSTEMS, + capacity_usage_figure, +) +from overflow.usage import ( + capacity_usage_by_system, + combine_usage, + enzyme_usage_table, + read_annotation, + usage_by_system, +) + + +def analyse( + conditions: Optional[list[str]] = None, + models_dir: Path = MODELS_DIR, + suffix: str = "_ribosome", + verbose: bool = True, +) -> dict[str, pd.DataFrame]: + """Per-enzyme usage for each condition, from its model's own optimum.""" + all_conditions = load_conditions() + tables: dict[str, pd.DataFrame] = {} + + for name in conditions or CONDITION_ORDER: + condition = all_conditions[name] + path = models_dir / f"ecModel_P_{name}{suffix}.yml" + model = load_ec_model(str(path), adapter=build_adapter(condition)) + solution = model.optimize() + if solution.status != "optimal": + raise RuntimeError(f"{name}: {path.name} has no solution ({solution.status})") + tables[name] = enzyme_usage_table(model, solution.fluxes) + if verbose: + used = tables[name]["absUse"] > 0 + print( + f"[{name}] growth {solution.fluxes[BIO_RXN]:.5f}, protein " + f"{solution.fluxes[POOL_RXN]:.1f} mg/gDW, " + f"{int(used.sum())} of {len(used)} enzymes carrying usage", + flush=True, + ) + return tables + + +def main(argv: Optional[list[str]] = None) -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("conditions", nargs="*", default=list(CONDITION_ORDER)) + parser.add_argument("--models-dir", type=Path, default=MODELS_DIR) + parser.add_argument("--results-dir", type=Path, default=RESULTS_DIR) + parser.add_argument("--suffix", default="_ribosome", + help="model file suffix, '' for the models without a ribosome") + parser.add_argument("--solver") + args = parser.parse_args(argv) + + if args.solver: + import cobra + + cobra.Configuration().solver = args.solver + + tables = analyse(args.conditions, args.models_dir, args.suffix) + usage = combine_usage(tables) + capacity = capacity_usage_by_system(usage, read_annotation()) + + out = args.results_dir / "enzymeUsage" + out.mkdir(parents=True, exist_ok=True) + usage.to_csv(out / "enzymeUsages.tsv", sep="\t", index=False) + capacity.to_csv(out / "capUsage.tsv", sep="\t", index=False) + medians = usage_by_system(capacity) + medians.to_csv(out / "systemMedians.tsv", sep="\t", index=False) + + capacity_usage_figure(capacity, SELECTED_SYSTEMS, out / "selectedSystemUsage.pdf") + capacity_usage_figure(capacity, SUPPLEMENT_SYSTEMS, out / "supplementSystemUsage.pdf") + + print() + print("median capacity usage per system (%)") + print(medians.to_string(index=False)) + return 0 + + +if __name__ == "__main__": # pragma: no cover + raise SystemExit(main()) diff --git a/src/overflow/plots.py b/src/overflow/plots.py new file mode 100644 index 0000000..1d7d76e --- /dev/null +++ b/src/overflow/plots.py @@ -0,0 +1,108 @@ +"""Figures summarising enzyme capacity usage. + +Each panel is one pathway and each box one condition, so a pathway that +tightens as the carbon-to-nitrogen ratio falls shows up as boxes walking +upward. Colour repeats the panel title rather than carrying identity of +its own; the numbers behind every panel are written beside the figure as +a table. +""" +from __future__ import annotations + +from pathlib import Path +from typing import Optional, Sequence + +import matplotlib +import pandas as pd + +matplotlib.use("Agg") +import matplotlib.pyplot as plt # noqa: E402 + +from overflow.config import CONDITION_ORDER # noqa: E402 + +#: Pathways shown in the main figure and in the supplement. +SELECTED_SYSTEMS = ("Glycolysis", "TCA cycle", "ETC", "Ribosome") +SUPPLEMENT_SYSTEMS = ( + "PP shunt", + "THF cycle", + "Nitrogen metabolism", + "Amino acid metabolism", +) + +#: Categorical slots one to four, validated against both surfaces. +PALETTE = ("#2a78d6", "#eb6834", "#1baf7a", "#eda100") + +INK = "#0b0b0b" +MUTED = "#52514e" +GRID = "#e6e5e0" +SURFACE = "#fcfcfb" + + +def capacity_usage_figure( + capacity: pd.DataFrame, + systems: Sequence[str] = SELECTED_SYSTEMS, + path: Optional[Path | str] = None, + conditions: Optional[Sequence[str]] = None, + width: float = 7.2, + height: float = 2.6, +): + """One panel per pathway, one box per condition. + + ``capacity`` is the table from + :func:`overflow.usage.capacity_usage_by_system`. + """ + if conditions is None: + conditions = [c for c in CONDITION_ORDER if c in capacity.columns] + + missing = [s for s in systems if s not in set(capacity["GOterm"])] + if missing: + raise ValueError(f"no enzymes annotated to {missing}") + + figure, axes = plt.subplots( + 1, len(systems), figsize=(width, height), sharey=True, + facecolor=SURFACE, layout="constrained", + ) + axes = [axes] if len(systems) == 1 else list(axes) + + for index, (axis, system) in enumerate(zip(axes, systems)): + rows = capacity[capacity["GOterm"] == system] + values = [rows[c].dropna().to_numpy() for c in conditions] + colour = PALETTE[index % len(PALETTE)] + + drawn = axis.boxplot( + values, + widths=0.55, + showfliers=True, + patch_artist=True, + medianprops=dict(color=INK, linewidth=1.1), + flierprops=dict( + marker="o", markersize=2.2, markerfacecolor="none", + markeredgecolor=colour, markeredgewidth=0.6, + ), + ) + for box in drawn["boxes"]: + box.set(facecolor=colour, alpha=0.18, edgecolor=colour, linewidth=0.9) + for part in ("whiskers", "caps"): + for line in drawn[part]: + line.set(color=colour, linewidth=0.9) + + axis.set_title(system, fontsize=8, color=INK, pad=6) + axis.set_xticks(range(1, len(conditions) + 1)) + axis.set_xticklabels(conditions, rotation=90, fontsize=7, color=MUTED) + axis.tick_params(axis="y", labelsize=7, colors=MUTED, length=3) + axis.tick_params(axis="x", length=0) + axis.set_facecolor(SURFACE) + axis.set_axisbelow(True) + axis.yaxis.grid(True, color=GRID, linewidth=0.5) + for side in ("top", "right"): + axis.spines[side].set_visible(False) + for side in ("left", "bottom"): + axis.spines[side].set(color=GRID, linewidth=0.6) + if index == 0: + axis.set_ylabel("Capacity usage (%)", fontsize=8, color=INK) + + axes[0].set_ylim(-4, 104) + + if path is not None: + figure.savefig(path, facecolor=SURFACE) + plt.close(figure) + return figure diff --git a/src/overflow/usage.py b/src/overflow/usage.py new file mode 100644 index 0000000..fdad040 --- /dev/null +++ b/src/overflow/usage.py @@ -0,0 +1,119 @@ +"""Summarising how much of each enzyme the model uses. + +Two quantities per enzyme: how much of it the flux distribution needs +(absolute usage, mg/gDW) and how much of what was available that is +(capacity usage). An enzyme at full capacity is one the condition is +pressing against; one well below it is being carried. +""" +from __future__ import annotations + +from pathlib import Path +from typing import TYPE_CHECKING, Mapping, Optional, Sequence + +import numpy as np +import pandas as pd +from geckopy import enzyme_usage + +from overflow.config import ANNOTATION_DATA, CONDITION_ORDER + +if TYPE_CHECKING: # pragma: no cover - typing only + from geckopy import EcModel + +ID_COLUMNS = ["protID", "geneID", "protName"] + + +def enzyme_usage_table(model: "EcModel", fluxes: Mapping[str, float]) -> pd.DataFrame: + """Per-enzyme usage for one condition.""" + result = enzyme_usage(model, fluxes) + gene_of = dict(zip(model.ec.enzymes, model.ec.genes)) + + def name_of(gene: str) -> str: + try: + return model.genes.get_by_id(gene).name or gene + except KeyError: + return gene + + genes = [gene_of.get(p, "") for p in result.prot_id] + return pd.DataFrame( + { + "protID": list(result.prot_id), + "geneID": genes, + "protName": [name_of(g) for g in genes], + "capUse": np.asarray(result.cap_usage, dtype=float), + "absUse": np.asarray(result.abs_usage, dtype=float), + "UB": np.asarray(result.ub, dtype=float), + } + ) + + +def combine_usage(tables: Mapping[str, pd.DataFrame]) -> pd.DataFrame: + """One wide table over conditions, in the published column layout.""" + conditions = [c for c in CONDITION_ORDER if c in tables] + base = tables[conditions[0]][ID_COLUMNS].copy() + for quantity in ("capUse", "absUse", "UB"): + for condition in conditions: + table = tables[condition].set_index("protID")[quantity] + base[f"{quantity}_{condition}"] = base["protID"].map(table) + return base + + +def read_annotation(path: Path | str = ANNOTATION_DATA) -> pd.DataFrame: + """Read the assignment of proteins to systems.""" + table = pd.read_csv(path, sep="\t") + return table.rename( + columns={ + "Entry": "protID", + "Gene names (ordered locus )": "geneID", + "Gene names (primary )": "protName", + "system": "system", + } + )[["protID", "geneID", "protName", "system"]] + + +def capacity_usage_by_system( + usage: pd.DataFrame, + annotation: Optional[pd.DataFrame] = None, + conditions: Optional[Sequence[str]] = None, + decimals: int = 3, +) -> pd.DataFrame: + """Capacity usage as a percentage, for the annotated systems only. + + Enzymes the model never uses in any condition are dropped: a + capacity usage of zero everywhere says nothing about how the + condition allocates protein. + """ + if annotation is None: + annotation = read_annotation() + if conditions is None: + conditions = [ + c for c in CONDITION_ORDER if f"capUse_{c}" in usage.columns + ] + + columns = [f"capUse_{c}" for c in conditions] + table = usage[ID_COLUMNS + columns].copy() + table = table[table[columns].sum(axis=1) != 0] + table[columns] = table[columns] * 100 + + # A protein listed under more than one system keeps the first, as + # the published summary does. + systems = annotation.drop_duplicates(subset="protID", keep="first").set_index( + "protID" + )["system"] + table = table[table["protID"].isin(systems.index)] + table["GOterm"] = table["protID"].map(systems) + + table = table.rename(columns={f"capUse_{c}": c for c in conditions}) + table[list(conditions)] = table[list(conditions)].round(decimals) + return table.reset_index(drop=True) + + +def usage_by_system(capacity: pd.DataFrame, conditions: Optional[Sequence[str]] = None) -> pd.DataFrame: + """Median capacity usage per system and condition.""" + if conditions is None: + conditions = [c for c in CONDITION_ORDER if c in capacity.columns] + return ( + capacity.groupby("GOterm")[list(conditions)] + .median() + .round(2) + .reset_index() + ) diff --git a/tests/test_usage.py b/tests/test_usage.py new file mode 100644 index 0000000..ff7114e --- /dev/null +++ b/tests/test_usage.py @@ -0,0 +1,174 @@ +"""Summarising enzyme usage. + +The transformation from per-enzyme usage to capacity usage per system is +pure arithmetic on a table, so it can be checked against the published +output exactly: feed it the committed `enzymeUsages.txt` and it must +reproduce the committed `capUsage.txt`. +""" +import numpy as np +import pandas as pd +import pytest + +from overflow.config import ROOT +from overflow.usage import ( + capacity_usage_by_system, + combine_usage, + read_annotation, + usage_by_system, +) + +LEGACY = ROOT / "legacy_matlab" / "results" / "enzymeUsage" + + +@pytest.fixture(scope="module") +def legacy_usage(): + return pd.read_csv(LEGACY / "enzymeUsages.txt", sep="\t") + + +@pytest.fixture(scope="module") +def legacy_capacity(): + return pd.read_csv(LEGACY / "capUsage.txt", sep="\t") + + +@pytest.fixture(scope="module") +def annotation(): + return read_annotation() + + +def test_capacity_usage_reproduces_the_published_table( + legacy_usage, legacy_capacity, annotation +): + ours = capacity_usage_by_system(legacy_usage, annotation) + assert list(ours.columns) == list(legacy_capacity.columns) + assert list(ours["protID"]) == list(legacy_capacity["protID"]) + assert list(ours["GOterm"]) == list(legacy_capacity["GOterm"]) + for condition in ("CN4", "CN22", "CN38", "CN75", "hGR"): + assert ours[condition].to_numpy() == pytest.approx( + legacy_capacity[condition].to_numpy() + ) + + +def test_enzymes_used_in_no_condition_are_dropped(annotation): + usage = pd.DataFrame( + { + "protID": ["P00127", "P00128"], + "geneID": ["YFR033C", "YDR529C"], + "protName": ["QCR6", "QCR7"], + "capUse_CN4": [0.0, 0.5], + "capUse_CN22": [0.0, 0.25], + } + ) + result = capacity_usage_by_system(usage, annotation) + assert list(result["protID"]) == ["P00128"] + + +def test_capacity_usage_is_reported_as_a_percentage(annotation): + usage = pd.DataFrame( + { + "protID": ["P00127"], "geneID": ["YFR033C"], "protName": ["QCR6"], + "capUse_CN4": [0.4213], + } + ) + assert capacity_usage_by_system(usage, annotation)["CN4"].iloc[0] == pytest.approx(42.13) + + +def test_an_unannotated_protein_is_left_out(annotation): + usage = pd.DataFrame( + { + "protID": ["NOT_ANNOTATED"], "geneID": ["YXX000W"], "protName": ["X"], + "capUse_CN4": [0.5], + } + ) + assert capacity_usage_by_system(usage, annotation).empty + + +def test_a_protein_in_two_systems_keeps_the_first(annotation): + """ADE3 is listed under both amino acid metabolism and the THF cycle.""" + duplicated = annotation[annotation["protID"] == "P07245"] + assert len(duplicated) == 2 + usage = pd.DataFrame( + { + "protID": ["P07245"], "geneID": ["YGR204W"], "protName": ["ADE3"], + "capUse_CN4": [0.5], + } + ) + result = capacity_usage_by_system(usage, annotation) + assert result["GOterm"].iloc[0] == duplicated["system"].iloc[0] + + +def test_combining_conditions_keeps_the_published_column_layout(): + def one(cap, abs_, ub): + return pd.DataFrame( + { + "protID": ["P1"], "geneID": ["Y1"], "protName": ["N1"], + "capUse": [cap], "absUse": [abs_], "UB": [ub], + } + ) + + wide = combine_usage({"CN4": one(0.5, 1.0, 2.0), "hGR": one(0.25, 1.0, 4.0)}) + assert list(wide.columns) == [ + "protID", "geneID", "protName", + "capUse_CN4", "capUse_hGR", "absUse_CN4", "absUse_hGR", "UB_CN4", "UB_hGR", + ] + assert wide["capUse_hGR"].iloc[0] == 0.25 + + +def test_combining_conditions_matches_proteins_by_identifier(): + """The conditions need not list their enzymes in the same order.""" + first = pd.DataFrame( + {"protID": ["P1", "P2"], "geneID": ["Y1", "Y2"], "protName": ["A", "B"], + "capUse": [0.1, 0.2], "absUse": [1.0, 2.0], "UB": [10.0, 10.0]} + ) + second = first.iloc[::-1].reset_index(drop=True).copy() + second["capUse"] = [0.9, 0.8] + wide = combine_usage({"CN4": first, "hGR": second}) + row = wide.set_index("protID").loc["P1"] + assert row["capUse_CN4"] == 0.1 and row["capUse_hGR"] == 0.8 + + +def test_system_medians_summarise_the_table(legacy_usage, annotation): + capacity = capacity_usage_by_system(legacy_usage, annotation) + medians = usage_by_system(capacity) + assert set(medians["GOterm"]) == set(capacity["GOterm"]) + assert (medians[["CN4", "hGR"]] <= 100.0).all().all() + + +# --- figures ---------------------------------------------------------- + +def test_the_figure_draws_one_panel_per_system(legacy_usage, annotation, tmp_path): + from overflow.plots import SELECTED_SYSTEMS, capacity_usage_figure + + capacity = capacity_usage_by_system(legacy_usage, annotation) + figure = capacity_usage_figure(capacity, SELECTED_SYSTEMS, tmp_path / "f.pdf") + assert [axis.get_title() for axis in figure.axes] == list(SELECTED_SYSTEMS) + assert (tmp_path / "f.pdf").stat().st_size > 1000 + + +def test_every_panel_actually_has_boxes(legacy_usage, annotation, tmp_path): + """A panel whose system matched nothing would still draw, empty.""" + from overflow.plots import SELECTED_SYSTEMS, capacity_usage_figure + + capacity = capacity_usage_by_system(legacy_usage, annotation) + figure = capacity_usage_figure(capacity, SELECTED_SYSTEMS, tmp_path / "f.pdf") + for axis, system in zip(figure.axes, SELECTED_SYSTEMS): + assert len(axis.patches) == 5, f"{system} should have one box per condition" + + +def test_a_system_with_no_enzymes_is_refused(legacy_usage, annotation, tmp_path): + from overflow.plots import capacity_usage_figure + + capacity = capacity_usage_by_system(legacy_usage, annotation) + with pytest.raises(ValueError, match="no enzymes annotated"): + capacity_usage_figure(capacity, ["Photosynthesis"], tmp_path / "f.pdf") + + +def test_the_axis_covers_the_full_percentage_range(legacy_usage, annotation, tmp_path): + """Capacity usage is a fraction of what was available, so panels have + to share a 0-100 axis or they cannot be compared.""" + from overflow.plots import SELECTED_SYSTEMS, capacity_usage_figure + + capacity = capacity_usage_by_system(legacy_usage, annotation) + figure = capacity_usage_figure(capacity, SELECTED_SYSTEMS, tmp_path / "f.pdf") + for axis in figure.axes: + low, high = axis.get_ylim() + assert low <= 0 and high >= 100