From 21ccad1f25181dce8c1f6c513a8eb8c76c726511 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 20 Sep 2026 13:37:44 +0200 Subject: [PATCH 1/3] feat: compare the results against the published ones Regenerates every side-by-side table from the committed MATLAB outputs, so the comparison is reproducible rather than quoted. Uptakes are compared unsigned: GECKO 2 splits reversible reactions, so its glucose uptake is a positive flux on r_1714_REV where this pipeline has a negative flux on r_1714. --- README.md | 10 +++ src/overflow/compare.py | 177 ++++++++++++++++++++++++++++++++++++++++ tests/test_compare.py | 95 +++++++++++++++++++++ 3 files changed, 282 insertions(+) create mode 100644 src/overflow/compare.py create mode 100644 tests/test_compare.py diff --git a/README.md b/README.md index 0fa5161..b65fae2 100644 --- a/README.md +++ b/README.md @@ -127,6 +127,16 @@ means rather than an average over the interior of the flux space. `--no-min-flux skips the second step, which lets the sampler wander into high-flux routes: the pentose phosphate pathway then carries several times the published flux. +## Comparing against the published results + +```bash +python -m overflow.compare +``` + +Writes `results/COMPARISON.md`: predicted exchange rates and median capacity +usage per system, this pipeline beside the MATLAB one, read from +`legacy_matlab/results/`. + ## Tests ```bash diff --git a/src/overflow/compare.py b/src/overflow/compare.py new file mode 100644 index 0000000..d0a9af6 --- /dev/null +++ b/src/overflow/compare.py @@ -0,0 +1,177 @@ +"""Comparing this pipeline's results with the published ones. + +The MATLAB implementation's outputs are committed under +``legacy_matlab/results``. Where a quantity means the same thing in both, +it is put side by side here rather than described in prose. +""" +from __future__ import annotations + +import argparse +from pathlib import Path +from typing import Optional + +import numpy as np +import pandas as pd + +from overflow.config import CONDITION_ORDER, RESULTS_DIR, ROOT, load_conditions +from overflow.ngam import relative_residual +from overflow.usage import usage_by_system + +LEGACY = ROOT / "legacy_matlab" / "results" + +#: Exchange reactions as the two model formats name them. GECKO 2 split +#: every reversible reaction, so an uptake is a separate reaction there. +EXCHANGES = { + "glucose": ("r_1714_REV", "r_1714"), + "CO2": ("r_1672", "r_1672"), + "oxygen": ("r_1992_REV", "r_1992"), + "ethanol": ("r_1761", "r_1761"), +} + + +def _legacy_fluxes(condition: str) -> dict[str, float]: + fluxes: dict[str, float] = {} + path = LEGACY / "modelSimulation" / f"allFluxes_{condition}.txt" + for line in path.read_text().splitlines(): + parts = line.split("\t") + if len(parts) < 4: + continue + try: + fluxes[parts[0]] = float(parts[3]) + except ValueError: + continue + return fluxes + + +def _new_fluxes(condition: str, results_dir: Path) -> dict[str, float]: + table = pd.read_csv( + results_dir / "modelSimulation" / f"allFluxes_{condition}.tsv", sep="\t" + ) + return dict(zip(table["rxnID"], table["flux"])) + + +def exchange_comparison(results_dir: Path = RESULTS_DIR) -> pd.DataFrame: + """Predicted exchange rates and how far each is from the measurement.""" + conditions = load_conditions() + rows = [] + for name in CONDITION_ORDER: + condition = conditions[name] + legacy, new = _legacy_fluxes(name), _new_fluxes(name, results_dir) + measured = np.array([condition.glucose, condition.co2, condition.oxygen]) + legacy_rates = np.array([abs(legacy.get(EXCHANGES[k][0], 0.0)) + for k in ("glucose", "CO2", "oxygen")]) + new_rates = np.array([abs(new.get(EXCHANGES[k][1], 0.0)) + for k in ("glucose", "CO2", "oxygen")]) + rows.append( + { + "condition": name, + "glucose_measured": condition.glucose, + "glucose_gecko2": legacy_rates[0], "glucose_gecko4": new_rates[0], + "CO2_measured": condition.co2, + "CO2_gecko2": legacy_rates[1], "CO2_gecko4": new_rates[1], + "oxygen_measured": condition.oxygen, + "oxygen_gecko2": legacy_rates[2], "oxygen_gecko4": new_rates[2], + "ethanol_measured": condition.byproducts["Ethanol"], + "ethanol_gecko2": legacy.get("r_1761", 0.0), + "ethanol_gecko4": new.get("r_1761", 0.0), + "residual_gecko2": relative_residual(legacy_rates, measured), + "residual_gecko4": relative_residual(new_rates, measured), + } + ) + return pd.DataFrame(rows).round(4) + + +def usage_comparison(results_dir: Path = RESULTS_DIR) -> pd.DataFrame: + """Median capacity usage per system, side by side.""" + legacy = usage_by_system( + pd.read_csv(LEGACY / "enzymeUsage" / "capUsage.txt", sep="\t") + ).set_index("GOterm") + new = usage_by_system( + pd.read_csv(results_dir / "enzymeUsage" / "capUsage.tsv", sep="\t") + ).set_index("GOterm") + + rows = [] + for system in sorted(set(legacy.index) | set(new.index)): + row: dict[str, object] = {"system": system} + for condition in CONDITION_ORDER: + row[f"{condition}_gecko2"] = ( + legacy.loc[system, condition] if system in legacy.index else np.nan + ) + row[f"{condition}_gecko4"] = ( + new.loc[system, condition] if system in new.index else np.nan + ) + rows.append(row) + return pd.DataFrame(rows) + + +def report(results_dir: Path = RESULTS_DIR) -> str: + """A written comparison of the two implementations.""" + exchanges = exchange_comparison(results_dir) + lines = [ + "# This pipeline against the published one", + "", + "Generated by `python -m overflow.compare`. GECKO 2 values are read from", + "`legacy_matlab/results/`, which is what the MATLAB implementation wrote.", + "", + "## Exchange rates", + "", + "Measured / GECKO 2 / GECKO 4, mmol/gDW/h. The residual is the relative", + "deviation over glucose, CO2 and oxygen together.", + "", + "| condition | glucose | CO2 | oxygen | ethanol | residual G2 | residual G4 |", + "|---|---|---|---|---|---|---|", + ] + for row in exchanges.itertuples(): + lines.append( + f"| {row.condition} " + f"| {row.glucose_measured:.3f} / {row.glucose_gecko2:.3f} / {row.glucose_gecko4:.3f} " + f"| {row.CO2_measured:.3f} / {row.CO2_gecko2:.3f} / {row.CO2_gecko4:.3f} " + f"| {row.oxygen_measured:.3f} / {row.oxygen_gecko2:.3f} / {row.oxygen_gecko4:.3f} " + f"| {row.ethanol_measured:.3f} / {row.ethanol_gecko2:.3f} / {row.ethanol_gecko4:.3f} " + f"| {row.residual_gecko2:.3f} | {row.residual_gecko4:.3f} |" + ) + lines += [ + "", + f"Mean residual: GECKO 2 {exchanges.residual_gecko2.mean():.3f}, " + f"GECKO 4 {exchanges.residual_gecko4.mean():.3f}.", + "", + ] + + try: + usage = usage_comparison(results_dir) + except FileNotFoundError: + return "\n".join(lines) + + lines += [ + "## Median capacity usage per system (%)", + "", + "| system | " + " | ".join(CONDITION_ORDER) + " |", + "|---" * (len(CONDITION_ORDER) + 1) + "|", + ] + for row in usage.itertuples(index=False): + values = [] + for condition in CONDITION_ORDER: + old = getattr(row, f"{condition}_gecko2") + new = getattr(row, f"{condition}_gecko4") + values.append(f"{old:.1f} / {new:.1f}") + lines.append(f"| {row.system} | " + " | ".join(values) + " |") + lines.append("") + return "\n".join(lines) + + +def main(argv: Optional[list[str]] = None) -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--results-dir", type=Path, default=RESULTS_DIR) + parser.add_argument("--output", type=Path, default=None) + args = parser.parse_args(argv) + + text = report(args.results_dir) + destination = args.output or (args.results_dir / "COMPARISON.md") + destination.parent.mkdir(parents=True, exist_ok=True) + destination.write_text(text + "\n") + print(text) + return 0 + + +if __name__ == "__main__": # pragma: no cover + raise SystemExit(main()) diff --git a/tests/test_compare.py b/tests/test_compare.py new file mode 100644 index 0000000..95a0a06 --- /dev/null +++ b/tests/test_compare.py @@ -0,0 +1,95 @@ +"""The comparison against the published results.""" +import pandas as pd +import pytest + +from overflow.compare import LEGACY, exchange_comparison, report, usage_comparison +from overflow.config import CONDITION_ORDER, ROOT + + +@pytest.fixture(scope="module") +def results_dir(tmp_path_factory): + """A results directory holding the legacy outputs, reshaped as this + pipeline writes them. Lets the comparison be exercised without a run.""" + root = tmp_path_factory.mktemp("results") + (root / "modelSimulation").mkdir() + (root / "enzymeUsage").mkdir() + + # GECKO 2 splits reversible reactions, so an uptake lives on the + # _REV copy as a positive flux. Fold those back onto the base + # reaction as a negative flux, which is this pipeline's format. + reversed_uptakes = {"r_1714_REV": "r_1714", "r_1992_REV": "r_1992"} + + for condition in CONDITION_ORDER: + rows = {} + for line in (LEGACY / "modelSimulation" / f"allFluxes_{condition}.txt").read_text().splitlines(): + parts = line.split("\t") + if len(parts) < 4: + continue + try: + flux = float(parts[3]) + except ValueError: + continue + reaction = parts[0] + if reaction in reversed_uptakes: + rows[reversed_uptakes[reaction]] = { + "rxnID": reversed_uptakes[reaction], "rxnName": parts[1], + "equation": parts[2], "flux": -flux, + } + elif reaction not in rows: + rows[reaction] = { + "rxnID": reaction, "rxnName": parts[1], + "equation": parts[2], "flux": flux, + } + pd.DataFrame(list(rows.values())).to_csv( + root / "modelSimulation" / f"allFluxes_{condition}.tsv", sep="\t", index=False + ) + + pd.read_csv(LEGACY / "enzymeUsage" / "capUsage.txt", sep="\t").to_csv( + root / "enzymeUsage" / "capUsage.tsv", sep="\t", index=False + ) + return root + + +def test_the_comparison_covers_every_condition(results_dir): + table = exchange_comparison(results_dir) + assert list(table["condition"]) == list(CONDITION_ORDER) + + +def test_comparing_a_run_against_itself_gives_no_difference(results_dir): + """Fed the legacy fluxes as though they were new ones, the two + columns must agree -- otherwise the comparison is reading one of the + two formats wrongly.""" + table = exchange_comparison(results_dir) + for quantity in ("glucose", "CO2", "oxygen", "ethanol"): + assert table[f"{quantity}_gecko2"].to_numpy() == pytest.approx( + table[f"{quantity}_gecko4"].to_numpy(), abs=1e-6 + ) + assert table["residual_gecko2"].to_numpy() == pytest.approx( + table["residual_gecko4"].to_numpy() + ) + + +def test_uptakes_are_compared_as_magnitudes(results_dir): + """GECKO 2 splits reversible reactions, so its glucose uptake is a + positive flux on r_1714_REV while this pipeline's is a negative flux + on r_1714. Comparing them unsigned is the whole point.""" + table = exchange_comparison(results_dir) + assert (table["glucose_gecko2"] > 0).all() + assert (table["glucose_gecko4"] > 0).all() + + +def test_usage_comparison_lines_the_systems_up(results_dir): + table = usage_comparison(results_dir) + assert "ETC" in set(table["system"]) + for condition in CONDITION_ORDER: + assert table[f"{condition}_gecko2"].to_numpy() == pytest.approx( + table[f"{condition}_gecko4"].to_numpy(), nan_ok=True + ) + + +def test_the_report_names_what_it_compares(results_dir): + text = report(results_dir) + assert "Exchange rates" in text + assert "Mean residual" in text + for condition in CONDITION_ORDER: + assert condition in text From a534c5dddd230241df4482e5aec62f883ea6d329 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 20 Sep 2026 14:09:30 +0200 Subject: [PATCH 2/3] docs: cover rate fitting, the sampled budget and the run order --- README.md | 27 ++++++++++++++++++----- src/overflow/compare.py | 49 +++++++++++++++++++++++++++++++++++++++++ tests/test_compare.py | 21 ++++++++++++++++++ 3 files changed, 92 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index b65fae2..36f1e64 100644 --- a/README.md +++ b/README.md @@ -133,13 +133,30 @@ pentose phosphate pathway then carries several times the published flux. python -m overflow.compare ``` -Writes `results/COMPARISON.md`: predicted exchange rates and median capacity -usage per system, this pipeline beside the MATLAB one, read from -`legacy_matlab/results/`. +Writes `results/COMPARISON.md`: predicted exchange rates, median capacity usage +per system and the sampled ATP and redox budget, this pipeline beside the MATLAB +one, read from `legacy_matlab/results/`. + +## Running the whole analysis + +```bash +python -m overflow.build --fit-rates --rate-tolerance 0.08 --uptake-flex 1.08 +python -m overflow.build_ribosome +python -m overflow.analyze_usage +python -m overflow.run_sampling --procs 12 +python -m overflow.compare +``` + +Each step reads the models the previous one wrote, so `--models-dir` and +`--results-dir` keep a run self-contained. ## Tests ```bash -pytest # everything -pytest -m "not slow" # skip the tests that load a genome-scale model +pytest # everything except the full condition builds +pytest -m "not slow" # skip the tests that load a genome-scale model +pytest -m integration # build a condition end to end; minutes ``` + +The suite runs on GLPK by default, so it needs no commercial solver licence; +set `OVERFLOW_TEST_SOLVER` to use another. diff --git a/src/overflow/compare.py b/src/overflow/compare.py index d0a9af6..92350fa 100644 --- a/src/overflow/compare.py +++ b/src/overflow/compare.py @@ -104,6 +104,34 @@ def usage_comparison(results_dir: Path = RESULTS_DIR) -> pd.DataFrame: return pd.DataFrame(rows) +#: Rows of the sampling summary worth putting side by side. +BUDGET_ROWS = ( + "rGlu", "ETC_rATP", "glycolysis_rATP", "GAEC_rATP", "NGAM_rATP", + "GAEC+NGAM+Metabolism_rATP", "rPDH", "rIDH", "rMDHc", "rMDHm", "rNDE", +) + + +def budget_comparison(results_dir: Path = RESULTS_DIR) -> pd.DataFrame: + """The sampled ATP and redox budget, side by side.""" + legacy = pd.read_csv( + LEGACY / "randomSampling" / "selectedFluxes.txt", sep="\t" + ).set_index("Row") + new = pd.read_csv( + results_dir / "randomSampling" / "selectedFluxes.tsv", sep="\t" + ).set_index("Row") + + rows = [] + for name in BUDGET_ROWS: + if name not in legacy.index or name not in new.index: + continue + row: dict[str, object] = {"row": name} + for condition in CONDITION_ORDER: + row[f"{condition}_gecko2"] = float(legacy.loc[name, condition]) + row[f"{condition}_gecko4"] = float(new.loc[name, condition]) + rows.append(row) + return pd.DataFrame(rows) + + def report(results_dir: Path = RESULTS_DIR) -> str: """A written comparison of the two implementations.""" exchanges = exchange_comparison(results_dir) @@ -156,6 +184,27 @@ def report(results_dir: Path = RESULTS_DIR) -> str: values.append(f"{old:.1f} / {new:.1f}") lines.append(f"| {row.system} | " + " | ".join(values) + " |") lines.append("") + + try: + budget = budget_comparison(results_dir) + except FileNotFoundError: + return "\n".join(lines) + + lines += [ + "## Sampled ATP and redox budget", + "", + "mmol/gDW/h, GECKO 2 / GECKO 4.", + "", + "| | " + " | ".join(CONDITION_ORDER) + " |", + "|---" * (len(CONDITION_ORDER) + 1) + "|", + ] + for row in budget.itertuples(index=False): + values = [ + f"{getattr(row, f'{c}_gecko2'):.3g} / {getattr(row, f'{c}_gecko4'):.3g}" + for c in CONDITION_ORDER + ] + lines.append(f"| {row.row} | " + " | ".join(values) + " |") + lines.append("") return "\n".join(lines) diff --git a/tests/test_compare.py b/tests/test_compare.py index 95a0a06..e8a69b7 100644 --- a/tests/test_compare.py +++ b/tests/test_compare.py @@ -93,3 +93,24 @@ def test_the_report_names_what_it_compares(results_dir): assert "Mean residual" in text for condition in CONDITION_ORDER: assert condition in text + + +def test_the_sampled_budget_is_compared_when_a_run_is_present(results_dir): + """The ATP budget lines up only if both files are read with the same + row names; a mismatch would silently produce an empty table.""" + import shutil + + from overflow.compare import BUDGET_ROWS, budget_comparison + + sampling = results_dir / "randomSampling" + sampling.mkdir(exist_ok=True) + legacy = pd.read_csv(LEGACY / "randomSampling" / "selectedFluxes.txt", sep="\t") + legacy.to_csv(sampling / "selectedFluxes.tsv", sep="\t", index=False) + + table = budget_comparison(results_dir) + assert set(table["row"]) == set(BUDGET_ROWS) + for condition in CONDITION_ORDER: + assert table[f"{condition}_gecko2"].to_numpy() == pytest.approx( + table[f"{condition}_gecko4"].to_numpy() + ) + shutil.rmtree(sampling) From bb6ad0a25d7173148cf8fe41d3a7b05e180e60e5 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Fri, 25 Sep 2026 13:33:06 +0200 Subject: [PATCH 3/3] docs: say that the analysis runs with rate fitting, and why --- README.md | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/README.md b/README.md index 36f1e64..9887df0 100644 --- a/README.md +++ b/README.md @@ -139,6 +139,14 @@ one, read from `legacy_matlab/results/`. ## Running the whole analysis +The analysis is run with `--fit-rates`, implemented in `overflow.build` +(`relax_to_measured_rates`): it holds CO2, oxygen and the byproducts at their +measurements and raises the measured enzyme abundances by the least that makes +that feasible. Without it those rates are free, and on the earlier tutorial model +the enzyme-constrained solution disposed of surplus carbon through unmeasured +exits instead of respiring it, so the flux distribution did not describe the +measured physiology. `--uptake-flex 1.08` gives CN4 the glucose it needs. + ```bash python -m overflow.build --fit-rates --rate-tolerance 0.08 --uptake-flex 1.08 python -m overflow.build_ribosome