diff --git a/README.md b/README.md index 49a72d8..0fa5161 100644 --- a/README.md +++ b/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 diff --git a/src/overflow/atp.py b/src/overflow/atp.py new file mode 100644 index 0000000..49b6487 --- /dev/null +++ b/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) diff --git a/src/overflow/run_sampling.py b/src/overflow/run_sampling.py new file mode 100644 index 0000000..ad5fa62 --- /dev/null +++ b/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()) diff --git a/src/overflow/sampling.py b/src/overflow/sampling.py new file mode 100644 index 0000000..8ff280a --- /dev/null +++ b/src/overflow/sampling.py @@ -0,0 +1,385 @@ +"""Random sampling of the conventional model under the measured rates. + +No enzyme constraints here: the question is what the measured exchange +rates alone imply about the flux distribution, so the model is the plain +yeast-GEM with its biomass rescaled to the condition and every measured +rate held inside a narrow band around its measurement. + +Sampling draws vertices by maximising small random objectives, which is +what the published analysis did. It is not a uniform draw from the +interior; means over such draws are vertex means. +""" +from __future__ import annotations + +import signal +from contextlib import contextmanager +from dataclasses import dataclass, field +from typing import TYPE_CHECKING, Mapping, Optional + +import numpy as np +import pandas as pd +import raven_toolbox.analysis.sampling as raven_sampling +from cobra.exceptions import Infeasible +from cobra.flux_analysis import pfba +from raven_toolbox.analysis.sampling import random_sampling + +from overflow.biomass import ( + GAM_NO_POLYMERIZATION, + POLYMERIZATION_COST, + YEAST_BIOMASS, + apply_protein_content, +) +from overflow.config import BIO_RXN, GROWTH_RXN, NGAM_RXN, Condition + +if TYPE_CHECKING: # pragma: no cover - typing only + import cobra + +#: Exchange reactions held at their measurement, and the sign convention +#: each is measured with: -1 for an uptake, +1 for a secretion. +MEASURED_EXCHANGES = { + "r_1714": ("Glucose", -1), + "r_1992": ("Oxygen", -1), + "r_1672": ("CO2", +1), + "r_1808": ("Glycerol", +1), + "r_1634": ("Acetate", +1), + "r_1761": ("Ethanol", +1), + "r_1793": ("Formate", +1), + GROWTH_RXN: ("D", +1), +} + +#: Formate dehydrogenase, blocked so that formate leaves the cell rather +#: than being oxidised to CO2 inside it. +FORMATE_DEHYDROGENASE = "r_0445" + +FORMATE_EXCHANGE = "r_1793" + + +@contextmanager +def time_limit(seconds: float): + """Raise ``TimeoutError`` in the main thread once ``seconds`` have passed.""" + + def expire(signum, frame): + raise TimeoutError + + previous = signal.signal(signal.SIGALRM, expire) + signal.setitimer(signal.ITIMER_REAL, seconds) + try: + yield + finally: + signal.setitimer(signal.ITIMER_REAL, 0) + signal.signal(signal.SIGALRM, previous) + + +def use_fork_start_method() -> None: + """Have multiprocessing fork rather than use a forkserver. + + Python 3.14 made forkserver the default on Linux. Its helper process + does not survive in every environment -- on this cluster it dies and + takes the run with it, as a broken pipe or a reset connection partway + through sampling. Fork is what cobrapy's process pool was written + against and works here. + """ + import multiprocessing + + if "fork" not in multiprocessing.get_all_start_methods(): + return + if multiprocessing.get_start_method(allow_none=True) == "fork": + return + try: + multiprocessing.set_start_method("fork", force=True) + except RuntimeError: # already started a pool; nothing to do + pass + + +def tolerant_pfba(model: "cobra.Model", fraction_of_optimum: float = 1 - 1e-6): + """Minimise total flux at the optimum, keeping the plain solution if that is infeasible. + + Pinning the objective exactly at its optimum can be numerically + infeasible for a single draw; a small slack avoids most of those and + the plain solution covers the rest. + """ + try: + return pfba(model, fraction_of_optimum=fraction_of_optimum) + except Infeasible: + return model.optimize() + + +def measured_rate(condition: Condition, field_name: str) -> float: + """One measured rate, as a positive magnitude.""" + if field_name == "D": + return condition.d_rate + if field_name == "Glucose": + return condition.glucose + if field_name == "Oxygen": + return condition.oxygen + if field_name == "CO2": + return condition.co2 + return condition.byproducts.get(field_name, 0.0) + + +def sampling_bounds( + condition: Condition, + tolerance: float = 0.05, + minimum_secretion: float = 0.01, + include_formate: bool = True, +) -> dict[str, tuple[float, float]]: + """Bounds holding every measured rate within ``tolerance`` of itself. + + A byproduct that was not detected gets a small allowance rather than + being blocked, so that the sampler can say whether the model wants to + make it. + """ + bounds: dict[str, tuple[float, float]] = {} + for reaction_id, (field_name, sign) in MEASURED_EXCHANGES.items(): + if reaction_id == FORMATE_EXCHANGE and not include_formate: + continue + value = measured_rate(condition, field_name) + low = (1 - tolerance) * value + high = (1 + tolerance) * value + if sign < 0: + low, high = -high, -low + else: + high = high if high != 0 else minimum_secretion + bounds[reaction_id] = (low, high) + return bounds + + +def prepare_model( + model: "cobra.Model", + condition: Condition, + gam_base: float = GAM_NO_POLYMERIZATION, +) -> dict[str, float]: + """Set the condition's biomass and open maintenance. + + Biomass is rescaled to the measured protein content and maintenance + rebuilt from it, which is what this analysis fixes the ATP budget + against. Unlike the enzyme-constrained models there is no protein + budget here to make that unaffordable. + """ + composition = apply_protein_content( + model, condition.p_tot, gam="polymerization", gam_base=gam_base, + scale_protein=True, + ) + model.reactions.get_by_id(FORMATE_DEHYDROGENASE).bounds = (0.0, 0.0) + model.reactions.get_by_id(NGAM_RXN).bounds = (0.0, 1000.0) + return composition + + +def apply_bounds(model: "cobra.Model", bounds: Mapping[str, tuple[float, float]]) -> None: + for reaction_id, (low, high) in bounds.items(): + model.reactions.get_by_id(reaction_id).bounds = (low, high) + + +def constrain_maintenance(model: "cobra.Model", fraction: float = 0.95) -> float: + """Hold maintenance near the most the condition can support. + + Everything else is fixed by the measured rates, so whatever ATP the + model does not spend on growth it spends here; leaving maintenance + free would let it dissipate that instead through whichever internal + cycle happens to be cheapest. + """ + maintenance = model.reactions.get_by_id(NGAM_RXN) + maintenance.bounds = (0.0, 1000.0) + model.objective = NGAM_RXN + model.objective.direction = "max" + solution = model.optimize() + if solution.status != "optimal": + raise RuntimeError("no feasible solution under the measured rates") + highest = float(solution.fluxes[NGAM_RXN]) + maintenance.lower_bound = fraction * highest + return highest + + +def polymerization_breakdown(model: "cobra.Model", gam_base: float = GAM_NO_POLYMERIZATION) -> dict[str, float]: + """ATP spent polymerising each component of one gram of biomass.""" + from raven_toolbox.biomass.scale import sum_biomass + + fractions = sum_biomass(model, YEAST_BIOMASS) + costs = { + component: cost * fractions.get(component, 0.0) + for component, cost in POLYMERIZATION_COST.items() + } + costs["total"] = sum(costs.values()) + costs["GAEC_noPol"] = gam_base + costs["GAEC_total"] = costs["total"] + gam_base + return costs + + +def loopless_bounds( + model: "cobra.Model", + processes: Optional[int] = None, + seed: int = 0, + attempts: int = 4, + seconds: float = 900.0, +) -> pd.DataFrame: + """Flux range of every reaction that does not need a closed cycle. + + Ordinary variability analysis lets a reaction in a thermodynamically + infeasible cycle run to whatever arbitrary bound the model carries; + the loop-free range is what the reaction can do while the rest of + the network stays consistent. + """ + from cobra.flux_analysis import flux_variability_analysis + + if processes is not None: + import cobra + + cobra.Configuration().processes = processes + for attempt in range(attempts): + # cobra picks the cyclic reactions with random LP weights, and how + # long that takes depends on them: most seeds finish in a minute or + # two, an occasional one runs for hours. A fresh seed after a time + # limit gets past it. + np.random.seed(seed + attempt) + try: + with time_limit(seconds): + return flux_variability_analysis( + model, fraction_of_optimum=0.0, loopless="cycleFreeFlux" + ) + except TimeoutError: + continue + raise TimeoutError( + f"the loop-free variability analysis did not finish in {seconds:.0f} s " + f"in any of {attempts} attempts" + ) + + +def good_reactions_from(ranges: pd.DataFrame, tolerance: float = 1e-9) -> list[str]: + """Reactions that can carry flux without a closed cycle. + + These are the ones worth using as random objectives; a reaction that + only moves inside a cycle gives a degenerate solve every time it is + picked, and the sampler retries up to a hundred times before giving + up on that draw. + """ + return [ + reaction_id + for reaction_id, row in ranges.iterrows() + if max(abs(float(row["minimum"])), abs(float(row["maximum"]))) > tolerance + ] + + +def apply_loopless_bounds( + model: "cobra.Model", ranges: pd.DataFrame, tolerance: float = 1e-9 +) -> int: + """Tighten every reaction to its loop-free range. + + Every flux distribution free of closed cycles already lies inside + these bounds, so nothing thermodynamically sensible is excluded -- + what goes is the room a cycle needs to inflate a flux. Returns how + many reactions were tightened. + """ + tightened = 0 + for reaction_id, row in ranges.iterrows(): + reaction = model.reactions.get_by_id(reaction_id) + low = max(reaction.lower_bound, float(row["minimum"]) - tolerance) + high = min(reaction.upper_bound, float(row["maximum"]) + tolerance) + if low > high: # numerical crossing; leave the reaction alone + continue + if low > reaction.lower_bound or high < reaction.upper_bound: + reaction.bounds = (low, high) + tightened += 1 + return tightened + + +@dataclass +class SamplingResult: + """Sampled fluxes for one condition and one formate setting.""" + + condition: str + include_formate: bool + reactions: list[str] = field(default_factory=list) + mean: np.ndarray = field(default_factory=lambda: np.empty(0)) + sd: np.ndarray = field(default_factory=lambda: np.empty(0)) + n_samples: int = 0 + max_maintenance: float = 0.0 + loopless_tightened: int = 0 + + def to_frame(self) -> pd.DataFrame: + with np.errstate(divide="ignore", invalid="ignore"): + relative = np.where(self.mean != 0, self.sd / self.mean * 100, np.nan) + return pd.DataFrame( + {"rxnID": self.reactions, "mean": self.mean, "sd": self.sd, + "sd_percent": relative} + ) + + @property + def means(self) -> pd.Series: + return pd.Series(self.mean, index=self.reactions) + + +def sample_condition( + model: "cobra.Model", + condition: Condition, + n_samples: int = 5000, + include_formate: bool = True, + tolerance: float = 0.05, + seed: Optional[int] = None, + good_reactions: Optional[list[str]] = None, + n_proc: Optional[int] = None, + min_flux: bool = True, + replace_max_bound: bool = False, + loopless: bool = True, +) -> tuple[SamplingResult, list[str]]: + """Sample one condition, returning the summary and the good reactions. + + ``good_reactions`` is the loop-free set the sampler may use as random + objectives; it depends only on the network, so it is computed once + and handed back for reuse. + + ``replace_max_bound`` opens the arbitrary 1000 bounds to infinity, as + the published analysis did, so that a reaction in a loop cannot sit + at 1000 and be mistaken for flux. It is off here: the screening that + picks the random objectives runs on the finite bounds, so once they + are opened an objective can turn out unbounded, and this sampler + raises on that where RAVEN's returned no solution and moved on. + + ``min_flux`` minimises total flux within each draw, as the published + second pass did. Without it the sampled means drift into + high-flux routes: the pentose phosphate pathway and a reversed + NADP isocitrate dehydrogenase carry several times the published + flux. The step is tolerant of a draw whose optimum cannot be pinned + exactly. + """ + apply_bounds(model, sampling_bounds(condition, tolerance, include_formate=include_formate)) + raven_sampling.pfba = tolerant_pfba + highest = constrain_maintenance(model) + + # Without this the malate dehydrogenases cycle against each other at + # several hundred mmol/gDW/h, which is the arbitrary 1000 bound + # showing through rather than anything the cell does. + tightened = 0 + if loopless: + # One loop-free variability analysis gives both the bounds and + # the objective set, and the objective set has to come from this + # condition: a reaction that moves freely in one condition can be + # pinned in another, and picking it then wastes a whole draw. + ranges = loopless_bounds(model, n_proc, seed=seed or 0) + tightened = apply_loopless_bounds(model, ranges) + good_reactions = good_reactions_from(ranges) + + sampled = random_sampling( + model, + n_samples, + method="random_objective", + n_objectives=2, + replace_max_bound=replace_max_bound, + suppress_errors=True, + good_reactions=good_reactions, + min_flux=min_flux, + seed=seed, + n_proc=n_proc, + ) + samples = sampled.samples + result = SamplingResult( + condition=condition.name, + include_formate=include_formate, + reactions=list(samples.columns), + mean=samples.mean(axis=0).to_numpy(), + sd=samples.std(axis=0, ddof=1).to_numpy(), + n_samples=len(samples), + max_maintenance=highest, + loopless_tightened=tightened, + ) + returned = sampled.good_reactions + return result, list(returned) if returned is not None else (good_reactions or []) diff --git a/tests/test_sampling.py b/tests/test_sampling.py new file mode 100644 index 0000000..34a5d8e --- /dev/null +++ b/tests/test_sampling.py @@ -0,0 +1,300 @@ +"""Sampling the conventional model under the measured rates. + +The ATP accounting is arithmetic on a flux vector, so it is checked +against the published summary: fed the committed sampling means it must +reproduce the committed `selectedFluxes.txt`. +""" +import numpy as np +import pandas as pd +import pytest + +from overflow.atp import atp_budget, cofactor_turnover +from overflow.config import CONDITION_ORDER, ROOT, load_conditions +from overflow.sampling import ( + FORMATE_DEHYDROGENASE, + FORMATE_EXCHANGE, + MEASURED_EXCHANGES, + measured_rate, + sampling_bounds, +) +from routes import two_route_model + +LEGACY = ROOT / "legacy_matlab" / "results" / "randomSampling" + +#: yeast-GEM 8.3.4 split glucose phosphorylation over two reactions. +LEGACY_GLYCOLYSIS_CONSUMING = ("r_4235", "r_0886", "r_0534") + + +@pytest.fixture(scope="module") +def legacy_means(): + table = pd.read_csv(LEGACY / "allFluxes.txt", sep="\t") + return { + condition: dict(zip(table["ID"], table[f"{condition}_AVERAGE"])) + for condition in CONDITION_ORDER + } + + +@pytest.fixture(scope="module") +def legacy_selected(): + return pd.read_csv(LEGACY / "selectedFluxes.txt", sep="\t").set_index("Row") + + +# --- the bounds ------------------------------------------------------- + +def test_uptakes_are_negative_and_secretions_positive(conditions): + bounds = sampling_bounds(conditions["CN75"]) + low, high = bounds["r_1714"] + assert low < high < 0, "glucose is taken up" + low, high = bounds["r_1761"] + assert 0 < low < high, "ethanol is secreted" + + +def test_every_measured_rate_is_held_within_five_percent(conditions): + condition = conditions["CN75"] + for reaction_id, (field_name, sign) in MEASURED_EXCHANGES.items(): + value = measured_rate(condition, field_name) + if value == 0: + continue + low, high = bounds = sampling_bounds(condition)[reaction_id] + assert abs(low) == pytest.approx(0.95 * value if sign > 0 else 1.05 * value) + assert abs(high) == pytest.approx(1.05 * value if sign > 0 else 0.95 * value) + + +def test_an_undetected_byproduct_gets_a_small_allowance(conditions): + """Blocking it would prevent the sampler from saying whether the + model wants to make it, which is the question the run asks.""" + low, high = sampling_bounds(conditions["CN4"])["r_1761"] + assert low == 0.0 + assert high == pytest.approx(0.01) + + +def test_formate_can_be_left_out_of_the_measurement_set(conditions): + with_formate = sampling_bounds(conditions["CN22"], include_formate=True) + without = sampling_bounds(conditions["CN22"], include_formate=False) + assert FORMATE_EXCHANGE in with_formate + assert FORMATE_EXCHANGE not in without + assert set(with_formate) - {FORMATE_EXCHANGE} == set(without) + + +def test_growth_is_held_at_the_dilution_rate(conditions): + low, high = sampling_bounds(conditions["hGR"])["r_2111"] + assert low == pytest.approx(0.95 * 0.29) + assert high == pytest.approx(1.05 * 0.29) + + +# --- the ATP accounting, against the published summary --------------- + +@pytest.mark.parametrize("condition", CONDITION_ORDER) +def test_atp_budget_reproduces_the_published_summary( + legacy_means, legacy_selected, condition +): + means = legacy_means[condition] + published = legacy_selected[condition] + gam = float(published["GAEC_rATP"]) / means["r_4041"] + + ours = atp_budget( + means, + gam=gam, + growth_rate=load_conditions()[condition].d_rate, + glycolysis_consuming=LEGACY_GLYCOLYSIS_CONSUMING, + ) + + # The published table is written to four significant figures, so a + # quantity derived from several of them carries their rounding. + for row in ( + "rGlu", "ETC_rATP", "ETC_YATP_glu", "ETC_YATP_mu", + "glycolysis_rATP", "glycolysis_YATP_glu", "glycolysis_YATP_mu", + "GAEC_rATP", "GAEC_YATP_glu", "NGAM_rATP", "NGAM_YATP_glu", + "Metabolism_rATP", "GAEC+NGAM+Metabolism_rATP", + "rPDH", "rIDH", "rMDHc", "rMDHm", "rNDE", + ): + assert float(ours[row]) == pytest.approx(float(published[row]), rel=2e-3), row + + +def test_the_budget_balances(legacy_means, legacy_selected): + """What respiration and glycolysis make is what growth, maintenance + and the rest of metabolism spend.""" + means = legacy_means["CN4"] + gam = float(legacy_selected.loc["GAEC_rATP", "CN4"]) / means["r_4041"] + budget = atp_budget(means, gam=gam, growth_rate=0.1, + glycolysis_consuming=LEGACY_GLYCOLYSIS_CONSUMING) + made = budget["ETC_rATP"] + budget["glycolysis_rATP"] + spent = budget["GAEC+NGAM+Metabolism_rATP"] + assert made == pytest.approx(spent) + + +def test_yields_are_rates_divided_by_uptake_and_growth(legacy_means, legacy_selected): + means = legacy_means["CN38"] + gam = float(legacy_selected.loc["GAEC_rATP", "CN38"]) / means["r_4041"] + budget = atp_budget(means, gam=gam, growth_rate=0.1, + glycolysis_consuming=LEGACY_GLYCOLYSIS_CONSUMING) + assert budget["ETC_YATP_glu"] == pytest.approx(budget["ETC_rATP"] / budget["rGlu"]) + assert budget["ETC_YATP_mu"] == pytest.approx(budget["ETC_rATP"] / 0.1) + + +def test_a_condition_with_no_glucose_does_not_divide_by_zero(): + budget = atp_budget({"r_1166": 0.0, "r_4041": 0.1}, gam=50.0, growth_rate=0.1) + assert np.isnan(budget["ETC_YATP_glu"]) + assert np.isfinite(budget["ETC_YATP_mu"]) + + +# --- cofactor turnover ------------------------------------------------ + +@pytest.mark.slow +def test_cofactor_turnover_counts_production_only(conv_model, legacy_means): + turnover = cofactor_turnover(conv_model, legacy_means["CN4"], "NAD", "rNAD") + assert (turnover >= 0).all() + assert turnover["rNAD[tot]"] == pytest.approx( + turnover.drop("rNAD[tot]").sum() + ) + + +@pytest.mark.slow +def test_cofactor_turnover_is_reported_per_compartment(conv_model, legacy_means): + turnover = cofactor_turnover(conv_model, legacy_means["hGR"], "NAD", "rNAD") + compartments = {k for k in turnover.index if k != "rNAD[tot]"} + assert len(compartments) >= 2, "NAD turns over in more than one compartment" + assert all(k.startswith("rNAD[") for k in compartments) + + +# --- the loop-free screening cache ----------------------------------- + +def test_the_loop_free_screening_round_trips(tmp_path): + """Screening every reaction for loop involvement is the dominant + cost of a sampling run and depends only on the network, so it is + kept rather than repeated.""" + from overflow.run_sampling import load_good_reactions, save_good_reactions + + assert load_good_reactions(tmp_path, "free") is None + save_good_reactions(tmp_path, "free", ["r_2", "r_1"]) + assert load_good_reactions(tmp_path, "free") == ["r_1", "r_2"] + assert load_good_reactions(tmp_path, "full") is None + + +def test_no_cache_directory_means_no_caching(tmp_path): + from overflow.run_sampling import load_good_reactions, save_good_reactions + + save_good_reactions(None, "free", ["r_1"]) + assert load_good_reactions(None, "free") is None + save_good_reactions(tmp_path, "free", []) + assert load_good_reactions(tmp_path, "free") is None + + +# --- holding reactions to their loop-free range ---------------------- + +def test_loopless_bounds_tighten_only_what_a_cycle_inflates(): + """A reaction whose ordinary range reaches the model's arbitrary + bound, but whose loop-free range is a few units, is held to the + latter; a reaction already inside its loop-free range is left be.""" + import cobra + import pandas as pd + + from overflow.sampling import apply_loopless_bounds + + model = cobra.Model("m") + a = cobra.Reaction("a", lower_bound=-1000.0, upper_bound=1000.0) + b = cobra.Reaction("b", lower_bound=0.0, upper_bound=3.0) + model.add_reactions([a, b]) + + ranges = pd.DataFrame( + {"minimum": [-2.079, 0.0], "maximum": [2.908, 5.0]}, index=["a", "b"] + ) + tightened = apply_loopless_bounds(model, ranges) + + assert tightened == 1 + assert a.lower_bound == pytest.approx(-2.079, abs=1e-6) + assert a.upper_bound == pytest.approx(2.908, abs=1e-6) + assert b.bounds == (0.0, 3.0), "already inside its loop-free range" + + +def test_loopless_bounds_never_widen_a_reaction(): + """The loop-free range of a reaction can exceed a bound the + condition imposes; the condition wins.""" + import cobra + import pandas as pd + + from overflow.sampling import apply_loopless_bounds + + model = cobra.Model("m") + r = cobra.Reaction("r", lower_bound=-1.0, upper_bound=1.0) + model.add_reactions([r]) + apply_loopless_bounds( + model, pd.DataFrame({"minimum": [-50.0], "maximum": [50.0]}, index=["r"]) + ) + assert r.bounds == (-1.0, 1.0) + + +def test_a_crossed_range_leaves_the_reaction_alone(): + """Numerical noise can put a loop-free minimum above its maximum; + that must not produce an infeasible reaction.""" + import cobra + import pandas as pd + + from overflow.sampling import apply_loopless_bounds + + model = cobra.Model("m") + r = cobra.Reaction("r", lower_bound=0.0, upper_bound=1.0) + model.add_reactions([r]) + apply_loopless_bounds( + model, pd.DataFrame({"minimum": [5.0], "maximum": [2.0]}, index=["r"]) + ) + assert r.bounds == (0.0, 1.0) + + +def test_a_variability_analysis_that_overruns_is_retried_with_a_new_seed(monkeypatch): + import time + + import cobra.flux_analysis as flux_analysis + import numpy as np + + from overflow.sampling import loopless_bounds + + seeds = [] + + def analysis(model, **kwargs): + seeds.append(int(np.random.get_state()[1][0])) + if len(seeds) == 1: + time.sleep(5) + return "ranges" + + monkeypatch.setattr(flux_analysis, "flux_variability_analysis", analysis) + assert loopless_bounds(object(), seed=7, attempts=3, seconds=0.2) == "ranges" + assert len(seeds) == 2 and seeds[0] != seeds[1] + + +def test_a_variability_analysis_that_never_finishes_is_reported(monkeypatch): + import time + + import cobra.flux_analysis as flux_analysis + + from overflow.sampling import loopless_bounds + + monkeypatch.setattr( + flux_analysis, "flux_variability_analysis", lambda model, **kw: time.sleep(5) + ) + with pytest.raises(TimeoutError, match="2 attempts"): + loopless_bounds(object(), attempts=2, seconds=0.1) + + +def test_a_draw_that_cannot_be_pinned_keeps_the_plain_solution(monkeypatch): + from cobra.exceptions import Infeasible + + import overflow.sampling as sampling + + model = two_route_model() + + def refuse(*args, **kwargs): + raise Infeasible("numerically infeasible") + + monkeypatch.setattr(sampling, "pfba", refuse) + solution = sampling.tolerant_pfba(model) + assert solution.status == "optimal" + assert solution.fluxes["demand"] == pytest.approx(1.0) + + +def test_the_tolerant_step_minimises_total_flux_when_it_can(): + import overflow.sampling as sampling + + solution = sampling.tolerant_pfba(two_route_model()) + assert solution.fluxes["direct"] == pytest.approx(0.0, abs=1e-5) + assert solution.fluxes["step1"] == pytest.approx(1.0, abs=1e-5)