klamt-lab / straindesign

StrainDesign is a python package for the computational design of metabolic networks and based on COBRApy
Apache License 2.0
36 stars 6 forks source link

OptCouple iMM904 #23

Closed EmanuelGoncalves closed 6 months ago

EmanuelGoncalves commented 6 months ago

Hi,

I'm testing some strain optimization using iMM904 model and managed to run OptKnock but run into some problems with OptCouple (Python 3.10.12, cobrapy 0.26.2, straindesign 1.11).

Am I missing anything?

Thank you!

import cobra
import straindesign as sd

rxnBiomass = "BIOMASS_SC5_notrace"
rxnTarget = "EX_etoh_e"
rxnGlucose = "EX_glc__D_e"
rxnOxygen = "EX_o2_e"

# !wget http://bigg.ucsd.edu/static/models/iMM904.xml
model = cobra.io.read_sbml_model("iMM904.xml")

medium_changed = model.medium.copy()
medium_changed[rxnGlucose] = 10.0
medium_changed[rxnOxygen] = 1.0
model.medium = medium_changed

reactions_to_exclude = set.union(
    set(r.id for r in model.reactions if len(r.compartments) > 1),
    set(r.id for r in model.reactions if len(r.metabolites) == 1),
    set(cobra.flux_analysis.find_blocked_reactions(model)),
)

reactions_to_test = {
    r.id: 1 for r in model.reactions if r.id not in reactions_to_exclude
}

module_optcouple = sd.SDModule(
    model,
    sd.names.OPTCOUPLE,
    prod_id=f"{rxnTarget}",
    inner_objective=f"{rxnBiomass}",
)

Output:

WARNING:root:  Removing reaction bounds when larger than the cobra-threshold of 1000.
ERROR:root:FVA problem not feasible.
WARNING:root:  Removing reaction bounds when larger than the cobra-threshold of 1000.
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[18], line 1
----> 1 sols = sd.compute_strain_designs(
      2     model,
      3     M=1000,
      4     max_cost=1,
      5     ko_cost=reactions_to_test,
      6     sd_modules=module_optcouple,
      7     solution_approach=sd.names.BEST,
      8 )

File /opt/miniconda3/envs/cobrapy/lib/python3.10/site-packages/straindesign/compute_strain_designs.py:391, in compute_strain_designs(model, **kwargs)
    388 logging.info("  Model size: " + str(len(cmp_model.reactions)) + " reactions, " + str(len(cmp_model.metabolites)) + " metabolites")
    389 logging.info("  " + str(len(cmp_ko_cost) + len(cmp_ki_cost) - len(essential_kis)) + " targetable reactions")
--> 391 sd_milp = SDMILP(cmp_model, sd_modules, **kwargs_milp)
    393 kwargs_computation = {}
    394 if MAX_SOLUTIONS in kwargs:

File /opt/miniconda3/envs/cobrapy/lib/python3.10/site-packages/straindesign/strainDesignMILP.py:89, in SDMILP.__init__(self, model, sd_modules, **kwargs)
     87 def __init__(self, model: Model, sd_modules: List[SDModule], **kwargs):
     88     # Construct problem
---> 89     SDProblem.__init__(self, model, sd_modules, **kwargs)
     90     # Build MILP object from constructed problem
     91     MILP_LP.__init__(self,
     92                      c=self.c,
     93                      A_ineq=self.A_ineq,
   (...)
    101                      M=self.M,
    102                      solver=self.solver)

File /opt/miniconda3/envs/cobrapy/lib/python3.10/site-packages/straindesign/strainDesignProblem.py:176, in SDProblem.__init__(self, model, sd_modules, *args, **kwargs)
    174 logging.info('Constructing strain design MILP for solver: ' + self.solver + '.')
    175 for i in range(len(sd_modules)):
--> 176     self.addModule(sd_modules[i])
    178 # Assign knock-ins/outs correctly by taking into account z_inverted
    179 # invert *(-1) rows in z_map_constr_eq, z_map_constr_ineq, z_map_vars
    180 # where there are "knock-ins" / additions
    181 # make knock-in/out matrix
    182 z_kos_kis = [
    183     1 if (not self.z_non_targetable[i]) and (not self.z_inverted[i]) else -1 if self.z_inverted[i] else 0
    184     for i in range(0, self.num_z)
    185 ]

File /opt/miniconda3/envs/cobrapy/lib/python3.10/site-packages/straindesign/strainDesignProblem.py:368, in SDProblem.addModule(self, sd_module)
    366 if MIN_GCP in sd_module:
    367     A_ineq_q = sparse.vstack((A_ineq_q, sparse.lil_matrix(c_p + [-c for c in c_b])), format='csr')
--> 368     b_ineq_q = b_ineq_q + [-sd_module[MIN_GCP]]
    369     z_map_constr_ineq_q = sparse.hstack((z_map_constr_ineq_q, sparse.csc_matrix((self.num_z, 1))))
    370 # 10. reassign lb, ub where possible

TypeError: bad operand type for unary -: 'NoneType'
VonAlphaBisZulu commented 6 months ago

Hi Emanuel,

Thank you for your request. It's great that you find OptKnock solutions to problems in the genome-scale iMM904. Your problem statement is technically correct. The only thing missing is the "min_gcp" parameter. In your example, you are looking for the best solution ("sd.names.BEST"), so you can set that value to 0.0. If you use "sd.names.ANY" to find non-optimal solutions faster, you may want to put a positive value here.

Nevertheless, I now set a default value of min_gcp that is 0.0 because it makes a lot of sense. It's intuitive to start the search without a specific minimum in mind, just like you did.

On a sidenote: GLPK is probably the weakest of the four solvers, and I would usually recommend gurobi or cplex for models with 250+ reactions. More complex problems can easily run multiple days, even with stronger solvers.

I updated it in the last push 2693072293c012835bb172b190f3ea9e7f5b1a9c. Please test the updated version and if your problem persists, please reopen this issue.