Add new reactions with GPR

Question

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?

Answer 1

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.