Skills › Scientific Agent Skills › 13c-metabolic-flux › mfa.py
mfa.py
84 lines · skills/13c-metabolic-flux/scripts/mfa.py at 5f2bc03
1#!/usr/bin/env python3
2"""Check, simulate, fit, and profile steady-state 13C metabolic flux models."""
3from __future__ import annotations
5import argparse
6import hashlib
7import importlib.metadata
8import json
9from pathlib import Path
10import sys
12sys.dont_write_bytecode = True
15def main():
16 parser = argparse.ArgumentParser(description=__doc__)
17 sub = parser.add_subparsers(dest="command", required=True)
18 for name in ("check", "simulate", "fit"):
19 p = sub.add_parser(name)
20 p.add_argument("--model", type=Path, required=True)
21 p.add_argument("--data", type=Path, required=True, help="Experiment/measurement JSON")
22 p.add_argument("--output", type=Path, help="JSON output; defaults to stdout")
23 if name == "simulate":
24 p.add_argument("--fluxes", type=Path, required=True, help="JSON mapping of reaction IDs to fluxes")
25 if name == "fit":
26 p.add_argument("--starts", type=int, default=8)
27 p.add_argument("--seed", type=int, default=2026)
28 p.add_argument("--maxiter", type=int, default=1000)
29 p.add_argument("--profile", action="append", default=[], help="Reaction to profile; repeat for multiple reactions")
30 p.add_argument("--profile-points", type=int, default=21)
31 p.add_argument("--profile-starts", type=int, default=4)
32 p.add_argument("--confidence", type=float, default=.95)
33 args = parser.parse_args()
34 from _mfa_model import InputError, Network, Observations, read_json, require
35 from _mfa_fit import FluxSpace, diagnostics, fit, profile
36 import numpy as np
37 try:
38 inputs = [args.model, args.data]
39 if args.command == "simulate":
40 inputs.append(args.fluxes)
41 require(args.output is None or all(args.output.resolve() != path.resolve() for path in inputs),
42 "Output must not overwrite an input file")
43 network = Network(read_json(args.model))
44 observations = Observations(read_json(args.data), network)
45 space = FluxSpace(network)
46 versions = {}
47 for package in ("mfapy", "numpy", "scipy"):
48 try:
49 versions[package] = importlib.metadata.version(package)
50 except importlib.metadata.PackageNotFoundError:
51 versions[package] = "not installed"
52 result = dict(command=args.command, model=network.spec["name"], versions=versions,
53 input_sha256={str(path): hashlib.sha256(path.read_bytes()).hexdigest() for path in inputs})
54 if args.command == "check":
55 result.update(status="input_contract_valid", free_flux_dimensions=space.dimension,
56 independent_measurements=observations.size, warnings=observations.warnings,
57 note="Checks declared mapping consistency and feasibility, not biological correctness or isotope steady state")
58 elif args.command == "simulate":
59 vector = network.check_flux(read_json(args.fluxes))
60 result["predictions"] = {exp["id"]: {key: value.tolist() for key, value in network.predict(vector, exp).items()}
61 for exp in observations.spec["experiments"]}
62 else:
63 fitted = fit(observations, args.starts, args.seed, args.maxiter)
64 result.update(diagnostics(observations, fitted))
65 result["settings"] = dict(starts=args.starts, seed=args.seed, maxiter=args.maxiter,
66 profile_points=args.profile_points, profile_starts=args.profile_starts,
67 confidence=args.confidence)
68 result["profiles"] = [profile(observations, fitted, rid, args.profile_points,
69 args.profile_starts, args.seed, args.maxiter, args.confidence)
70 for rid in dict.fromkeys(args.profile)]
71 text = json.dumps(result, indent=2, allow_nan=False) + "\n"
72 if args.output:
73 args.output.write_text(text, encoding="utf-8")
74 else:
75 print(text, end="")
76 return 0
77 except (InputError, ValueError, TypeError, KeyError, OSError, np.linalg.LinAlgError,
78 FloatingPointError, ZeroDivisionError) as exc:
79 print(f"mfa: {exc}", file=sys.stderr)
80 return 2
83if __name__ == "__main__":
84 raise SystemExit(main())