Phase 7: comparison against the published results

README.md

@@ -127,9 +127,44 @@ 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, 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
+
+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
+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.

src/overflow/compare.py

@@ -0,0 +1,226 @@
+"""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)
+
+
+#: 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)
+    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("")
+
+    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)
+
+
+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())

tests/test_compare.py

@@ -0,0 +1,116 @@
+"""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
+
+
+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)