VonAlphaBisZulu
commented
1 year ago Hi Tyler,
The Yeast GEM is one of the largest models out there, with some very complex gene rules. Honestly, strain design will be challenging, but in my opinion possible for ethanol or other small alcohols.
For guidance, check out the genome scale strain design chapter.
Here my advice:
- You will need Gurobi or CPLEX for this task. SCIP and GPLK cannot handle problems of that size and complexity
- YeastGEM has plenty of exchange reactions that can be, but in practice aren't used (such as taurocholate export). You want to close these, if possible, and only keep potentially real fermentation products open, like lactate, acetate etc. If you wish to be extra careful, you can open up the exchanges for pyruvate and other primary metabolites again, once the first computation runs smoothly.
- Verify your strain design goals. Generate plots of your suppressed and protected region. Until you found a first, arbitrary solution, you want to relax the search criteria as much as possible. In your case, you demanded a possible growth rate of 0.1, but the maximal growth rate in the default yeastGEM is only 0.082... so you'll have to change that.
- For your first computation, stick with reaction-based intervention strategies. If you eventually proceed with gene-based strain designs and face problems, reduce the number and complexity of the GPR (gene-protein-reaction) rules of your model. You'll find examples in the manual (same chapter as before). Again, pick the most relaxed settings for your first run. This means high max costs, only a single solution to return and plenty of time, e.g. 24h. You can tighten the constraints and parameters, once your first computation works.
Here is a script that I scribbled down and that for me. It took about 25 minutes to compute on my computer:
import straindesign as sd
import cobra
import matplotlib.pyplot as plt
import logging
logging.basicConfig(level=logging.INFO)
cobra.Configuration().solver = 'cplex'
# 1) Load model and, if available, set the solver to cplex or gurobi. This will increase the chances
# of sucess enormously
yeastGEM = cobra.io.read_sbml_model('C:/Users/phili/Downloads/yeast-GEM.xml')
yeastGEM.solver = 'cplex'
# r_1761: ethanol
# r_1714: glucose
# r_2111: growth
## 2) Block unnecessary exchange reactions
exchange_reacs = [r for r in yeastGEM.reactions if all(s<0 for s in r.metabolites.values())]
# protect main fermentation products
fermentation_products = {
"(S)-lactate exchange",
"(R)-lactate exchange",
"carbon dioxide exchange",
"acetate exchange",
"ethanol exchange",
"formate exchange",
"H+ exchange",
"water exchange",
"methanol exchange",
"oxygen exchange",
"succinate exchange",
"diphosphate exchange",
"growth",
"glycolaldehyde exchange"}
# shut all exchange fluxes
for r in exchange_reacs:
if r.name not in fermentation_products:
r.upper_bound = 0.0
# 3) Verify that the pathways are operational
sol = sd.fba(yeastGEM,obj='r_1761',obj_sense='max')
print(f"Maximum possible ethanol synthesis rate: {sol.fluxes['r_1761']}.")
print(f"Glucose uptake rate: {sol.fluxes['r_1714']}.")
sol = sd.fba(yeastGEM,obj='r_2111',obj_sense='max')
print(f"Maximum possible growth rate: {sol.fluxes['r_2111']}.")
min_etoh_yield = 0.25
min_growth = 0.01
# Plot suppressed and protected regions
module_suppress = sd.SDModule(yeastGEM,sd.names.SUPPRESS,constraints=f'r_1761 + {min_etoh_yield} r_1714 <= 0')
module_protect = sd.SDModule(yeastGEM,sd.names.PROTECT, constraints=f'r_2111 >= {min_growth}')
# Wild-type plot
datapoints, triang, plot1 = sd.plot_flux_space(yeastGEM,
('r_2111',('r_1761','-r_1714')),
show=False);
_, _, plot2 = sd.plot_flux_space(yeastGEM,
(f'r_2111',('r_1761','-r_1714')),
constraints=f'r_2111 >= {min_growth}',
show=False);
plot2.set_facecolor('#70AD47')
plot2.set_edgecolor('#70AD47')
# pGCP design plot
_, _, plot3 = sd.plot_flux_space(yeastGEM,
('r_2111',('r_1761','-r_1714')),
# The sign of the glucose exchange reaction is flipped since
# reaction is defined in the direction of secretion.
constraints=f'r_1761 + {min_etoh_yield} r_1714 <= 0',
show=False);
plot3.set_facecolor('#ED7D31')
plot3.set_edgecolor('#ED7D31')
# adjust axes limits and show plot
plot3.axes.set_xlim(0, 1.05*max([a[0] for a in datapoints]))
plot3.axes.set_ylim(0, 1.05*max([a[1] for a in datapoints]))
plt.show()
# 4) Compute reaction-based strain designs with relaxed settings
# possible knockout of O2
ko_cost = {r.id : 1.0 for r in yeastGEM.reactions if r.genes}
ko_cost.update({'r_1992': 1.0})
# Compute strain designs
sols = sd.compute_strain_designs(yeastGEM,
sd_modules = [module_suppress,module_protect],
max_solutions = 1,
max_cost = 25,
ko_cost = ko_cost,
time_limit = 86400,
solution_approach = sd.names.ANY)
# 5) Print results
print("Knockouts")
print({yeastGEM.reactions.get_by_id(n).name for n in {k for k in sols.get_reaction_sd()[0].keys()}})
# Wild-type plot
datapoints, triang, plot1 = sd.plot_flux_space(yeastGEM,
('r_2111',('r_1761','-r_1714')),
show=False);
_, _, plot2 = sd.plot_flux_space(yeastGEM,
('r_2111',('r_1761','-r_1714')),
constraints=f'r_2111>={min_growth}',
show=False);
plot2.set_facecolor('#70AD47')
plot2.set_edgecolor('#70AD47')
# pGCP design plot
_, _, plot3 = sd.plot_flux_space(yeastGEM,
('r_2111',('r_1761','-r_1714')),
# The sign of the glucose exchange reaction is flipped since
# reaction is defined in the direction of secretion.
constraints=f'r_1761 + {min_etoh_yield} r_1714 <= 0',
show=False);
plot3.set_facecolor('#ED7D31')
plot3.set_edgecolor('#ED7D31')
# plotting designed strain
interventions = [[{s:1.0},'=',0.0] for s,v in sols.reaction_sd[0].items() if v < 1]
_, _, plot4 = sd.plot_flux_space(yeastGEM,
('r_2111',('r_1761','-r_1714')),
# The sign of the glucose exchange reaction is flipped since
# reaction is defined in the direction of secretion.
constraints=interventions,
show=False);
plot4.set_facecolor('#FFC000')
plot4.set_edgecolor('#FFC000')
# adjust axes limits and show plot
plot4.axes.set_xlim(0, 1.05*max([a[0] for a in datapoints]))
plot4.axes.set_ylim(0, 1.05*max([a[1] for a in datapoints]))
plt.show()Keep in mind that there can be particularities in the yeastGEM that make strain design with that model difficult, and that I am not sure were fixed. For instance, growth was impossible under anaerobic conditions. You might get strange results that circumvent knocking out oxygen uptake by introducing a number of surrogate knockouts, e.g. at the ATP synthase or similar. A model must contain all stoichiometric possibilities, in order for StrainDesign to work, and your computation will only be as good as your input model. If the model says growth without O2 is not possible, it won't return O2 as a knockout.
In my case, the script above found the knockout set:
{'aldehyde dehydrogenase (acetaldehyde, NADP)', 'pyruvate dehydrogenase', 'ATP synthase', 'acetaldehyde dehydrogenase', 'aldehyde dehydrogenase (acetylaldehyde, NAD)'}. Not the most exciting one, I've seen in my life, but it makes sense, given the model limitations.
Axis are y: ethanol yield, x: growth rate
Good luck with your computations
tylerautera
Hi Strain Design Team,
I am working with the Chalmers Yeast Model, and I have been unable to create a scenario where there are solutions found despite several attempts on heterologous pathways/targets and native targets. I am wondering if I am not setting up the model correctly, or if there is something else going on.
Here is a link to the model: Chalmer's Yeast GEM
As a simple example, I was attempting to get SUCP solutions for ethanol, here are the modules I used:
As for the other parameters for compute_strain_design:
It would be great to get some guidance on what might be driving the issue here.