12. Using COBRApy with MASSpy

This notebook example demonstrates how to convert COBRApy objects into their equivalent MASSpy objects, and highlights some of the differences between them.

import cobra.test

from mass import MassMetabolite, MassModel, MassReaction

12.1. Converting COBRA to MASS

Converting COBRApy objects into their MASSpy equivalents is a simple process. It only requires the user to instantiate the MASSpy object using the COBRApy object.

# Get some COBRA objects
cobra_model = cobra.test.create_test_model("textbook")
cobra_metabolite = cobra_model.metabolites.get_by_id("atp_c")
cobra_reaction = cobra_model.reactions.get_by_id("PGI")

Set parameter Username

12.1.1. Metabolite to MassMetabolite

To convert a cobra.Metabolite to a mass.MassMetabolite:

mass_metabolite = MassMetabolite(cobra_metabolite)
mass_metabolite

MassMetabolite identifier	atp_c
Name	ATP
Memory address	0x07febc9889f70
Formula	C10H12N5O13P3
Compartment	c
Initial Condition	None
In 0 reaction(s)

Note that converted metabolites do not retain any references to the previously associated cobra.Reaction or cobra.Model.

for metabolite in [cobra_metabolite, mass_metabolite]:
    print("Number of Reactions: {0}; Model: {1}".format(len(metabolite.reactions), metabolite.model))

Number of Reactions: 13; Model: e_coli_core
Number of Reactions: 0; Model: None

However, all attributes that the mass.MassMetabolite object inherits from the cobra.Metabolite object are preserved:

for attr in ["id", "name", "formula", "charge", "compartment"]:
    print("Identical '{0}': {1}".format(
        attr, getattr(cobra_metabolite, attr) == getattr(mass_metabolite, attr)))

Identical 'id': True
Identical 'name': True
Identical 'formula': True
Identical 'charge': True
Identical 'compartment': True

12.1.2. Reaction to MassReaction

To convert a cobra.Reaction to a mass.MassReaction:

mass_reaction = MassReaction(cobra_reaction)
mass_reaction

Reaction identifier	PGI
Name	glucose-6-phosphate isomerase
Memory address	0x07febc9889fd0
Subsystem
Kinetic Reversibility	True
Stoichiometry	g6p_c <=> f6p_c D-Glucose 6-phosphate <=> D-Fructose 6-phosphate
GPR	b4025
Bounds	(-1000.0, 1000.0)

Upon conversion of a reaction, all associated cobra.Metabolite objects are converted to mass.MassMetabolite objects.

for metabolite in mass_reaction.metabolites:
    print(metabolite, type(metabolite))

g6p_c <class 'mass.core.mass_metabolite.MassMetabolite'>
f6p_c <class 'mass.core.mass_metabolite.MassMetabolite'>

If there are genes present, they are copied from one reaction to another in order to create a new cobra.Gene object for the MassReaction.

print(cobra_reaction.genes)
print(mass_reaction.genes)

frozenset({<Gene b4025 at 0x7febb95c2b80>})
frozenset({<Gene b4025 at 0x7feba0e93eb0>})

All other references to COBRApy objects are removed.

print(cobra_reaction.model)
print(mass_reaction.model)

e_coli_core
None

All attributes that the mass.MassReaction object inherits from the cobra.Reaction object are preserved upon conversion.

for attr in ["id", "name", "subsystem", "bounds", "compartments", "gene_reaction_rule"]:
    print("Identical '{0}': {1}".format(
        attr, getattr(cobra_reaction, attr) == getattr(mass_reaction, attr)))

Identical 'id': True
Identical 'name': True
Identical 'subsystem': True
Identical 'bounds': True
Identical 'compartments': True
Identical 'gene_reaction_rule': True

12.1.3. Model to MassModel

To convert a cobra.Model to a mass.MassModel:

mass_model = MassModel(cobra_model)
mass_model

Name	e_coli_core
Memory address	0x07feba0e93f70
Stoichiometric Matrix	72x95
Matrix Rank	67
Number of metabolites	72
Initial conditions defined	0/72
Number of reactions	95
Number of genes	137
Number of enzyme modules	0
Number of groups	0
Objective expression	1.0*Biomass_Ecoli_core - 1.0*Biomass_Ecoli_core_reverse_2cdba
Compartments	cytosol, extracellular

During conversion, the original cobra.Model remains untouched, while a new mass.MassModel is created using the equivalent mass objects. All references to the original cobra.Model are updated with references to the newly created mass.MassModel.

print("All MassMetabolites: {0}".format(
    all([isinstance(met, MassMetabolite)
         for met in mass_model.metabolites])))
print("All MassReactions: {0}".format(
    all([isinstance(rxn, MassReaction)
         for rxn in mass_model.reactions])))

All MassMetabolites: True
All MassReactions: True

12.2. Differences between COBRA and MASS

Although there are several similarities between COBRApy and MASSpy, there are some key differences in behavior that are worth highlighting.

12.2.1. COBRA vs. MASS reactions

There are some key differences between cobra.Reaction and mass.MassReaction objects. They are summarized below:

12.2.1.1. reversible vs. reversibility attributes

One key difference observed is how a reaction direction is determined. A cobra.Reaction utilizes the lower and upper bound values to determine the reversibility attribute.

print(cobra_reaction.reaction)
print(cobra_reaction.bounds)
print(cobra_reaction.reversibility)

g6p_c <=> f6p_c
(-1000.0, 1000.0)
True

Changing the reaction bounds affects the direction a reaction can proceed:

for header, bounds in zip(["Both Directions", "Forward Direction", "Reverse Direction"],
                          [(-1000, 1000), (0, 1000), (-1000, 0)]):

    print("\n".join((header, "-" * len(header))))
    cobra_reaction.bounds = bounds
    print(cobra_reaction.reaction)
    print(cobra_reaction.bounds)
    print("Reversibility: {0}\n".format(cobra_reaction.reversibility))

Both Directions
---------------
g6p_c <=> f6p_c
(-1000, 1000)
Reversibility: True

Forward Direction
-----------------
g6p_c --> f6p_c
(0, 1000)
Reversibility: False

Reverse Direction
-----------------
g6p_c <-- f6p_c
(-1000, 0)
Reversibility: False

Although MassReaction objects still have the reversibility attribute based on reaction bounds, the reaction rate equation is based on the reversible attribute. Additionally, the displayed reaction arrow for a reaction string now depends on the reversible attribute, rather than the reversibility attribute.

Therefore, even if the flux is constrained to proceed in one direction by the bounds, the kinetic rate expression still accounts for a reverse rate.

for header, bounds in zip(["Forward Direction (Flux)", "Reverse Direction (Flux)", "Both Directions (Flux)"],
                          [(0, 1000), (-1000, 0), (-1000, 1000)]):

    print("\n".join((header, "-" * len(header))))
    mass_reaction.bounds = bounds
    print(mass_reaction.reaction)
    print(mass_reaction.bounds)
    print("Reversibility: {0}".format(mass_reaction.reversibility))
    print("Reversible (Kinetic): {0}".format(mass_reaction.reversible))
    print("Rate: {0}\n".format(mass_reaction.rate))

Forward Direction (Flux)
------------------------
g6p_c <=> f6p_c
(0, 1000)
Reversibility: False
Reversible (Kinetic): True
Rate: kf_PGI*(g6p_c(t) - f6p_c(t)/Keq_PGI)

Reverse Direction (Flux)
------------------------
g6p_c <=> f6p_c
(-1000, 0)
Reversibility: False
Reversible (Kinetic): True
Rate: kf_PGI*(g6p_c(t) - f6p_c(t)/Keq_PGI)

Both Directions (Flux)
----------------------
g6p_c <=> f6p_c
(-1000, 1000)
Reversibility: True
Reversible (Kinetic): True
Rate: kf_PGI*(g6p_c(t) - f6p_c(t)/Keq_PGI)

Changing the reversible attribute affects the kinetic rate expression for the reaction, but it does not affect the reaction bounds.

for header, reversible in zip(["Both Directions (Kinetics)", "Forward Direction (Kinetics)"], [True, False]):
    print("\n".join((header, "-" * len(header))))
    mass_reaction.reversible = reversible
    print(mass_reaction.reaction)
    print(mass_reaction.bounds)
    print("Reversibility: {0}".format(mass_reaction.reversibility))
    print("Reversible (Kinetic): {0}".format(mass_reaction.reversible))
    print("Rate: {0}\n".format(mass_reaction.rate))

Both Directions (Kinetics)
--------------------------
g6p_c <=> f6p_c
(-1000, 1000)
Reversibility: True
Reversible (Kinetic): True
Rate: kf_PGI*(g6p_c(t) - f6p_c(t)/Keq_PGI)

Forward Direction (Kinetics)
----------------------------
g6p_c --> f6p_c
(-1000, 1000)
Reversibility: True
Reversible (Kinetic): False
Rate: kf_PGI*g6p_c(t)

To obtain the reaction in the reverse direction instead of the forward direction, the MassReaction.reverse_stoichiometry() method can be used. Setting inplace=False produces a new reaction, while setting inplace=True modifies the existing reaction. Setting reverse_bounds=True switches the lower and upper bound values.

mass_reaction_rev = mass_reaction.reverse_stoichiometry(inplace=False)
print(mass_reaction_rev.reaction)
print(mass_reaction_rev.bounds)
print("Reversibility: {0}".format(mass_reaction_rev.reversibility))
print("Reversible (Kinetic): {0}".format(mass_reaction_rev.reversible))
print("Rate: {0}\n".format(mass_reaction_rev.rate))

f6p_c --> g6p_c
(-1000, 1000)
Reversibility: True
Reversible (Kinetic): False
Rate: kf_PGI*f6p_c(t)

12.2.1.2. flux vs. steady_state_flux attributes

Another difference observed between cobra.Reaction and mass.MassReaction is how flux values are stored. When a model is optimized for FBA, the flux attribute of the reaction reflects the solution directly produced by the solver.

cobra_model = cobra.test.create_test_model("textbook")
cobra_model.optimize()
cobra_reaction = cobra_model.reactions.get_by_id("PGI")
cobra_reaction.flux

4.860861146496812

MassModel objects retain their ability to be optimized for FBA. Consequently, the ability to retrieve a solution for a reaction using the flux attribute is also retained.

cobra_model = cobra.test.create_test_model("textbook")
mass_model = MassModel(cobra_model)
mass_model.optimize()
mass_reaction = mass_model.reactions.get_by_id("PGI")
mass_reaction.flux

4.860861146496812

The value of the flux attribute is not the same as the steady_state_flux attribute, which is used in various mass methods:

print(mass_reaction.steady_state_flux)

None

To set steady_state_flux attributes for all reactions based on the optimization solutions, the MassModel.set_steady_state_fluxes_from_solver() method is used.

mass_model.set_steady_state_fluxes_from_solver()
# Display for first 10 reactions
for reaction in mass_model.reactions[:10]:
    print(reaction.id, reaction.steady_state_flux)

ACALD 0.0
ACALDt 0.0
ACKr 0.0
ACONTa 6.007249575350327
ACONTb 6.007249575350327
ACt2r 0.0
ADK1 0.0
AKGDH 5.064375661482088
AKGt2r 0.0
ALCD2x 0.0