GPR represntation in pyCobra

Anonymous
2011-09-23
2013-05-30

  • Anonymous
    2011-09-23

    I have a question regarding how you treat (and plan to treat) GPRs in pyCobra.  I am writing a script where I need to use the delete_model_genes function to restrict flux through reactions as well as delete the genes.  I want to actually remove all traces of the gene from the model.

    I am using the command: gene.remove_from_model(model) to do this, however it doesn't seem to be working (looks like the remove_from_model code is inherited from the metabolite class which tries to remove my gene from the model's metabolites…)

    Also, I noticed your comments on switching to binary representations of GPRs in the future.  How are you planning to do this?  I want to be able to write/plan my code with that in mind.

    Thanks!

     
  • Daniel Hyduke
    Daniel Hyduke
    2011-09-24

    I hadn't written the remove_from_model function for cobra.Gene b/c I hadn't made a decision on what to do with the gene association representation and parsing those strings can be annoying.  I've updated the Gene class to have a remove_from_model function.

    You'll note that cobra.Gene has a functional attribute which is starting to come into play when mutating a gene into an inactive state; currently, only is boolean but that will possible change to a continuous representation.

    P.S. I'll probably change the word delete to mutate and then allow the user to specify the type of mutation: complete, inactivation, decreased efficiency, …. But that'll be a halloween project.

    I've updated the code, here's what happens with the salmonella.pickle in test/data

               
            the_reaction = cobra_model.reactions
            print the_reaction
            print the_reaction.gene_reaction_rule
            print repr(the_reaction._genes)
            the_gene = the_reaction._genes.keys()
            print the_gene
            the_gene.remove_from_model(cobra_model)
            print the_reaction.gene_reaction_rule
            print repr(the_reaction._genes)
            print 'lower bound: %f'%the_reaction.lower_bound
            print 'upper bound: %f'%the_reaction.upper_bound
            the_gene = the_reaction._genes.keys()
            print the_gene
            the_gene.remove_from_model(cobra_model)
            print the_reaction.gene_reaction_rule
            print repr(the_reaction._genes)
            print 'lower bound: %f'%the_reaction.lower_bound
            print 'upper bound: %f'%the_reaction.upper_bound

    3OAS140
    ( STM2378  or  STM1197 )
    {'STM1197': 1.0, 'STM2378': 1.0}
    STM1197
    ( STM2378  or  False )
    {'STM2378': 1.0}
    lower bound: 0.000000
    upper bound: 1000.000000
    STM2378
    ( False  or  False )
    {}
    lower bound: 0.000000
    upper bound: 0.000000

     

  • Anonymous
    2011-09-25

    Perfect! 

    Is there code to re-parse the gene_reaction_rule to remove the "FALSE" text.  Essentially I am looking for a way to update the gene_reaction_rule based on genes left in the model (I will be renaming them too, so would want the updated gene.name to be reflected in the gene_reaction_rule text). 

    I suppose this would be the opposite of what reaction.parse_gene_association() does, in effect using the Model.genes to update the reactions.parse_gene_association rule

     
  • Daniel Hyduke
    Daniel Hyduke
    2011-09-25

    At the moment, you should do this yourself as it brings up the issue of what to do for an enzyme activity that requires the gene that you've just removed from the model.  In reality, the loss of an essential subunit from a heteromer would render the activity null; however, if you're dealing with an organism that has that activity then a different locus may have taken over the role.

    If you want to remove the False statement, look at the Gene.remove_from_model code and manipulation.delete.delete_model_genes code for starters on parsing the string.  If you want to construct a gene association that uses Gene.name instead of Gene.id then look at the Gene.remove_from_model (see below).  Just be warned that functions that interact with Reaction.gene_reaction_rule expect Gene.id not Gene.name in the text string.
    the_reaction = cobra_model.reactions
    the_gene_reaction_relation = the_reaction.gene_reaction_rule
    for other_gene in the_reaction._genes:
        other_gene_re = re.compile('(^|(?<=( |\()))%s(?=( |\)|$))'%re.escape(other_gene.id))
        the_gene_reaction_relation = other_gene_re.sub(other_gene.name, the_gene_reaction_relation)