Phase 6: random sampling

README.md

@@ -110,6 +110,23 @@ 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.
 
+## Random sampling
+
+```bash
+python -m overflow.run_sampling --procs 16
+```
+
+Samples the conventional model under the measured rates, with and without
+formate in the measured set, and writes the sampled means and standard
+deviations, the byproducts the model secretes when it is not told to, and the
+ATP and redox budget.
+
+Sampling draws vertices by maximising small random objectives and then minimising
+total flux at each, as the published analysis did; means over such draws are vertex
+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.
+
 ## Tests
 
 ```bash

src/overflow/atp.py

@@ -0,0 +1,140 @@
+"""Where the ATP comes from and what turns over the redox cofactors.
+
+The sampled flux distributions are summarised into a budget: how much
+ATP respiration makes, how much substrate-level phosphorylation makes,
+and how much growth and maintenance spend. Expressed per unit of glucose
+and per unit of growth, so conditions with different uptake rates can be
+compared.
+"""
+from __future__ import annotations
+
+from typing import TYPE_CHECKING, Mapping, Optional, Sequence
+
+import numpy as np
+import pandas as pd
+
+from overflow.config import BIO_RXN, NGAM_RXN
+
+if TYPE_CHECKING:  # pragma: no cover - typing only
+    import cobra
+
+#: ATP made by the respiratory chain and by substrate-level
+#: phosphorylation in the TCA cycle.
+RESPIRATION = ("r_0226", "r_1022")
+
+#: ATP made in glycolysis, and the steps that spend it. yeast-GEM merged
+#: glucokinase (r_4235 in 8.3.4) into hexokinase, so r_0534 now carries
+#: both and is counted once.
+GLYCOLYSIS_PRODUCING = ("r_0892", "r_0962")
+GLYCOLYSIS_CONSUMING = ("r_0886", "r_0534")
+
+GLUCOSE_TRANSPORT = "r_1166"
+
+#: Reactions whose flux the published summary reports directly.
+TURNOVER = {
+    "rPDH": "r_0961",
+    "rIDH": "r_0658",
+    "rMDHc": "r_0714",
+    "rMDHm": "r_0713",
+    "rNDE": "r_0770",
+}
+
+
+def _sum(means: Mapping[str, float], reactions: Sequence[str]) -> float:
+    return float(sum(float(means.get(r, 0.0)) for r in reactions))
+
+
+def atp_budget(
+    means: Mapping[str, float],
+    gam: float,
+    growth_rate: float,
+    glucose_rate: Optional[float] = None,
+    respiration: Sequence[str] = RESPIRATION,
+    glycolysis_producing: Sequence[str] = GLYCOLYSIS_PRODUCING,
+    glycolysis_consuming: Sequence[str] = GLYCOLYSIS_CONSUMING,
+    glucose_transport: str = GLUCOSE_TRANSPORT,
+    turnover: Mapping[str, str] = TURNOVER,
+) -> pd.Series:
+    """The ATP budget implied by one flux distribution.
+
+    ``gam`` is the growth-associated ATP cost written into the biomass
+    reaction. Rates are mmol ATP/gDW/h; yields are per mmol glucose and
+    per unit of growth.
+    """
+    if glucose_rate is None:
+        glucose_rate = float(means.get(glucose_transport, 0.0))
+
+    respiration_rate = _sum(means, respiration)
+    glycolysis = _sum(means, glycolysis_producing) - _sum(means, glycolysis_consuming)
+    growth_cost = gam * float(means.get(BIO_RXN, 0.0))
+    maintenance = float(means.get(NGAM_RXN, 0.0))
+    metabolism = respiration_rate + glycolysis - growth_cost - maintenance
+
+    rows: dict[str, float] = {"rGlu": glucose_rate}
+    for name, rate in (
+        ("ETC", respiration_rate),
+        ("glycolysis", glycolysis),
+        ("GAEC", growth_cost),
+        ("NGAM", maintenance),
+        ("Metabolism", metabolism),
+    ):
+        rows[f"{name}_rATP"] = rate
+        rows[f"{name}_YATP_glu"] = rate / glucose_rate if glucose_rate else np.nan
+        rows[f"{name}_YATP_mu"] = rate / growth_rate if growth_rate else np.nan
+
+    spent = growth_cost + maintenance + metabolism
+    rows["GAEC+NGAM+Metabolism_rATP"] = spent
+    rows["GAEC+NGAM+Metabolism_YATP_glu"] = spent / glucose_rate if glucose_rate else np.nan
+    rows["GAEC+NGAM+Metabolism_YATP_mu"] = spent / growth_rate if growth_rate else np.nan
+
+    for name, reaction in turnover.items():
+        rows[name] = float(means.get(reaction, 0.0))
+
+    return pd.Series(rows)
+
+
+def cofactor_turnover(
+    model: "cobra.Model",
+    means: Mapping[str, float],
+    metabolite_name: str,
+    label: str,
+) -> pd.Series:
+    """Production of one cofactor per compartment.
+
+    Only production counts: consumption is its mirror image, and summing
+    both would report zero everywhere.
+    """
+    metabolites = [m for m in model.metabolites if m.name == metabolite_name]
+    rows: dict[str, float] = {}
+    total = 0.0
+    for metabolite in metabolites:
+        produced = 0.0
+        for reaction in metabolite.reactions:
+            flux = float(means.get(reaction.id, 0.0))
+            rate = reaction.metabolites[metabolite] * flux
+            if rate > 0:
+                produced += rate
+        if produced == 0:
+            continue
+        rows[f"{label}[{metabolite.compartment}]"] = produced
+        total += produced
+    rows[f"{label}[tot]"] = total
+    return pd.Series(rows)
+
+
+def selected_fluxes(
+    model: "cobra.Model",
+    means: Mapping[str, float],
+    gam: float,
+    growth_rate: float,
+    polymerization: Optional[Mapping[str, float]] = None,
+) -> pd.Series:
+    """The published summary of one condition's sampled fluxes."""
+    parts = [
+        atp_budget(means, gam, growth_rate),
+        cofactor_turnover(model, means, "NAD", "rNAD"),
+        cofactor_turnover(model, means, "NADP(+)", "rNADP"),
+    ]
+    if polymerization:
+        parts.append(pd.Series({f"GAECpol_{k}": v for k, v in polymerization.items()}))
+    return pd.concat(parts)

src/overflow/run_sampling.py

@@ -0,0 +1,227 @@
+"""Random sampling of every condition.
+
+Two passes per condition: one with formate excluded from the measured
+set, which asks the model which byproducts it would make if it were not
+told, and one with every measured rate applied.
+
+``python -m overflow.run_sampling``.
+"""
+from __future__ import annotations
+
+import argparse
+import json
+import time
+from pathlib import Path
+from typing import Optional
+
+import pandas as pd
+
+from overflow.adapter import load_gem
+from overflow.atp import selected_fluxes
+from overflow.biomass import GAM_NO_POLYMERIZATION, current_gam
+from overflow.config import (
+    BIO_RXN,
+    CONDITION_ORDER,
+    RESULTS_DIR,
+    load_conditions,
+)
+from overflow.sampling import (
+    polymerization_breakdown,
+    prepare_model,
+    sample_condition,
+    use_fork_start_method,
+)
+
+#: Mean flux above which a byproduct counts as predicted [mmol/gDW/h].
+DETECTION = 0.001
+
+
+def load_good_reactions(cache: Optional[Path], name: str) -> Optional[list[str]]:
+    """The cached loop-free reaction set, if one was kept."""
+    if cache is None:
+        return None
+    path = cache / f"goodReactions_{name}.json"
+    if not path.is_file():
+        return None
+    return json.loads(path.read_text())
+
+
+def save_good_reactions(cache: Optional[Path], name: str, reactions: list[str]) -> None:
+    if cache is None or not reactions:
+        return
+    cache.mkdir(parents=True, exist_ok=True)
+    (cache / f"goodReactions_{name}.json").write_text(json.dumps(sorted(reactions)))
+
+
+def run(
+    conditions: Optional[list[str]] = None,
+    n_samples: int = 5000,
+    seed: Optional[int] = 0,
+    n_proc: Optional[int] = None,
+    cache: Optional[Path] = None,
+    replace_max_bound: bool = False,
+    min_flux: bool = True,
+    loopless: bool = True,
+    verbose: bool = True,
+) -> dict[str, dict]:
+    """Sample each condition twice and summarise.
+
+    Which reactions are free of loops depends on the network and not on
+    the condition, so that screening is done once and reused; ``cache``
+    keeps it across runs, where it is otherwise the dominant cost.
+    """
+    all_conditions = load_conditions()
+    results: dict[str, dict] = {}
+    # With loop-free bounding the objective set is derived per condition
+    # from that condition's own analysis, so a cached one would be ignored.
+    good_free = None if loopless else load_good_reactions(cache, "free")
+    good_full = None if loopless else load_good_reactions(cache, "full")
+    if verbose and (good_free or good_full):
+        print(
+            f"reusing cached loop-free sets: "
+            f"{len(good_free or [])} / {len(good_full or [])} reactions",
+            flush=True,
+        )
+
+    for name in conditions or CONDITION_ORDER:
+        condition = all_conditions[name]
+        started = time.time()
+
+        model = load_gem()
+        prepare_model(model, condition, gam_base=GAM_NO_POLYMERIZATION)
+        gam = current_gam(model)
+        polymerization = polymerization_breakdown(model)
+
+        free, good_free = sample_condition(
+            model.copy(), condition, n_samples=n_samples, include_formate=False,
+            seed=seed, good_reactions=good_free, n_proc=n_proc,
+            replace_max_bound=replace_max_bound, min_flux=min_flux,
+            loopless=loopless,
+        )
+        full, good_full = sample_condition(
+            model, condition, n_samples=n_samples, include_formate=True,
+            seed=seed, good_reactions=good_full, n_proc=n_proc,
+            replace_max_bound=replace_max_bound, min_flux=min_flux,
+            loopless=loopless,
+        )
+
+        summary = selected_fluxes(
+            model, full.means, gam=gam, growth_rate=condition.d_rate,
+            polymerization=polymerization,
+        )
+        if not loopless:
+            save_good_reactions(cache, "free", good_free)
+            save_good_reactions(cache, "full", good_full)
+        results[name] = {
+            "free": free, "full": full, "summary": summary,
+            "gam": gam, "model": model,
+        }
+        if verbose:
+            print(
+                f"[{name}] {full.n_samples} samples, GAM {gam:.2f}, "
+                f"maintenance up to {full.max_maintenance:.2f}; "
+                f"{full.loopless_tightened} reactions held to their loop-free range; "
+                f"growth {full.means.get(BIO_RXN, float('nan')):.5f} "
+                f"[{time.time() - started:.0f}s]",
+                flush=True,
+            )
+    return results
+
+
+def alternative_exchanges(results: dict[str, dict], model, threshold: float = DETECTION) -> pd.DataFrame:
+    """Byproducts the model secretes when it is not told to.
+
+    Taken from the pass with formate left out of the measured set.
+    """
+    exchanges = [r for r in model.reactions if len(r.metabolites) == 1]
+    rows = []
+    for reaction in exchanges:
+        means = {
+            name: float(result["free"].means.get(reaction.id, 0.0))
+            for name, result in results.items()
+        }
+        if max(means.values()) <= threshold:
+            continue
+        rows.append({"rxnID": reaction.id, "rxnName": reaction.name, **means})
+    return pd.DataFrame(rows).sort_values("rxnID")
+
+
+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("--results-dir", type=Path, default=RESULTS_DIR)
+    parser.add_argument("--samples", type=int, default=5000)
+    parser.add_argument("--seed", type=int, default=0)
+    parser.add_argument("--procs", type=int, default=None)
+    parser.add_argument(
+        "--cache", type=Path, default=None,
+        help="directory to keep the loop-free reaction screening in",
+    )
+    parser.add_argument(
+        "--replace-max-bound", action="store_true",
+        help="open the arbitrary 1000 bounds to infinity, as the published "
+             "analysis did; the sampler raises if a random objective then "
+             "turns out unbounded",
+    )
+    parser.add_argument(
+        "--no-min-flux", dest="min_flux", action="store_false",
+        help="skip the minimisation of total flux within each draw that the "
+             "published second pass did",
+    )
+    parser.add_argument(
+        "--keep-loops", action="store_true",
+        help="do not tighten reactions to their loop-free range first",
+    )
+    parser.add_argument("--solver")
+    args = parser.parse_args(argv)
+
+    if args.procs:
+        use_fork_start_method()
+
+    if args.solver or args.procs:
+        import cobra
+
+        if args.solver:
+            cobra.Configuration().solver = args.solver
+        if args.procs:
+            # The loop-free screening is a flux variability analysis, which
+            # reads its process count from here rather than from the
+            # sampler's argument. Left unset it runs on one core and
+            # dominates the run.
+            cobra.Configuration().processes = args.procs
+
+    results = run(
+        args.conditions, n_samples=args.samples, seed=args.seed,
+        n_proc=args.procs, cache=args.cache,
+        replace_max_bound=args.replace_max_bound, min_flux=args.min_flux,
+        loopless=not args.keep_loops,
+    )
+
+    out = args.results_dir / "randomSampling"
+    out.mkdir(parents=True, exist_ok=True)
+
+    any_model = next(iter(results.values()))["model"]
+    alternative_exchanges(results, any_model).to_csv(
+        out / "altExchangeFlux.tsv", sep="\t", index=False
+    )
+
+    wide = pd.DataFrame({"rxnID": next(iter(results.values()))["full"].reactions})
+    wide["rxnName"] = [any_model.reactions.get_by_id(r).name for r in wide["rxnID"]]
+    for name, result in results.items():
+        frame = result["full"].to_frame().set_index("rxnID")
+        wide[f"{name}_AVERAGE"] = wide["rxnID"].map(frame["mean"])
+        wide[f"{name}_STDEV"] = wide["rxnID"].map(frame["sd"])
+        wide[f"{name}_STDEV_percent"] = wide["rxnID"].map(frame["sd_percent"])
+    wide.to_csv(out / "allFluxes.tsv", sep="\t", index=False)
+
+    summary = pd.DataFrame({name: r["summary"] for name, r in results.items()})
+    summary.index.name = "Row"
+    summary.to_csv(out / "selectedFluxes.tsv", sep="\t")
+
+    print()
+    print(summary.head(24).round(4).to_string())
+    return 0
+
+
+if __name__ == "__main__":  # pragma: no cover
+    raise SystemExit(main())