Phase 5: enzyme usage by system

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

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())

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

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()
+    )