cdanielmachado / kefir_paper

Supplementary Material for Blasche et al
2 stars 1 forks source link

Gapfilling of models #1

Closed franciscozorrilla closed 3 years ago

franciscozorrilla commented 3 years ago

Hi Daniel,

Could you provide further details regarding how the models were reconstructed? In particular, were the models gap-filled on this milk media?

I am specifically simulating the 4 L. lactis models with the script below.

First I simulated each model using the milk media metabolites as the environment. However, I noticed that some exchange reactions are unconstrained for metabolites not present in the milk media file. I manually set the bounds for these reactions to 0 for the second simulation. However, I find that only 1/4 models are able to grow after manually constraining exchange reactions of metabolites not present in the milk media.

import sys
import glob
import cobra
import reframed
from reframed import load_cbmodel
from reframed import Environment
from reframed import FBA

# Loop to read in models in hardcoded folder with .xml extension
for file in glob.iglob(r'*.xml'):

    # Use try to handle exceptions with reading SMBL model
    try:

        # Load model
        model = load_cbmodel(file, flavor='cobra')

        # SIM 1: Use Milk media env
        env = Environment.from_compounds(["h2o","o2","co2","ca2","cl","cobalt2","cu2","fe2","fe3","h","k","mg2","mn2","mobd","na1","nh4","ni2","pi","so4","zn2","ala__L","asn__L","asp__L","glu__L","gln__L","gly","his__L","ile__L","leu__L","lys__L","orn","phe__L","peamn","pro__L","ser__L","thr__L","trp__L","tyr__L","val__L","lcts","glc__D","gal","gal_bD","cit","lac__D","lac__L","for","ac","oxa","pydx","cbl1","thm","pnto__R","fol","ribflv","nac","btn","but","caproic","octa","dca","ddca","ttdca","ptdca","hdca","ocdca","arach","ttdcea","hdcea","ocdcea","lnlc","arachd","ade","gua","ins","thymd","ura","xan","4abut"])

        # SIM 1: Solve FBA
        sol = FBA(model, constraints=env)

        # SIM 1: Open text file with .exf (exchange fluxes) extension
        sys.stdout = open(file+"_milk.exf", "w")

        # SIM 1: Get exchange fluxes
        exf = sol.show_values(pattern="R_EX", sort=True)

        # SIM 2: Constrain uptake of nutrients not in milk media
        model.reactions.R_EX_glyald_e.lb = 0
        model.reactions.R_EX_glyc3p_e.lb = 0
        model.reactions.R_EX_gthrd_e.lb = 0
        model.reactions.R_EX_h2o2_e.lb = 0
        model.reactions.R_EX_acald_e.lb = 0
        model.reactions.R_EX_ala__D_e.lb = 0
        model.reactions.R_EX_mal__L_e.lb = 0
        #model.reactions.R_EX_2pglyc_e.lb = 0

        # SIM 2: Run FBA
        sol = FBA(model)

        # SIM 2: Open text file with .exf (exchange fluxes) extension
        sys.stdout = open(file+"_milkm.exf", "w")

        # SIM 2: Get exchange fluxes
        exf = sol.show_values(pattern="R_EX", sort=True)

    # Catch cobra.io.sbml.CobraSBMLError and continue loop    
    except cobra.io.sbml.CobraSBMLError:

        # Print message with model ID for log
        print("The model",file,"raised an SBML error and could not be read in by COBRA")

        #Nothing to see here fellas
        pass

Best wishes, Francisco

cdanielmachado commented 3 years ago

Hi Francisco,

It is possible that the milk formulation used for model reconstruction and then later for simulations was not exactly the same. I changed the formulation a few times after several discussions with Sonja and Ahmad about what should be present in milk.

Here are two other formulations that you might want to try: milk v1 milk v2.

By the way, the way you are creating the environment, it will only update the exchange reactions associated with the compounds, and leave all other exchanges in the model as they are. If you want to make sure that the model uses exactly what is in the medium you can use one of these two options:

Environment.empty(model, inplace=True)  # reset all exchanges to zero
env = Environment.from_compounds(cpds)
sol = FBA(model, constraints=env)

or

env = Environment.from_compounds(cpds)
env.apply(model, exclusive=True) ## exclusive will block anything not specified in the medium
sol = FBA(model)

Also, I am not really sure why you are reading the model with reframed and using cobrapy to catch exceptions 😅

franciscozorrilla commented 3 years ago

Thanks Daniel! I will use those methods for ensuring that models only uptake what is present in the simulation media.

However, assuming that all models were initialized with a given simulation media within the CarveMe call, then I shouldn't need to worry about exchange reactions for metabolites not in the simulation media in my script simulations, correct?

Ah yes, that must be some vestigial code from when I switched from cobrapy to reframed 🧐

cdanielmachado commented 3 years ago

Yes, in fact you can do the following to check which compounds are in the medium used to initialise the models: Environment.from_model(model).get_compounds()