A few things:

1. Should linear MOMA be implemented?
2. Should a modifiable LP be generated for things like succesive deletions?
3. Code needs more extensive unit testing

Right now we use essentially the native pickling function. We should discuss the merits and drawbacks of a more finely tuned serialization format.

Major pros:

• more robust to changes in cobrapy code so we don't accidentally break pickles, which is a potential source for bugs.
• More space efficient. Some ME model pickles are huge.

Cons:

• Need to make sure that new features are handled well. Not doing so would be a source of bugs.
• Need to be careful to not make incorrect assumptions about how cobrapy is used when making pickles.

Unknown:

• performance (could go either way depending on implementation)

I tried to get media from ecoli_pickle with cobrapy(GitHub HEAD), but it seems
it does not contain "media_compositions".
Is there any way to set "media_compositions" to ecoli model?

In [2]: sal_model = cobra.test.create_test_model(cobra.test.salmonella_pickle) 

In [3]: sal_model.media_compositions 
Out[27]: defaultdict(<type 'dict'>, {'LB': {'EX_ser__L_e': -5.0, 
'EX_dcyt_e': -5.0, 'EX_hg2_e': -1000.0, 'EX_ins_e': -5.0, 'EX_cd2_e': 
-1000.0, 'EX_nac_e': -5.0, 'EX_tungs_e': -1000.0, 'EX_glu__L_e': -5.0, 
'EX_trp__L_e': -5.0, 'EX_h_e': -100.0, 'EX_mobd_e': -1000.0, 
'EX_val__L_e': -5.0, 'EX_cobalt2_e': -1000.0, 'EX_so4_e': -1000.0, 
'EX_co2_e': -1000.0, 'EX_k_e': -1000.0, 'EX_cu2_e': -1000.0, 
'EX_zn2_e': -1000.0, 'EX_na1_e': -1000.0, 'EX_cl_e': -1000.0, 
'EX_leu__L_e': -5.0, 'EX_arg__L_e': -5.0, 'EX_pnto__R_e': -5.0, 
'EX_lys__L_e': -5.0, 'EX_ala__L_e': -5.0, 'EX_cbl1_e': -0.01, 
'EX_thr__L_e': -5.0, 'EX_fe3_e': -1000.0, 'EX_adn_e': -5.0, 'EX_pi_e': 
-1000.0, 'EX_thymd_e': -5.0, 'EX_mn2_e': -1000.0, 'EX_phe__L_e': -5.0, 
'EX_ura_e': -5.0, 'EX_dad__2_e': -5.0, 'EX_h2o_e': -100.0, 
'EX_aso3_e': -1000.0, 'EX_hxan_e': -5.0, 'EX_glc__D_e': -5.0, 
'EX_uri_e': -5.0, 'EX_his__L_e': -5.0, 'EX_o2_e': -18.5, 
'EX_pro__L_e': -5.0, 'EX_asp__L_e': -5.0, 'EX_gly_e': -5.0, 
'EX_fe2_e': -1000.0, 'EX_ca2_e': -1000.0, 'EX_mg2_e': -1000.0, 
'EX_cysi__L_e': -5.0, 'EX_tyr__L_e': -5.0, 'EX_met__L_e': -5.0, 
'EX_ile__L_e': -5.0}, 'Host': {'EX_dcyt_e': -5.0, 'EX_ins_e': -5.0, 
'EX_cobalt2_e': -10.0, 'EX_dmso_e': -1000.0, 'EX_so4_e': -1000.0, 
'EX_melib_e': -5.0, 'EX_glcn_e': -5.0, 'EX_uri_e': -5.0, 
'EX_rib__D_e': -5.0, 'EX_tungs_e': -10.0, 'EX_malt_e': -5.0, 
'EX_glu__L_e': -5.0, 'EX_arab__L_e': -5.0, 'EX_12ppd__S_e': -5.0, 
'EX_h_e': -100.0, 'EX_mobd_e': -10.0, 'EX_mnl_e': -5.0, 'EX_mg2_e': 
-10.0, 'EX_acnam_e': -5.0, 'EX_co2_e': -1000.0, 'EX_thm_e': -5.0, 
'EX_glcr_e': -5.0, 'EX_k_e': -10.0, 'EX_cu2_e': -10.0, 'EX_na1_e': 
-1000.0, 'EX_cbl1_e': -0.01, 'EX_cl_e': -1000.0, 'EX_sbt__D_e': -5.0, 
'EX_arg__L_e': -5.0, 'EX_pnto__R_e': -5.0, 'EX_lys__L_e': -5.0, 
'EX_fuc__L_e': -5.0, 'EX_ala__L_e': -5.0, 'EX_csn_e': -5.0, 
'EX_gal_e': -5.0, 'EX_crn_e': -5.0, 'EX_no3_e': -1000.0, 
'EX_thr__L_e': -5.0, 'EX_fe3_e': -10.0, 'EX_adn_e': -5.0, 'EX_pi_e': 
-1000.0, 'EX_glcur_e': -5.0, 'EX_thymd_e': -5.0, 'EX_mn2_e': -10.0, 
'EX_etha_e': -5.0, 'EX_galt_e': -5.0, 'EX_ura_e': -5.0, 'EX_dad__2_e': 
-5.0, 'EX_gsn_e': -5.0, 'EX_h2o_e': -100.0, 'EX_hxan_e': -5.0, 
'EX_glc__D_e': -5.0, 'EX_galct__D_e': -5.0, 'EX_his__L_e': -5.0, 
'EX_met__L_e': -5.0, 'EX_o2_e': -18.5, 'EX_asn__L_e': -5.0, 
'EX_pro__L_e': -5.0, 'EX_acgam_e': -5.0, 'EX_so3_e': -1000.0, 
'EX_cys__L_e': -5.0, 'EX_rmn_e': -5.0, 'EX_ca2_e': -10.0, 
'EX_val__L_e': -5.0, 'EX_man_e': -5.0, 'EX_zn2_e': -10.0, 
'EX_leu__L_e': -5.0, 'EX_fru_e': -5.0, 'EX_fe2_e': -10.0, 
'EX_alltn_e': -5.0, 'EX_galctn__D_e': -5.0, 'EX_ile__L_e': -5.0}, 
'M9': {'EX_cobalt2_e': -1000.0, 'EX_tungs_e': -1000.0, 'EX_h_e': 
-100.0, 'EX_mobd_e': -1000.0, 'EX_mg2_e': -1000.0, 'EX_so4_e': 
-1000.0, 'EX_co2_e': -1000.0, 'EX_k_e': -1000.0, 'EX_fe3_e': -1000.0, 
'EX_na1_e': -1000.0, 'EX_cl_e': -1000.0, 'EX_cu2_e': -1000.0, 
'EX_nh4_e': -1000.0, 'EX_fe2_e': -1000.0, 'EX_cbl1_e': -0.01, 
'EX_pi_e': -1000.0, 'EX_mn2_e': -1000.0, 'EX_h2o_e': -100.0, 
'EX_glc__D_e': -5.0, 'EX_o2_e': -18.5, 'EX_ca2_e': -1000.0, 
'EX_zn2_e': -1000.0}, 'LPM': {'EX_ca2_e': -0.005, 'EX_cit_e': -0.0005, 
'EX_h_e': -100, 'EX_mobd_e': -0.005, 'EX_o2_e': -18.5, 'EX_cobalt2_e': 
-0.005, 'EX_mg2_e': -0.008, 'EX_thm_e': -2.9649835443413292e-08, 
'EX_so4_e': -1.0, 'EX_glyc_e': -41.04680204148116, 'EX_co2_e': -18.5, 
'EX_pi_e': -0.337, 'EX_k_e': -6.0, 'EX_cu2_e': -0.005, 'EX_cl_e': 
-5.016, 'EX_zn2_e': -0.005, 'EX_h2o_e': -1000, 'EX_nh4_e': -15.0, 
'EX_mn2_e': -0.005, 'EX_fe3_e': -0.005}, 'MgM': {'EX_ser__L_e': 
-0.3115, 'EX_cobalt2_e': -0.01, 'EX_glu__L_e': -0.802, 'EX_glyc_e': 
-38.0, 'EX_h_e': -100, 'EX_mobd_e': -0.01, 'EX_val__L_e': -0.3, 
'EX_so4_e': -8.0, 'EX_co2_e': -18.5, 'EX_k_e': -7.002685538757439, 
'EX_fe3_e': -0.01, 'EX_na1_e': -0.551983923305118, 'EX_cl_e': 
-5.645853297492152, 'EX_leu__L_e': -0.4385, 'EX_arg__L_e': -0.1135, 
'EX_nh4_e': -0.3545, 'EX_tyr__L_e': -0.0035, 'EX_lys__L_e': -0.3245, 
'EX_ala__L_e': -0.22, 'EX_thr__L_e': -0.2085, 'EX_pi_e': -1.0, 
'EX_mn2_e': -0.01, 'EX_phe__L_e': -0.0455, 'EX_h2o_e': -1000, 
'EX_his__L_e': -0.097, 'EX_o2_e': -18.5, 'EX_pro__L_e': -0.6185, 
'EX_asp__L_e': -0.3285, 'EX_gly_e': -0.1555, 'EX_cys__L_e': -0.0145, 
'EX_cu2_e': -0.01, 'EX_ca2_e': -0.0020709616248315785, 'EX_mg2_e': 
-0.008098333676198314, 'EX_zn2_e': -0.01, 'EX_met__L_e': -0.1125, 
'EX_ile__L_e': -0.196}}) 

In [4]: eco_model = cobra.test.create_test_model(cobra.test.ecoli_pickle) 

In [5]: eco_model.media_compositions 
--------------------------------------------------------------------------- 
AttributeError                            Traceback (most recent call last) 
<ipython-input-29-5df76f76ba9e> in <module>() 
----> 1 eco_model.media_compositions 

AttributeError: 'Model' object has no attribute 'media_compositions' 

In [6]: 

I am getting a segmentation fault when

from cobra.test.unit_tests import create_test_model, test_pickle

model = create_test_model(test_pickle)
reaction = model.reactions.get_by_id('biomass_iRR1083_metals')
reaction.metabolites

You can find the crash report here

The version of iAF692 exported from Bigg appears to differ from the version published in http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1681478/. The difference is -
boundaryCondition="false" in the former and "true" in the later.

When the "true" conditions are imported using create_cobra_model_from_sbml_file(), the metabolites are not added to the model. If this is incorrect, the code at line 121 needs to be updated.

I am working on yeast7 consensus model right now http://sourceforge.net/projects/yeast/files/
And when I load it with cobrapy and optimize I get very high flux for objective - smt like 800+ while in matlab I get expected (from readme) 0.14.
I guess that not everything from the SBML file was loaded or was loaded correctly - for instance objective function had to be set manually.

I installed the cobra.py and i am able to import that into python. However I could not install the PyGLPK from their website.
I installed it from the synaptic package manager and i am able to import it into pthon. im using ubuntu 12.04. i am new to both ubuntu(linux) and python. please suggest me something. im attaching the screen shot when i ran test_all() function from cobra.

stderr should not get printed out

The following

from cobra import Model, Reaction, Metabolite
from cobra.io import write_sbml_model
cobra_model = Model('example_cobra_model')  # Best practise: SBML compliant IDs

reaction = Reaction('3OAS140')
reaction.name = '3 oxoacyl acyl carrier protein synthase n C140 '
reaction.subsystem = 'Cell Envelope Biosynthesis'
reaction.lower_bound = -100.  # This is the default
reaction.upper_bound = -10.  # This is the default
reaction.objective_coefficient = 0. # this is the default

# Create the metabolites
ACP_c = Metabolite('ACP_c', formula='C11H21N2O7PRS',
    name='acyl-carrier-protein', compartment='c')
omrsACP_c = Metabolite('3omrsACP_c', formula='C25H45N2O9PRS',
    name='3-Oxotetradecanoyl-acyl-carrier-protein', compartment='c')
co2_c = Metabolite('co2_c', formula='CO2', name='CO2', compartment='c')
malACP_c = Metabolite('malACP_c', formula='C14H22N2O10PRS',
    name='Malonyl-acyl-carrier-protein', compartment='c')
h_c = Metabolite('h_c', formula='H', name='H', compartment='c')
ddcaACP_c = Metabolite('ddcaACP_c', formula='C23H43N2O8PRS',
    name='Dodecanoyl-ACP-n-C120ACP', compartment='c') 

reaction.add_metabolites({malACP_c: -1.0,
                          h_c: -1.0,
                          ddcaACP_c: -1.0,
                          co2_c: 1.0,
                          ACP_c: 1.0,
                          omrsACP_c: 1.0})

# A string representation of the reaction
print(reaction.reaction)

cobra_model.add_reaction(reaction)
print cobra_model.reactions[0].reversibility
write_sbml_model(cobra_model, 'reaction_reversibility_in_SBML_fail.xml')
print open('reaction_reversibility_in_SBML_fail.xml', 'r').read()

will write the reaction in the wrong direction in the resulting SBML file (malACP_c + ddcaACP_c + h_c --> co2_c + ACP_c + 3omrsACP_c).

<reaction id="R_3OAS140" name="3 oxoacyl acyl carrier protein synthase n C140 " reversible="false" fast="false">
<notes>
  <html xmlns="http://www.w3.org/1999/xhtml">
    <p>GENE_ASSOCIATION: </p>
    <p>SUBSYSTEM: Cell Envelope Biosynthesis</p>
  </html>
</notes>
<listOfReactants>
  <speciesReference species="M_malACP_c" stoichiometry="1" constant="true"/>
  <speciesReference species="M_ddcaACP_c" stoichiometry="1" constant="true"/>
  <speciesReference species="M_h_c" stoichiometry="1" constant="true"/>
</listOfReactants>
<listOfProducts>
  <speciesReference species="M_co2_c" stoichiometry="1" constant="true"/>
  <speciesReference species="M_ACP_c" stoichiometry="1" constant="true"/>
  <speciesReference species="M_3omrsACP_c" stoichiometry="1" constant="true"/>
</listOfProducts>
</reaction>

It should be co2_c + ACP_c + 3omrsACP_c --> malACP_c + ddcaACP_c + h_c as irreversible SBML reactions can only carry positive fluxes.

However, labeling a reaction as irreversible is interpreted as an assertion
that the rate expression will not have negative values during a simulation.

The Systems Biology Markup Language (SBML): Language Specification for Level 3 Version 1 Core

Hi,
I tried to install the latest version via pip and I get a compiler error on the cython file. When I run cobra after this installation, only glpk is in cobra.solvers (not cglpk).
Downloading cobra-0.3.0.tar.gz (2.3Mb): 2.3Mb downloaded
Running setup.py egg_info for package cobra

Error compiling Cython file:
------------------------------------------------------------
...
# Do not want to print out to terminal. Even when not verbose, output
# will be redirected through the hook.
glp_term_out(GLP_OFF)


cdef int hook(void *info, const char *s):
                                    ^
------------------------------------------------------------

cobra/solvers/cglpk.pyx:92:37: Expected ')', found '*'
Compiling cobra/solvers/cglpk.pyx because it changed.
Cythonizing cobra/solvers/cglpk.pyx

Do you also see this problem? I was getting this error before, but I needed QSOpt to work, so I had written my own interface to that solver using the regular glpk bindings.

I have noticed some strange behaviour when adding a reaction that
contains the same metabolite on both sides:

from cobra import Model, Reaction, Metabolite

model = Model("test")

A = Metabolite("A", formula="C1", name="Substance A", compartment="c")
B = Metabolite("B", formula="C2", name="Substance B", compartment="c")
C = Metabolite("C", formula="C3", name="Substance C", compartment="c")

reaction = Reaction("new_reaction")
reaction.name = "Test reaction"

reaction.add_metabolites({A: -1.0,
                          B: -1.0,
                          A: 1.0,
                          C: 1.0})

print reaction.reaction
model.add_reaction(reaction)

output is

B --> C + A

instead of

A + B --> C + A

I'm working with cobra 0.3.0.dev.

Hi,

I think I found two bugs in the assess_objective function:

line 48:

sink_reaction = Reaction('test_sink_' + the_component)

throws a TypeError because it tries to concatenate a string and a metabolite object,

sink_reaction = Reaction('test_sink_' + str(the_component))

solves it.

line 53:

sink_reaction.add_metabolites(the_component, tmp_coefficient)

throws an AttributeError, because the argument has to be a dict.

sink_reaction.add_metabolites({the_component: tmp_coefficient})

solves it.

Hope it helps.

Kind Regards,
Michael

The version of iAF692 exported from Bigg appears to differ from the version published in http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1681478/. The difference is -
boundaryCondition="false" in the former and "true" in the later.

When the "true" conditions are imported using create_cobra_model_from_sbml_file(), the metabolites are not added to the model. If this is incorrect, the code at line 121 needs to be updated.

For example, to determine if a solver supports QP, test for the availability of the function quadratic_component, and to test for MIP, test for a boolean flag such as MIP_AVAILABLE

When no solver is passed in, either gurobi or cplex should be used (if available) instead of glpk if a quadratic objective is given

cobra.io.save_json_model() cannot handle most numpy objects. This pops up in a confusing way if you try to serialize reaction reversibilities. Turns out the following returns numpy.bool_ if the bounds are numpy integers:

Reaction.py@131

        return self.lower_bound < 0 and self.upper_bound > 0

I guess the only sensible solution is to subclass the JSON encoder to handle the numpy cases.

In manipulation/modify.py module, when converting irreversible model to reversible model by the fuction "revert_to_reversible(cobra_model)", I think the cobra_model.solution (if availible) should be also updated. Something like this worked for me:

if not cobra_model.solution is None:
for reverse in reverse_reactions:
forward = reverse.reflection
cobra_model.solution.x_dict[forward.id] -= /cobra_model.solution.x_dict[reverse.id]
del cobra_model.solution.x_dict[reverse.id]

This patch is usefull to me to eliminate "futile cycle" in the solution after computing pFBA.

Hope it helps
cheers
miguel

Why is

pickle.loads(pickle.dumps(iJO))  # (1 loops, best of 3: 1.88 s per loop)

is slower than

json_model = from_json(to_json(iJO))  # (1 loops, best of 3: 380 ms per loop)

Could one define __getstate__ and __setstate__ in cobra.core.Model.Model to use to_json and from_json? Speeding up serialization will help a lot with parallelizing cobra calculations.

The code below shows the error I get when I try to convert a model to irreversible and then revert it back.

In [76]: model_sbml = cobra.io.read_sbml_model("yeast_7.11_cobra.xml")

In [77]: len(model_sbml.reactions)
Out[77]: 3490

In [78]: modify.convert_to_irreversible(model_sbml)

In [79]: len(model_sbml.reactions)
Out[79]: 4813

In [80]: modify.revert_to_reversible(model_sbml)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-80-92bae93b42f2> in <module>()
----> 1 modify.revert_to_reversible(model_sbml)

/home/melomcr/Research/cobrapy/cobra/manipulation/modify.pyc in revert_to_reversible(cobra_model, update_solution)
    155         for reverse in reverse_reactions:
    156             forward = reverse.reflection
--> 157             x_dict[forward.id] -= x_dict.pop(reverse.id)
    158         cobra_model.solution.x = [x_dict[r.id] for r in cobra_model.reactions]
    159 

TypeError: 'NoneType' object has no attribute '__getitem__'

That portion of the code should only be called if the revert_to_reversible function was called with update_solution=True (which is the default), AND if the Solution is not None.
It seems the problem is that "Solution" is never evaluated as "None", even if the model has just been loaded, as in the example below. Since the object is there, it won't be "None". That is problem number (1).

In [94]: model_sbml = cobra.io.read_sbml_model("yeast_7.11_cobra.xml")

In [95]: model_sbml.solution is None
Out[95]: False

To solve problem (1), the check for a solution should be changed to

cobra_model.solution.status is not "NA" 

This solves the problem in my fork.

Even if I do optimize the model and have a valid solution, the function still crashes, which is problem number (2). Since they are both in the same function and give the same output, I decided to place them both here, in the same issue.
If I convert to irreversible first then optimize and then revert to reversible, everything runs ok.
If I optimize first then convert to irreversible and then revert to reversible, I get the following error.

In [130]: model_sbml = cobra.io.read_sbml_model("yeast_7.11_cobra.xml")

In [131]: modify.convert_to_irreversible(model_sbml)

In [132]: model_sbml.optimize(solver="gurobi")
Out[132]: <Solution 0.07 at 0x1a12cc10>

In [133]: modify.revert_to_reversible(model_sbml)

In [134]: model_sbml = cobra.io.read_sbml_model("yeast_7.11_cobra.xml")

In [135]: model_sbml.optimize(solver="gurobi")
Out[135]: <Solution 0.07 at 0x3325ac90>

In [136]: modify.convert_to_irreversible(model_sbml)

In [137]: modify.revert_to_reversible(model_sbml)
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-137-92bae93b42f2> in <module>()
----> 1 modify.revert_to_reversible(model_sbml)

/home/melomcr/Research/cobrapy/cobra/manipulation/modify.pyc in revert_to_reversible(cobra_model, update_solution)
    155         for reverse in reverse_reactions:
    156             forward = reverse.reflection
--> 157             x_dict[forward.id] -= x_dict.pop(reverse.id)
    158         cobra_model.solution.x = [x_dict[r.id] for r in cobra_model.reactions]
    159 

KeyError: 'r_0003_reverse'

In [138]:

I am not sure, however, what is the problem with the line 157. I did check and the 'r_0003_reverse' reaction is there. Does anyone have any ideas?

All tests were done using Yeast model 7.11 that can be found here.

Ali,

I recall you had code or wanted code that'd dress the model variables with the simulation results after a simulation to make it easier to access results; e.g., reaction.x vs model.solution.x_dict[reaction.id]. My main concern about the dressing procedure was that if done at the wrong time it could induce an unnecessary performance penalty in instances where one's running many simulations and doesn't need all of the results unless the procedure was optimized for the different solvers.

I was thinking about adding a function to core.Solution that will effectively dress the variables and was wondering if you've already implemented this or if i should do so?

def dress_results(self):
[setattr(k, 'x', v) for k, v in zip(self.reactions, self.x])];
[setattr(k, 'y', v) for k, v in zip(self.metabolites, self.y])];

hi everyone,

I have noticed a small glitch between the implementations of the cobra toolbox in python and Matlab. After removing several genes, metabolites and reactions from a yeast model ( always taking into account that it should be possible to create biomass ) I am able to get a solution (FBA) in python but not in Matlab.

What I do more specifically:
I take a standard yeast model (SBML file) and import it into Matlab. I modify the lower bounds to match my medium and remove all dead end reactions using the standard functions.
next, I export this model as an SBML file and import it into python.
And then I use this model to determine which reaction should be removed, and remove them using the code below.
Next, I perform an FBA on both the model before deleting the reactions, metabolites and genes and after deleting. Finally,I write out an SBML file

# Original model
model = cobra.io.sbml.create_cobra_model_from_sbml_file(os.path.join('data','iMM904_NADcorrected_pappEtAl2011_adj_and_remove.xml'))
model.optimize(solver='glpk')
print 'optimal f full model '+str(model.solution.f)

# Remove
for reaction in removed_reactions: model.reactions.__getattr__(reaction[2:]).remove_from_model(model)
for gene in removed_genes: model.genes.__getattr__(gene).remove_from_model(model)
for metabolite in removed_metabolites: model.metabolites.__getattr__(metabolite[2:]).remove_from_model(model)

# New model calculations
model.optimize(solver='glpk')
print 'optimal f partial model ' + str(model.solution.f)
FVA_result = cobra.flux_analysis.variability.flux_variability_analysis(model, fraction_of_optimum=1.0, number_of_processes=15)

# Write to file
write_cobra_model_to_sbml_file(model, os.path.join('data','iMM904_examo_%s.xml' % description))

Outputs:
optimal f full model 0.287998303658
optimal f partial model 0.0635426076692

as you can see both the full model and the model in which I deleted some reactions, metabolites and genes are able to create biomass. The partial model is only able to do it less efficiently.

I import the SBML file in Matlab using the code below:

initCobraToolbox

reduced_model = readCbModel(['C:\Master_CNCZ_new\trunk\CMBI\src\EXAMO_scripts\data\iMM904_examo_WT_MDM34.xml']);
old_model = readCbModel(['C:\Master_CNCZ_new\trunk\CMBI\src\EXAMO_scripts\data\iMM904_NADcorrected_pappEtAl2011_adj_and_remove.xml']);

FBAsolution_old = optimizeCbModel(old_model);
FBAsolution_red = optimizeCbModel(reduced_model);
FBAsolution_old.f
FBAsolution_red.f

Outputs:
ans =
0.2880
ans =
0

as you can see the partial model gives the exact same result in both python in Matlab. However, the model in which I deleted some genes, reactions and metabolites is not able to generate a solution in Matlab. If I try to do an FVA, it generates an error in Matlab but works fine in python.

Am I doing something wrong in exporting the SBML file from python, or could there be some other cause?

Any help would be greatly appreciated,

cheers,
Rob

here's what I'm using:
Matlab R2013a
Python 2.7
cobra.version: 0.2.1

Initilization Matlab:
Warning: LP solver Gurobi not useable: gurobi.m not in Matlab path

In changeCobraSolver at 138
In initCobraToolbox at 42
In fa_fba_analysis at 1
Warning: LP solver Gurobi not useable: gurobi_mex not in Matlab path
In changeCobraSolver at 133
In initCobraToolbox at 42
In fa_fba_analysis at 1
Warning: LP solver CPLEX through Tomlab not usable: tomRun.m not in Matlab path
In changeCobraSolver at 113
In initCobraToolbox at 42
In fa_fba_analysis at 1
LP solver set to glpk successful
Warning: MILP solver Gurobi not useable: gurobi.m not in Matlab path
In changeCobraSolver at 175
In initCobraToolbox at 55
In fa_fba_analysis at 1
Warning: MILP solver Gurobi not useable: gurobi_mex not in Matlab path
In changeCobraSolver at 170
In initCobraToolbox at 55
In fa_fba_analysis at 1
Warning: MILP solver CPLEX through Tomlab not usable: tomRun.m not in Matlab path
In changeCobraSolver at 160
In initCobraToolbox at 55
In fa_fba_analysis at 1
MILP solver set to glpk successful
Warning: QP solver Gurobi not useable: gurobi.m not in Matlab path
In changeCobraSolver at 231
In initCobraToolbox at 67
In fa_fba_analysis at 1
Warning: QP solver Gurobi not useable: gurobi_mex not in Matlab path
In changeCobraSolver at 224
In initCobraToolbox at 67
In fa_fba_analysis at 1
Warning: QP solver CPLEX through Tomlab not usable: tomRun.m not in Matlab path
In changeCobraSolver at 195
In initCobraToolbox at 67
In fa_fba_analysis at 1
Warning: qpng solver has not been fully tested - results may not be correct
In changeCobraSolver at 205
In initCobraToolbox at 67
In fa_fba_analysis at 1
QP solver set to qpng successful
Warning: MIQP solver gurobi not usable: gurobi.m not in Matlab path
In changeCobraSolver at 262
In initCobraToolbox at 78
In fa_fba_analysis at 1
Warning: MIQP solver gurobi not usable: gurobi_mex not in Matlab path
In changeCobraSolver at 255
In initCobraToolbox at 78
In fa_fba_analysis at 1
Warning: MIQP solver CPLEX through Tomlab not usable: tomRun.m not in Matlab path
In changeCobraSolver at 248
In initCobraToolbox at 78
In fa_fba_analysis at 1
MIQP solver set failed
CB map output set to svg successful

It would be nice to have functions to rename reaction ids and names.

Currently, if the user changes the reaction name by doing reaction.id = new_id, the old_id stays in the dictlist.

The user has to know this and manually update the index with:
cobra_model.reactions._generate_index() . One alternative to adding functions is to improve documentation on the Reaction object.

I can do a pull request if people think either is a good idea.

Josh

I'm trying to get set up with the Yeast model, but it seems the yeast models do not produce biomass when simulating using cobrapy. Its probably something stupid.

I got the latest yeast model from http://sourceforge.net/projects/yeast/files/yeast_7.1.zip/download

(I tried others as well)

then I ran:
from cobra.io.sbml import create_cobra_model_from_sbml_file

yeast_model = 'yeast_7.1_cobra.xml'
cobra_model = create_cobra_model_from_sbml_file(yeast_model, print_time=True)
cobra_model.optimize(solver='glpk')
print '\nSimulated growth rate is %1.3f' % cobra_model.solution.f

and got the output:
"Simulated growth rate is nan"

Their readme states that the answer should have been:
maximal growth rate: 0.13975 flux units

Any ideas? If not, I'll dig into it in a bit. Don't have a MATLAB license to check the answer using the original cobra toolbox.

Josh

Double deletion should take keyword arguments for solver parameters. This will require modifying the deletion workers.

The double deletion should also first run all single deletions and then use those to avoid submitting to the workers pairs for which one of the deletions has already been computed to be lethal.

Hi,

moma does not run on python3:

--> 226 for reaction_1, reaction_2 in common_dict.iteritems():

seems like iteritems() has been renamed to items()

Best regards

In the latest version of COBRA (0.3.0), cobra.flux_analysis.variability.find_blocked_reactions() fails. I get the following error message:

  File "/usr/local/lib/python2.7/dist-packages/cobra/flux_analysis/variability.py", line 81, in find_blocked_reactions
    if not the_reactions:
UnboundLocalError: local variable 'the_reactions' referenced before assignment

Passing the_reactions to the function as a keyword argument does not resolve the problem. It appears that the name of the_reactions has been changed to reaction_list without fully updating the body of the function.

OS X 10.9.2
Mac Pro 2014
Apple LLVM version 5.1 (clang-503.0.38) (based on LLVM 3.4svn)
Target: x86_64-apple-darwin13.1.0

cc -fno-strict-aliasing -fno-common -dynamic -arch x86_64 -arch i386 -g -Os -pipe -fno-common -fno-strict-aliasing -fwrapv -mno-fused-madd -DENABLE_DTRACE -DMACOSX -DNDEBUG -Wall -Wstrict-prototypes -Wshorten-64-to-32 -DNDEBUG -g -fwrapv -Os -Wall -Wstrict-prototypes -DENABLE_DTRACE -arch x86_64 -arch i386 -pipe -I/System/Library/Frameworks/Python.framework/Versions/2.7/include/python2.7 -c cobra/solvers/cglpk.c -o build/temp.macosx-10.9-intel-2.7/cobra/solvers/cglpk.o

clang: error: unknown argument: '-mno-fused-madd' [-Wunused-command-line-argument-hard-error-in-future]

clang: note: this will be a hard error (cannot be downgraded to a warning) in the future

simple workaround
export CFLAGS=-Qunused-arguments
export CPPFLAGS=-Qunused-arguments
before building

Currently, DictList cannot be serialized with dill. Consequently, cobrapy models etc. cannot be serialized with dill. Of course this is not a critical issue. It would just be nice to have support for dill as it can serialize many things (like closures, partials etc.) that pickle can't.

The following code should reproduce the problem (I am currently on cobrapy version 0.3.0.dev-130-g3776829, OS X 10.9.1, Python 2.7.5).

import pickle
import dill

# DictList seems to be undillable
from cobra.core.DictList import DictList
from cobra.core.Object import Object

dict_list = DictList()
dict_list.append(Object('test1'))

# Works
repickled = pickle.loads(pickle.dumps(dict_list))
print repickled.get_by_id('test1')

# Does not work
redilled= dill.loads(dill.dumps(dict_list))
print redilled.get_by_id('test1')

It produces the the following exception

Traceback (most recent call last):
  File "/Users/niko/Arbejder/Dev/cobrapy/cobra/oven/phantomas1234/dill_support.py", line 27, in <module>
    redilled= dill.loads(dill.dumps(dict_list))
  File "/Users/niko/.virtualenvs/meteng/lib/python2.7/site-packages/dill/dill.py", line 135, in loads
    return load(file)
  File "/Users/niko/.virtualenvs/meteng/lib/python2.7/site-packages/dill/dill.py", line 128, in load
    obj = pik.load()
  File "/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/pickle.py", line 858, in load
    dispatch[key](self)
  File "/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/pickle.py", line 1182, in load_append
    list.append(value)
  File "/Users/niko/Arbejder/Dev/cobrapy/cobra/core/DictList.py", line 111, in append
    self._check(the_id)
  File "/Users/niko/Arbejder/Dev/cobrapy/cobra/core/DictList.py", line 30, in _check
    if id in self._dict:
  File "/Users/niko/Arbejder/Dev/cobrapy/cobra/core/DictList.py", line 271, in __getattr__
    (attr))
AttributeError: DictList has no attribute or entry _dict

Hi,
I'm trying to read an SBML file using a python script that calls cobrapy. I'm getting the following error:

/usr/local/lib/python2.7/dist-packages/cobra-0.2.0b1-py2.7.egg/cobra/io/init.py:8: UserWarning: cobra.io.sbml will not be functional: ImportError dynamic module does not define init function (initlibsbml)
warn("cobra.io.sbml will not be functional: ImportError %s" % error)
Traceback (most recent call last):
File "ReadSBML.py", line 11, in
from cobra.io.sbml import create_cobra_model_from_sbml_file
File "/usr/local/lib/python2.7/dist-packages/cobra-0.2.0b1-py2.7.egg/cobra/io/sbml.py", line 24, in
from libsbml import SBMLDocument, SpeciesReference, KineticLaw, Parameter
ImportError: dynamic module does not define init function (initlibsbml)

What can I do to solve this?

Thanks!

Often I run into a case where I have a list of cobrapy objects, e.g. reactions, and the repr looks like this:

['FRD2', 'LDH_D2', 'PTA2', 'FRD3']

This is identical to a list of strings with reaction ids. So, let's say I want to find a unique set of reactions. I will see the following strange behavior:

print a
['FRD2', 'LDH_D2', 'PTA2', 'FRD3']
print b
['PTAr', 'LDH_D', 'LDH_D2', 'FRD2', 'FRD3', 'PTA2']
print set(b).difference(a)
['PTAr', 'LDH_D', 'LDH_D2', 'FRD2', 'FRD3', 'PTA2']

In this case, a is a list of Reaction objects and b is a list of reaction_id strings. A solution to this would be something like the following for a reaction repr:

[model.reactions.FRD2,
 model.reactions.LDH_D2, 
 model.reactions.PTA2,
 model.reactions.FRD3]

The documentation supports this kind of display: https://docs.python.org/2/reference/datamodel.html#object.__repr__

The following demonstrates

import cobra.test
model1 = cobra.io.read_sbml_model(cobra.test.salmonella_sbml)
model = cobra.Model("test")
model += model1

This bug probably effects other functions too.

The following commits exhibit the bug:
816271d
ab9d295

The following commits work fine:
a824f8a (speedup branch)
82d580a (0.2.1)
01ec76b

I think this was not caught in unit tests because the pickles were not being loaded properly.

When retrieving the genes associated with a reaction, the frozenset contains the genes 'AND' and 'OR'. In the source code for the class Reaction there is this:

def gene_reaction_rule(self, new_rule):
self._gene_reaction_rule = new_rule
gene_names = set((re.compile(' {2,}').sub(' ', re.compile('(| and| or|+|)').sub('', self._gene_reaction_rule))).split(' ' ))

which clearly doesn't separate the upper case bolean rules. Changing the regular expression from
'(| and| or|+|)'
to
'(| and| or| AND| OR|+|)'
Would solve the issue.

You can recreate this bug easily:

import cobra
import cobra.io
import cobra.test
m = cobra.io.load_matlab_model(cobra.test.ecoli_mat)
m.genes.get_by_id('b0432').remove_from_model(m)

And the trace:

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-5-591696bb1aa3> in <module>()
----> 1 m.genes.get_by_id('b0432').remove_from_model(m)

/Users/zaking/sharedrepos/cobrapy_master/cobra/core/Gene.pyc in remove_from_model(self, model, make_dependent_reactions_nonfunctional)
     74             the_reaction.gene_reaction_rule = the_gene_re.sub(gene_state,
     75                                                               the_reaction.gene_reaction_rule)
---> 76             the_reaction._genes.remove(self)
     77             #Now deactivate the reaction if its gene association evaluates to False
     78             the_gene_reaction_relation = the_reaction.gene_reaction_rule

KeyError: <Gene b0432 at 0x114d72d90>

Somewhere in the last set up updates MOMA appears to have stopped functioning.
from cobra.test import create_test_model

model = create_test_model()
from cobra.flux_analysis import single_deletion
single_deletion(model, ['STM0247'], method='moma', solver='gurobi')
single_deletion(model, ['STM0247'], method='moma', solver='gurobi')

cplex 12.5
gurobi 5.5
0.3.0.dev-118-g09707d2
OS X 10.9.1
Python 2.7.5 (default, Aug 25 2013, 00:04:04)
[GCC 4.2.1 Compatible Apple LLVM 5.0 (clang-500.0.68)] on darwin

Could somebody confirm that the failure happens on another OS?

To ensure the long-term viability of exported JSON models, it would be good to define a schema for them. json-schema lets you define these easily, and validators are available in a number of languages, including Python and JavaScript.

Escher will use json-schema for defining the Escher map file format. We can also consider .escher and .cobra filename extensions to clarify which file is which.

I can help out with this.

When running a very large number of simulations using different upper bounds for reactions in the model, some times GLPK cannot find a solution. In the rare event is happens (one every few hundred thousand cases) cobrapy hangs and crashes some time later. It took me a while to find the issue because the only error output would be from MPI saying one of the ranks terminated unexpectedly. Using C-GLPK, this rare event produces an exception (the output is copied below), so using a try-except I could treat such cases. I am not familiar with the inner workings of GLPK, but I believe some similar treatment should be creaded to handle such cases.

Traceback (most recent call last):
File "/home/melomcr/cobrapy/cobrapy_mpi_v2.py", line 650, in
resultVmaxGR = prepare_and_optimize_model(SimPerCore, model, proteins, v_info_dict)
File "/home/melomcr/cobrapy/cobrapy_mpi_v2.py", line 315, in prepare_and_optimize_model
localModel.optimize(solver='cglpk') # 'gurobi' or 'cglpk' or 'glpk'
File "/home/melomcr/Programas/cobrapy-master/cobra/core/Model.py", line 301, in optimize
*_kwargs)
File "/home/melomcr/Programas/cobrapy-master/cobra/solvers/init.py", line 110, in optimize
the_solution = solver_function.solve(cobra_model, *_kwargs)
File "cglpk.pyx", line 186, in cobra.solvers.cglpk.GLP.solve (cobra/solvers/cglpk.c:3392)
File "cglpk.pyx", line 170, in cobra.solvers.cglpk.GLP.solve_problem (cobra/solvers/cglpk.c:3155)
AssertionError

Hi,
I have a problem translating a piece of COBRA code from matlab to python. I believe it's rather trivial but honestly I'm a bit confused with the result I obtain:

The matlab code is taken from the example 4 of the supplementary material of this paper: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3108565/
The xml file can be found here:
http://systemsbiology.ucsd.edu/sites/default/files/Attachments/Images/downloads/Ecoli_core/ecoli_core_model.xml

Matlab code:

model = readCbModel('ecoli_core_model.xml')
model = changeObjective(model,'Biomass_Ecoli_core_w_GAM');
FBAsolution = optimizeCbModel(model,'max');
fprintf('Flux Biomass : %u\n',FBAsolution.f);

model = changeRxnBounds(model,'EX_o2(e)',-17,'b');
growthRates = zeros(21,1);
for i = 0:20
model = changeRxnBounds(model,'EX_glc(e)',-i,'b');
FBAsolution = optimizeCbModel(model,'max');
growthRates(i+1) = FBAsolution.f;
end
disp(growthRates)
plot(growthRates)

% This works properly

Python code:

from cobra.io.sbml import create_cobra_model_from_sbml_file

solver = 'glpk' #Change to 'gurobi' or 'cplex'
sbml_file = 'ecoli_core_model.xml'
obj_function = 'Biomass_Ecoli_core_w_GAM'

cobra_model = create_cobra_model_from_sbml_file(sbml_file, print_time=True)

Set the objective function

obj_rxn = cobra_model.reactions.get_by_id(obj_function)
obj_rxn.objective_coefficient = 1
cobra_model.optimize(solver=solver)

model = changeRxnBounds(model,'EX_o2(e)',-17,'b');

R1 = cobra_model.reactions.get_by_id("EX_o2_e")
R1.lower_bound = -17
R1.upper_bound = -17

growthRates = zeros(21,1);

growthRates = ()
R2 = cobra_model.reactions.get_by_id("EX_glc_e")

for i = 0:20

for i in range(0,20):
R2.lower_bound = -i
R2.upper_bound = -i
cobra_model.optimize(solver=solver)
growthRates.append(cobra_model.solution.f)

print growthRates

-> This, however doesn't work, the model cannot be optimized.

Can you also reproduced this ?
What I am doing wrong here ?

I'd really appreciate a little bit of help.
Thanks a lot
Regards,
Michael

https://github.com/opencobra/cobrapy/blob/master/cobra/manipulation/delete.py#L94 fails for reactions that contain gene IDs that overlap. For example,

(AV727381.1) or (1537.1) or (H94367.1) or (29796.1) or (10975.1) or (7381.1) or (7384.1) or (7385.1) or (7386.1) or (7388.1) or (27089.1)

will evaluate to

(AV72True) or (True) or (True) or (True) or (True) or (True) or (True) or (True) or (False) or (True) or (True)

because of 7381.1 is a substring of AV727381.1 and eval will choke on that.

Instead of using string based hacks, I would suggest using something more robust (like sympy) to capture the logic of the GPRs. For example, I've used the following in the past:

gpr_rule = sympify(rule.replace('or', '|').replace('and', '&'))
gpr_rule.subs({gene1: True, gene2: False, ...})

Using sympy would also have the benefit of being able to do things like

gpr_rule.subs({Or: Max, And: Min})

if you want to translate gene expression values in reaction-based expression values. Maybe there is a more specialized logic package out there so one wouldn't have to install complete sympy.

Posted in this discussion thread: https://groups.google.com/forum/#!topic/cobra-pie/-rMDWuSpZy0

If an already irreversible model is passed in, the optimize_minimal_flux fails.

There are really two issues that cause this:

• Even if already_irreversible is True, the functions attempts to run revert_to_reversible
• model.remove_reactions fails when given an empty list

Perhaps you could add
python-dev
as neccessary packages on debian-like systems in the install instructions (1)?
Otherwise pip install cobra --pre does not compile:
cobra/solvers/cglpk.c:8:22: fatal error: pyconfig.h: No such file or directory

The following sample code demonstrates the bug

import cobra.test
model = cobra.test.create_test_model()
pgi = model.reactions.get_by_id("PGI")
print "before deletion:", pgi.reaction
model.remove_reactions([pgi])
print "after deletion:", pgi.reaction

outputs:

before deletion: 1.0 g6p_c <=> 1.0 f6p_c
after deletion: 2.0 g6p_c <=> 2.0 f6p_c

Hi there,

Lately I've been working with some macromolecules that were used as pools for several species, and that have formulas that look like:

C36.90H69.37N1.66O4.50PS

As the element counts are non-integers, the formulas are not parsed and an error is thrown.
I just wanted to propose a workaround I rigged to avoid that, if you are interested. It is just a modification of the matching regex and a tweak in the function core.Formula.parse_composition():

element_re = re.compile("([A-Z][a-z]?)([0-9.]+[0-9.]?|(?=[A-Z])?)"

def parse_composition(self):
        """Breaks the chemical formula down by element."""
        tmp_formula = f.replace("*", "")
        composition = {}
        parsed = element_re.findall(tmp_formula)
        for (element, count) in parsed:
            if count == '':
                count = 1
            else:
                try:
                    count = float(count)
                except:
                    warn("failed to parse %s (in formula %s)" % (count, self.formula))
                    self.elements = {}
                    return
            if element in composition:
                composition[element] += count
            else:
                composition[element] = count
        self.elements = composition

the regex, for the input above, returns:

[('C', '36.90'),
 ('H', '69.37'),
 ('N', '1.66'),
 ('O', '4.50'),
 ('P', ''),
 ('S', '')]

Anyway let me know if this is useful.

7 characters is enough to get the commit but 4 is not

If I run flux_variability_analysis on a model that does not have a solution it fails with the following error.

/cobrapy/cobra/flux_analysis/variability.pyc in flux_variability_analysis(cobra_model, reaction_list, fraction_of_optimum, solver, objective_sense, **solver_args)
39 for i, r in enumerate(cobra_model.reactions):
40 if r.objective_coefficient != 0:
---> 41 f = solution.x_dict[r.id]
42 new_bounds = (f * fraction_of_optimum, f)
43 solver.change_variable_bounds(lp, i, min(new_bounds), max(new_bounds))

TypeError: 'NoneType' object has no attribute 'getitem'

Building from source fails without this

When trying to set the reversibility of my reactions I get following error:

reaction.reversibility = True
AttributeError: can't set attribute

My code is listed below. I have a dictionary (Reactions) with all the neccesary information where I loop trough.

Can somebody help me with this?

Many Thanks,

Bart

for i in Reactions:
reaction = Reaction(Reactions[i][0])
reaction.reversibility = True
reaction.name= Reactions[i][0]
reaction.lower_bound = Reactions[i][7][0]
reaction.upper_bound = Reactions[i][7][1]
for q in range(len(Reactions[i][5][0])):
reaction.add_metabolites({cobra_model.metabolites.get_by_id(Reactions[i][5][0][q]):int(Reactions[i][5][1][q])})
for q in range(len(Reactions[i][6][0])):
reaction.add_metabolites({cobra_model.metabolites.get_by_id(Reactions[i][6][0][q]):int(Reactions[i][6][1][q])})
cobra_model.add_reaction(reaction)

In the latest master, I get this error when loading the RECON1 SBML file from BIGG:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-3-59ccb4b31860> in <module>()
----> 1 cobra.io.read_sbml_model('RECON1.xml')

/Users/zaking/sharedrepos/cobrapy_master/cobra/io/sbml.pyc in create_cobra_model_from_sbml_file(sbml_filename, old_sbml, legacy_metabolite, print_time, use_hyphens)
    262                         except ValueError, e:
    263                             print the_item
--> 264                             raise e
    265                     tmp_locus_id = tmp_row_dict['LOCUS']
    266                     if 'TRANSCRIPT' in tmp_row_dict:

ValueError: too many values to unpack

I see the same error using read_legacy_sbml.

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-6-b1b24138809c> in <module>()
----> 1 cobra.io.read_legacy_sbml('RECON1.xml')

/Users/zaking/sharedrepos/cobrapy_master/cobra/io/sbml.pyc in read_legacy_sbml(filename, use_hyphens)
    558 def read_legacy_sbml(filename, use_hyphens=False):
    559     """read in an sbml file and fix the sbml id's"""
--> 560     model = create_cobra_model_from_sbml_file(filename)
    561     for metabolite in model.metabolites:
    562         metabolite.id = fix_legacy_id(metabolite.id)

/Users/zaking/sharedrepos/cobrapy_master/cobra/io/sbml.pyc in create_cobra_model_from_sbml_file(sbml_filename, old_sbml, legacy_metabolite, print_time, use_hyphens)
    262                         except ValueError, e:
    263                             print the_item
--> 264                             raise e
    265                     tmp_locus_id = tmp_row_dict['LOCUS']
    266                     if 'TRANSCRIPT' in tmp_row_dict:

ValueError: too many values to unpack

With the current deepcopy function for cobra.core.DictList, the relational references between metabolites and genes and their respective reactions are no longer maintained in the Model after performing a deepcopy:

from cobra.test import create_test_model
from copy import deepcopy
a = create_test_model()
b = deepcopy(a)
metabolite_a = a.metabolites[0]
metabolite_b = b.metabolites.get_by_id(metabolite_a.id)
[x._reaction for x in [metabolite_a, metabolite_b]]

Out[7]: [set(['3OAS120', 'FA100ACPHi', 'AACPS8', 'EAR100y', 'EAR100x']), set([])]

Reverting to the earlier code from the 0.2.0 release fixes the problem.

This bug is likely due to how the states for cobra.Objects are managed to speed up pickling / unpickling times.