Adding reactions from another model does not reset the reactions their metabolites refer to

[Midnighter]

From the mailing list:

Thank you! That is very helpful. If I add a reaction, the metabolites are also added (good!) but all reactions from the universal model associated with those metabolites are included in the metabolite description, even though they don't exist in the model. Is this the expected behavior?

For example, here I've only added EX_fru_e, none of the other reactions are in 'model'.

In [310]: model.add_reactions([universe.reactions.EX_fru_e])
In [311]: model.metabolites.fru_e.reactions
Out[311]:
frozenset({<Reaction FRUt2r at 0x7fa918d40630>,
           <Reaction KESTOTTRASEe at 0x7fa90f107048>,
           <Reaction KESTOPTASEe at 0x7fa90f107320>,
           <Reaction LEVANASE_2e at 0x7fa90f124550>,
           <Reaction INULINASEe at 0x7fa90f161d68>,
           <Reaction OLIGOFRUASEe at 0x7fa90f09c550>,
           <Reaction LEVANASE_3e at 0x7fa90f124588>,
           <Reaction EX_fru_e at 0x7fa911398390>,
           <Reaction LEVANASE_4e at 0x7fa90f1245c0>,
           <Reaction KESTOASEe at 0x7fa90f1071d0>})

[cdiener]

Yeah that is not handled correctly at all. I would make metabolite.reactions a dynamic property that just gets its reactions from the model with query. Another problem with that code example is that the universal model is broken afterwards (the reaction is stripped from it and there are dangling solver variables). This is part of the "can a reaction belong to two models" discussion that was never concluded.

[pstjohn]

Yeah that's an unfortunately leftover from cobrapy's early days. Probably the only reason we still have a Species class.

We should probably have add_reaction and add_metabolite check for an exisiting model instance, and if one is found, copy the incoming object.

All this code will have to get fixed:

$ git grep '\._reaction'
cobra/core/gene.py:        for the_reaction in list(self._reaction):
cobra/core/gene.py:        self._reaction.clear()
cobra/core/model.py:                new_met._reaction.add(new_reaction)
cobra/core/model.py:                new_gene._reaction.add(new_reaction)
cobra/core/model.py:                for the_reaction in list(x._reaction):
cobra/core/model.py:                for x in list(x._reaction):
cobra/core/model.py:                    model_metabolite._reaction.add(reaction)
cobra/core/model.py:                            model_metabolite._reaction.remove, reaction))
cobra/core/model.py:                    if reaction in met._reaction:
cobra/core/model.py:                        met._reaction.remove(reaction)
cobra/core/model.py:                            context(partial(met._reaction.add, reaction))
cobra/core/model.py:                        if remove_orphans and len(met._reaction) == 0:
cobra/core/model.py:                    if reaction in gene._reaction:
cobra/core/model.py:                        gene._reaction.remove(reaction)
cobra/core/model.py:                            context(partial(gene._reaction.add, reaction))
cobra/core/model.py:                        if remove_orphans and len(gene._reaction) == 0:
cobra/core/model.py:                met._reaction.clear()
cobra/core/model.py:                gene._reaction.clear()
cobra/core/model.py:                    met._reaction.add(rxn)
cobra/core/model.py:                    gene._reaction.add(rxn)
cobra/core/reaction.py:            g._reaction.add(self)
cobra/core/reaction.py:                    g._reaction.remove(self)
cobra/core/reaction.py:            x._reaction.add(self)
cobra/core/reaction.py:            x._reaction.add(self)
cobra/core/reaction.py:            x._reaction.add(self)
cobra/core/reaction.py:            x._reaction.add(self)
cobra/core/reaction.py:                metabolite._reaction.add(self)
cobra/core/reaction.py:                metabolite._reaction.remove(self)
cobra/core/reaction.py:        cobra_gene._reaction.add(self)
cobra/core/reaction.py:        cobra_gene._reaction.discard(self)
cobra/core/species.py:        self._reaction = set()
cobra/core/species.py:        return frozenset(self._reaction)
cobra/io/sbml.py:        if len(metabolite._reaction) == 1:
cobra/io/sbml.py:            if list(metabolite._reaction)[0].id.startswith("EX_"):
cobra/manipulation/delete.py:        if len(the_metabolite._reaction) == 0:
cobra/manipulation/delete.py:        potential_reactions.update(gene._reaction)
cobra/manipulation/modify.py:            recompute_reactions.update(old_gene._reaction)
cobra/manipulation/modify.py:                gene._reaction.add(reverse_reaction)
cobra/test/test_model.py:        assert dummy_reaction in actual_metabolite._reaction
cobra/test/test_model.py:        m2.genes.b4025._reaction.clear()

[Midnighter]

Yes, probably during the active time of GSoC, we will prioritize planning and implementing cobrapy 1.0. It's time to lose some baggage.

[ensakz]

Is there any progress on fixing this bug? or do you have any workaround suggestions for solving it?

Thanks in advance

[cdiener]

It's only Metabolite.reactions that is buggy. You can always loop through reactions yourself:

metab = ...
rxns = [rxn for rxn in model.reactions if metab in rxn.metabolites]

This will only give you reactions from the model.

[ensakz]

In my case, after I save the model using cobra.io.save_json_model() function, the reactions related to the metabolite of interest gets updated and after I load json model file again I can see reactions which are from my model of interest (not from another model) when I use model.metabolite.some_metabolite.reactions function. I use cobrapy 0.17.0 version though

[akaviaLab]

I think that this section in model.py:add_reactions()

            for gene in list(reaction._genes):
                # If the gene is not in the model, add it
                if not self.genes.has_id(gene.id):
                    self.genes += [gene]
                    gene._model = self

                    if context:
                        # Remove the gene later
                        context(partial(self.genes.__isub__, [gene]))
                        context(partial(setattr, gene, "_model", None))

                # Otherwise, make the gene point to the one in the model
                else:
                    model_gene = self.genes.get_by_id(gene.id)
                    if model_gene is not gene:
                        reaction._dissociate_gene(gene)
                        reaction._associate_gene(model_gene)

Can be replaced by update_genes_from_gpr()

Am I right?

[cdiener]

@akaviaLab I think you would know better than us at that point. Feel free to send a PR.