Runtime error when attempting to remove genes from model

[djinnome]

Hi folks,

If you do the following in the Docker image:

import pickle, cobrame
with open('/home/meuser/me_models/iJL1678b.pickle', 'rb') as f:
     me = pickle.load(f)
     me.remove_genes_from_model( ['b0002'] )

Then I get the following error:

Complex (ASPKINIHOMOSERDEHYDROGI-CPLX) removed from model
Traceback (most recent call last):
  File "/source/cobrame/cobrame/core/model.py", line 385, in remove_genes_from_model
    from optlang.interface import SymbolicParameter
ImportError: cannot import name 'SymbolicParameter'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/source/cobrame/cobrame/core/model.py", line 387, in remove_genes_from_model
    protein.remove_from_model(method='destructive')
  File "/usr/lib/python3.6/site-packages/cobra-0.5.11-py3.6-linux-x86_64.egg/cobra/core/Metabolite.py", line 142, in remove_from_model
    for x in self._reaction:
RuntimeError: Set changed size during iteration

Please advise

[djinnome]

A bit more context:

I would like to comprehensively knock out genes in the ME model, so if there is a workaround to perform this operation, that would be great.

CobraPy has remove_from_model(), which is actually called during the CobraME remove_genes_from_model() method, but a reasonable workaround for me would be just a CobraME version of the CobraPy gene.knock_out()

The problem is that the pickled ME object does not contain any genes

>>> me.genes
[]

Perhaps something as simple as this:

def knock_out_genes_from_model( self, gene_list ): 
   for gene in gene_list:
            # Find all complexes the gene product is part of and knock out the associated reactions
            protein = self.metabolites.get_by_id('protein_'+gene)
            for cplx in protein.complexes:
                print('Complex (%s) knocked out in model' % cplx.id)
                for rxn in cplx.metabolic_reactions:
                    rxn.bounds = (0,0)

[coltonlloyd]

Hi Jeremy,

The issue here is that there is a bug in COBRApy 0.5.11 that breaks the remove_from_model(method='destructive) call. I've known about this but forgot to address it so thanks for the reminder!

While I'm correcting this, you could avoid this problem by just changing the bounds of each gene's translation reaction (if the gene is translated into a protein) or, more generally, by changing the bounds of each reaction that the gene participates as a reactant in (shown in the example below). You can do this more effectively using a function like below (which was modified from https://github.com/coltonlloyd/cobrame_supplement/blob/master/Table_3_model_essentiality/essentiality.py):

import json
import ecolime
from cobrame.core.reaction import TranscriptionReaction
from qminospy.me1 import ME_NLP1

def compute_gene_essentiality_at_growth_rate(me, gr, out_location):
    me_nlp = ME_NLP1(me, growth_key='mu')
    expressions = me_nlp.compile_expressions()
    me_nlp.compiled_expressions = expressions

    hs = None

    all_genes = me.metabolites.query(re.compile("^RNA_b[0-9]"))

    results = {}
    for gene_RNA in list(all_genes):

        default_bounds = {}
        for r in gene_RNA.reactions:
            if not r.id.startswith("DM") and not \
                    isinstance(r, TranscriptionReaction):
                default_bounds[r] = (r.lower_bound, r.upper_bound)
                r.knock_out()
        x, status, hs = me_nlp.solvelp(gr, basis=hs)

        if status == 'optimal':
            results[gene_RNA.id] = 'NONESSENTIAL'
        else:
            results[gene_RNA.id] = 'ESSENTIAL'

        print("%s\t%s" % (gene_RNA.id.split("_")[1], str(status)))

        with open("%s/iJL1678b_essentiality_%.2f_gr.json" % (out_location, gr),
                  "w") as outfile:
            json.dump(results, outfile, indent=True)

        # Reset bounds
        for r in default_bounds:
            r.lower_bound = default_bounds[r][0]
            r.upper_bound = default_bounds[r][1]

If you use a small enough growth rate (I think I used .1 in manuscript), this will save you from having to run a full bisection for each gene. This also compiles all of the symbolic expressions upfront and reuses the previous simulations basis vector to speed up computations.

Hope this helps!

[djinnome]

Thanks for the quick response and the detailed solution. I will definitely try it out and let you know how it works.!