Sample conditions in parallel, and shorten the loop-search limit

README.md

@@ -128,6 +128,18 @@ 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.
 
+Conditions are independent, so they can be sampled as separate runs (for example one
+Slurm job each) and combined afterwards:
+
+```bash
+python -m overflow.run_sampling CN4 --procs 12 --results-dir runs/CN4   # one per condition
+python -m overflow.run_sampling --merge runs/CN4 runs/CN22 runs/CN38 runs/CN75 runs/hGR
+```
+
+The loop-free screening that precedes sampling has a five-minute limit and retries
+with a new seed, because cobra's search for cyclic reactions occasionally runs for
+hours.
+
 ## Comparing against the published results
 
 ```bash

src/overflow/run_sampling.py

@@ -5,6 +5,9 @@
 told, and one with every measured rate applied.
 
 ``python -m overflow.run_sampling``.
+
+Conditions are independent, so they can be sampled as separate runs, each into
+its own results directory, and combined afterwards with ``--merge``.
 """
 from __future__ import annotations
 
@@ -147,6 +150,40 @@ def alternative_exchanges(results: dict[str, dict], model, threshold: float = DE
     return pd.DataFrame(rows).sort_values("rxnID")
 
 
+TABLES = ("allFluxes", "selectedFluxes", "altExchangeFlux")
+
+
+def merge_sampling(sources: list[Path], destination: Path) -> None:
+    """Combine the tables of separate per-condition runs into one set.
+
+    Every source is a results directory holding ``randomSampling/``. In
+    ``altExchangeFlux`` a reaction is listed only where a condition secretes
+    it above the detection threshold, so a condition that does not is left
+    empty for that row.
+    """
+    def read(source: Path, table: str) -> pd.DataFrame:
+        return pd.read_csv(Path(source) / "randomSampling" / f"{table}.tsv", sep="\t")
+
+    fluxes = read(sources[0], "allFluxes")
+    summary = read(sources[0], "selectedFluxes")
+    alternative = read(sources[0], "altExchangeFlux")
+    for source in sources[1:]:
+        part = read(source, "allFluxes")
+        if list(part["rxnID"]) != list(fluxes["rxnID"]):
+            raise ValueError(f"{source} was sampled on a different set of reactions")
+        fluxes = fluxes.merge(part.drop(columns="rxnName"), on="rxnID")
+        summary = summary.merge(read(source, "selectedFluxes"), on="Row")
+        alternative = alternative.merge(
+            read(source, "altExchangeFlux").drop(columns="rxnName"), on="rxnID", how="outer"
+        )
+
+    out = Path(destination) / "randomSampling"
+    out.mkdir(parents=True, exist_ok=True)
+    fluxes.to_csv(out / "allFluxes.tsv", sep="\t", index=False)
+    summary.to_csv(out / "selectedFluxes.tsv", sep="\t", index=False)
+    alternative.sort_values("rxnID").to_csv(out / "altExchangeFlux.tsv", sep="\t", index=False)
+
+
 def main(argv: Optional[list[str]] = None) -> int:
     parser = argparse.ArgumentParser(description=__doc__)
     parser.add_argument("conditions", nargs="*", default=list(CONDITION_ORDER))
@@ -174,8 +211,17 @@ def main(argv: Optional[list[str]] = None) -> int:
         help="do not tighten reactions to their loop-free range first",
     )
     parser.add_argument("--solver", help=HELP)
+    parser.add_argument(
+        "--merge", nargs="+", type=Path, metavar="RESULTS_DIR",
+        help="do not sample: combine the randomSampling tables of these "
+             "per-condition runs into --results-dir",
+    )
     args = parser.parse_args(argv)
 
+    if args.merge:
+        merge_sampling(args.merge, args.results_dir)
+        return 0
+
     if args.procs:
         use_fork_start_method()
 

src/overflow/sampling.py

@@ -210,7 +210,7 @@ def loopless_bounds(
     processes: Optional[int] = None,
     seed: int = 0,
     attempts: int = 4,
-    seconds: float = 900.0,
+    seconds: float = 300.0,
 ) -> pd.DataFrame:
     """Flux range of every reaction that does not need a closed cycle.
 

tests/test_sampling.py

@@ -298,3 +298,54 @@ def test_the_tolerant_step_minimises_total_flux_when_it_can():
     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)
+
+
+def _run_dir(path, condition, flux, alternative):
+    folder = path / condition / "randomSampling"
+    folder.mkdir(parents=True)
+    pd.DataFrame(
+        {"rxnID": ["r1", "r2"], "rxnName": ["a", "b"],
+         f"{condition}_AVERAGE": flux, f"{condition}_STDEV": [0.1, 0.2]}
+    ).to_csv(folder / "allFluxes.tsv", sep="\t", index=False)
+    pd.DataFrame({"Row": ["rGlu", "ETC_rATP"], condition: flux}).to_csv(
+        folder / "selectedFluxes.tsv", sep="\t", index=False
+    )
+    pd.DataFrame({"rxnID": list(alternative), "rxnName": ["x"] * len(alternative),
+                  condition: list(alternative.values())}).to_csv(
+        folder / "altExchangeFlux.tsv", sep="\t", index=False
+    )
+    return path / condition
+
+
+def test_separate_condition_runs_merge_into_one_set(tmp_path):
+    from overflow.run_sampling import merge_sampling
+
+    first = _run_dir(tmp_path, "CN4", [1.0, 2.0], {"r_formate": 0.3})
+    second = _run_dir(tmp_path, "CN22", [3.0, 4.0], {"r_glycine": 0.2})
+    merge_sampling([first, second], tmp_path / "merged")
+    out = tmp_path / "merged" / "randomSampling"
+
+    fluxes = pd.read_csv(out / "allFluxes.tsv", sep="\t")
+    assert list(fluxes.columns) == ["rxnID", "rxnName", "CN4_AVERAGE", "CN4_STDEV",
+                                    "CN22_AVERAGE", "CN22_STDEV"]
+    assert list(fluxes["CN22_AVERAGE"]) == [3.0, 4.0]
+
+    summary = pd.read_csv(out / "selectedFluxes.tsv", sep="\t")
+    assert list(summary.columns) == ["Row", "CN4", "CN22"]
+
+    alternative = pd.read_csv(out / "altExchangeFlux.tsv", sep="\t").set_index("rxnID")
+    assert alternative.loc["r_formate", "CN4"] == pytest.approx(0.3)
+    assert pd.isna(alternative.loc["r_formate", "CN22"])
+    assert alternative.loc["r_glycine", "CN22"] == pytest.approx(0.2)
+
+
+def test_runs_on_different_reactions_are_not_merged(tmp_path):
+    from overflow.run_sampling import merge_sampling
+
+    first = _run_dir(tmp_path, "CN4", [1.0, 2.0], {"r_formate": 0.3})
+    second = _run_dir(tmp_path, "CN22", [3.0, 4.0], {"r_formate": 0.3})
+    table = pd.read_csv(second / "randomSampling" / "allFluxes.tsv", sep="\t")
+    table["rxnID"] = ["r1", "r9"]
+    table.to_csv(second / "randomSampling" / "allFluxes.tsv", sep="\t", index=False)
+    with pytest.raises(ValueError, match="different set of reactions"):
+        merge_sampling([first, second], tmp_path / "merged")