Search code examples
pythonsbmlcbmpy

Add new reactions with GPR

I have a model and want to add an entire new pathway to it. Some of the metabolites already exist in the model, others have to be created. I also have to add GPRs to the reactions using genes not yet present in the model.

I found the function addReaction, but always get an error when I use it:

import cbmpy

cmod = cbmpy.CBRead.readSBML3FBC('model.xml')

cmod.addReaction('R_foo')

AssertionError: ERROR: requires a Reaction object, not something of type <type 'str'>

Any ideas how I can pass a reaction object and add metabolites and a GPR?

Solution

  • You are looking for createReaction. The following will work (I use the model from this question):

    import cbmpy as cbm

mod = cbm.CBRead.readSBML3FBC('e_coli_core.xml')

mod.createReaction('R_foo')

    This will print

    Reaction "R_foo" bounds set to: -INF <= R_foo <= INF Add reagents with cmod.createReactionReagent(R_foo, metabolite, coefficient)

    So, by default, one adds reversible reactions (see below how to add an irreversible one) and it also tells you how to add reagents.

    Let's first assume that you add a reactions for which all reagents are already present in the model. Then you can use createReactionReagent to add reagents along with their stoichiometric factor like this:

    mod.createReactionReagent('R_foo', 'M_fum_c', -1.)
mod.createReactionReagent('R_foo', 'M_nh4_c', 5.)

    We can check whether the reaction was added correctly:

    mod.getReaction('R_foo').getStoichiometry()

    will return

    [(-1.0, 'M_fum_c'), (5.0, 'M_nh4_c')]

    You can then easily add a GPR to the reaction using createGeneProteinAssociation:

    mod.createGeneProteinAssociation('R_foo', 'gene_1 or gene_2')

    Again the check whether it worked as intended:

    mod.getGPRforReaction('R_foo').getAssociationStr()

    yields:

    '((gene_1 or gene_2))'

    If the genes are not present in the model, they will be added automatically:

    mod.getGeneIds()[-2:]

    will return

    ['gene_1', 'gene_2']

    As you want to add an entire pathway, we do the same now for a second reaction with a reagent not yet part of the model:

    # add an irreversible reaction
mod.createReaction('R_bar', reversible=False)
mod.createReactionReagent('R_bar', 'M_succ_c', -1.5)

    Let's assume, the metabolite you want to add is called A, then

    mod.createReactionReagent('R_bar', 'A', 1.0)

    will fail with

    AssertionError: Metabolite A does not exist

    That means that we first have to create it using createSpecies:

    mod.createSpecies('A', name='species A', compartment='c', charge=-2, chemFormula='C1H2O3')

    The arguments name till chemFormula are not required. Now you can call 

    mod.createReactionReagent('R_bar', 'A', 1.0)

    You can check this question on how to add additional annotation to the species or reactions.

    You might then also want to add all reactions belonging to this pathway into a group which makes accessing the reactions very easy later on:

    pw_reactions = ['R_foo', 'R_bar']

# create an empty group
mod.createGroup('my_new_pathway')

# we can only add objects to a group so we get the reaction object for each reaction in the pathway
reaction_objects = [mod.getReaction(ri) for ri in pw_reactions]

# add all the reaction objects to the group
new_pw.addMember(reaction_objects)

    Now you can access the members of this group using

    mod.getGroup('my_new_pathway').getMemberIDs()
# returns ['R_foo', 'R_bar']

    if you are interested in the reaction IDs or 

    mod.getGroup('my_new_pathway').getMembers()

    if you are interested in the reaction objects themselves.