ecFBA-epistasis

Genetic interactions of an enzyme-constrained metabolic model as a function of protein availability.

Code, data and figures for:

Competing constraints on protein availability and nutrient uptake reshape yeast genetic interactions

Almaas E, Kumelj T.

Data archive on Zenodo:

Repository structure

• src/ — simulation, analysis and figure scripts
• data/ — the models, the deletion simulations and every derived table, laid out as described in data/README.md
• figures/ — the figures, as written by the scripts in src/
• environment.yml — Python dependencies

Data: what is here and what is on Zenodo

Everything the scripts read is in data/, with one exception. The deletion simulations are stored as Parquet (data/simulations/*/deletions.parquet, one row per solved linear program, 203 million rows for ecYeastGEM-A); the interaction landscape computed from them is a DuckDB store (data/simulations/*/interactions.duckdb, 151 and 128 MB), which exceeds GitHub's file-size limit and is therefore not in this repository. It is in the Zenodo archive, together with a copy of everything else in data/ and a MANIFEST.sha256 of the whole set. Either

• download the archive and unpack it over data/, or
• rebuild the store from the Parquet that is here:

python src/epistasis_pipeline.py \
  --db-path data/simulations/ecYeastGEM-A/deletions.parquet \
  --out-db data/simulations/ecYeastGEM-A/interactions.duckdb \
  --outdir data/landscape/ecYeastGEM-A

This takes about 40 s and 12 GB of memory on a workstation and reproduces the archived store exactly: every table (singles, eps_pairs, pair_profiles, pool_summary, gene_degree, fba_summary, wt) row for row, with numerical differences below 1e-15. The same command with ecYeastGEM in the three paths rebuilds the store of the unmodified model. Only src/fig3_landscape.py and the pair analyses in section 5 below open the store; the other figures run from the CSV tables directly.

Setup

mamba env create -f environment.yml
mamba activate ecFBA

Python 3.10, reframed 1.5.3, Gurobi 10 (licence required), numpy, pandas, scipy, duckdb, pyarrow, matplotlib. All scripts are run from the repository root, for example

python src/fig3_landscape.py

Solver settings are fixed in src/model_scans.py: FeasibilityTol 1e-8, OptimalityTol 1e-9, NumericFocus 3.

The model

ecYeastGEM-A is ecYeastGEM v8.6.2 with three exchange reactions closed (acetaldehyde r_1631, (R,R)-2,3-butanediol r_1549 and diphosphate r_4527); see Methods. model_scans.load() applies the closure at load time and src/make_ecYeastGEM_A.py writes the adjusted model as SBML (data/models/ecYeastGEM-A.xml). Protein availability Ptot (g gDW-1) enters as the lower bound of the pool exchange reaction, -250 Ptot in mg gDW-1 (sigma f = 0.25).

Pipeline

1. Deletion simulations

src/KO_simulation_slurm_v3.py solves the wild type, all single deletions and all pairwise deletions of the 1159 genes at each protein pool value, for the enzyme-constrained model and for the FBA reference. The simulations were run on a compute cluster under Slurm; the three positional arguments are the task index, the number of tasks and a job name, and with 1 1 test the script runs on a single machine. The main sweep of ecYeastGEM-A (421 Ptot values, 0.04 to 0.88 in steps of 0.002) took about 5,800 CPU hours; the raw output is not archived, the tables built from it are (data/simulations/).

Checking the archive without a cluster. Any machine with the environment and a Gurobi licence can rerun a piece of the simulations and compare it with the archived results:

python src/KO_simulation_slurm_v3.py 1 1 test \
  --block-exchanges R_r_1631,R_r_1549,R_r_4527 \
  --pool-min 75 --pool-max 76 --pool-inc 0.5 --max-genes 40
python src/check_simulation.py

The first command simulates 40 genes, singles and all their pairs, at Ptot 0.300, 0.302 and 0.304 and in FBA, about a minute; --block-exchanges closes the three secretion routes in both models, which is what defines ecYeastGEM-A. The second joins its 2,680 solutions to the same genotypes and pools in data/simulations/ecYeastGEM-A/deletions.parquet and prints the largest difference in growth rate; it is 1e-9 h-1, solver tolerance, and the script exits with PASS. Any gene set, pool range or perturbation the simulation script accepts can be checked the same way. Without --block-exchanges the script simulates the unmodified ecYeastGEM, checked with python src/check_simulation.py --archive data/simulations/ecYeastGEM/deletions.parquet. src/wt_sweep_tight.py solves the wild type alone on a fine grid with reduced costs.

2. Simulation store and interaction landscape (workstation)

python src/build_simulation_store.py \
  --ecfba-glob 'results/KO_ecfba_simulation_<jobid>_*.csv' \
  --fba-file results/KO_fba_simulation_<jobid>_1.csv \
  --wt-from-rows --select-all-files --pool-min 10 --pool-max 220 \
  --outdir data/simulations/ecYeastGEM-A
python src/update_wildtype_rows.py      # fine wild-type sweep into the store
python src/epistasis_pipeline.py \
  --db-path data/simulations/ecYeastGEM-A/deletions.parquet \
  --out-db data/simulations/ecYeastGEM-A/interactions.duckdb \
  --outdir data/landscape/ecYeastGEM-A

The store step turns the raw output files into deletions.parquet and wildtype.csv; the landscape step computes W_i, W_ij and eps for every pair at every pool and writes interactions.duckdb (tables singles, eps_pairs, pair_profiles, pool_summary, gene_degree, fba_summary, wt) together with the CSV summaries in data/landscape/. The landscape step accepts deletions.parquet directly, which is how the store is rebuilt from this repository (see Data above).

3. Wild-type physiology and the phase plane

script	output	used in
model_scans.py	Ptot sweep, glucose-bound scan, one-constraint sweep	Fig 2, Figs A and C in S1 File; Table 1
plane_boundaries.py, plane_regions.py, extend_plane_chemostat.py	regime boundaries of the (nu_glc, Ptot) plane by bisection	Fig 2
fit_ptot_vanhoek.py	least-squares fit of Ptot to the chemostat data of van Hoek et al. (1998)	Fig B and Table A in S1 File
allocation_sweep.py, allocation_robustness.py	pool allocation by subsystem along Ptot, and under kcat perturbation	Fig C in S1 File
recruitment_order.py	order and protein cost of the respiratory modules	Results
check_scaling_invariance.py, one_constraint_doubles.py	homogeneity checks: joint scaling of Ptot and the bounds; pairs with the pool as the only binding constraint	Methods, Table 1
vemuri_aox_nox.py	heterologous oxidases and D_crit	Discussion

4. Single deletions

script	output	used in
ptdep_genes.py, gene_groups.py	the Ptot-dependent genes and their subsystem groups	Fig 3A
mutant_shadow_prices.py, shadow_price_check.py, shadow_price_classes.py	marginal value of protein in every deletion background, with the finite-difference and solver-method checks	Results, Methods
threshold_kcat.py	saturation thresholds of the respiratory deletions under kcat perturbation	Results
dpsat.py	saturation Ptot of every deletion strain	Fig D in S1 File
compare_models_singles.py	single-deletion fitness of ecYeastGEM-A against ecYeastGEM	Discussion

5. Pairs

script	output	used in
landscape_numbers.py	every landscape number quoted in Results	Results, Table 1, Table B in S1 File
sign_changes.py	the sign-changing pairs and their four groups	Fig 3, Fig E in S1 File
pair_classes.py	classification of every called pair at every Ptot	Results
analyze_pair_structure.py	subsystem assignment, pair classes, enrichment, strict pairs	Fig 4, Fig F in S1 File
alleviation_branches.py	alleviation of biosynthetic deletions by respiratory deletions	Results
pair_shadow_prices.py, epistasis_protein_price.py	marginal value of protein in the double mutant and the protein price of an interaction	Results
mechanism_table.py, flux_maps.py, mechanism_cases.py	flux-level mechanism of ADH1 x ADH5 and PRO2 x QCR2	Fig 5, Table C in S1 File
robustness_kcat.py, robustness_pairs.py, robustness_gpm1.py	turnover-number robustness of the strict pairs and of the GPM1 sign inversion	Discussion

analyze_pair_structure.py reads the ecYeastGEM store by default; for ecYeastGEM-A it is run as

PAIRSTRUCT_OUT=data/pairs/ecYeastGEM-A \
PAIRSTRUCT_DB=data/simulations/ecYeastGEM-A/interactions.duckdb \
PAIRSTRUCT_BD=data/analysis/ptdep_genes.csv \
  python src/analyze_pair_structure.py

Reproducing the figures

for f in src/fig*.py; do python "$f"; done

writes Fig1 to Fig5 and S1_Fig to S6_Fig as PDF and PNG to figures/. src/figstyle.py fixes the page width, fonts and the Okabe-Ito palette shared by all figures. The figure scripts read the tables in data/ and finish in a few minutes; only fig3_landscape.py opens the interaction store.

Please cite

If you use this repository, its data, or derived results in academic work, please cite:

Almaas E, Kumelj T. Competing constraints on protein availability and nutrient uptake reshape yeast genetic interactions.

Contact

Eivind Almaas eivind.almaas@ntnu.no Norwegian University of Science and Technology (NTNU), Trondheim, Norway

Licence

MIT; see LICENSE. ecYeastGEM and Yeast-GEM are distributed by their authors under their own licences and are not part of this repository.