Skip to content

[BUG] Gene equality outside of models #1434

Description

@oxinabox

Hi,
Several bugs around the handling of genes in reactions when not connected to models
which were found by @jedforrest.
Mostly they boil down to genes with same ID not always being treated as identical.
I think they probably should be, since e.g. gapfilling treats reactions of same ID as identical.
And once they are in a model then the same ID is treated as identical, for perventing duplicates being added when adding a series of reactions.

Declared seperately in GPR from string not equal

import ast
from cobra.core.gene import GPR
from cobra import Model, Reaction

def print_genes(genes):
    print(sorted(genes, key=lambda x: x.id))

# updating GPR by string
rxn1 = Reaction('rxn1')
rxn1.gene_reaction_rule = "(A and B) or C"
genes_abc = rxn1.genes
rxn1.gene_reaction_rule = "A and B"
genes_ab = rxn1.genes

Results in:

In [2]: genes_ab.issubset(genes_abc) # False!
Out[2]: False

In [3]: print_genes(genes_abc)
   ...: print_genes(genes_ab)
[<Gene A at 0x7e3f246f2a80>, <Gene B at 0x7e3f15f2bce0>, <Gene C at 0x7e3f16a8d460>]
[<Gene A at 0x7e3f98e1fe60>, <Gene B at 0x7e3f99cf0fb0>]

Created via GPR manipulation:

This one might arguably need model.repair called to fix it.
But this isa Reaction so there is no such model.

# updating GPR by AST Expression
rxn2 = Reaction('rxn2')
rxn2.gene_reaction_rule = "(A and B) or C"
genes_abc = rxn2.genes
sub_op = rxn2.gpr.body.values[0]  # '(A and B)'
rxn2.gpr = GPR(ast.Expression(sub_op))
genes_ab = rxn2.genes

similarly gives

In [5]: genes_ab.issubset(genes_abc) # False!
Out[5]: False

In [6]: print_genes(genes_abc)
   ...: print_genes(genes_ab)
[<Gene A at 0x7e3f1b5c33e0>, <Gene B at 0x7e3f1682f620>, <Gene C at 0x7e3f1930d310>]
[<Gene A at 0x7e3f16a8d460>, <Gene B at 0x7e3f15f2bce0>]

adding a reaction to a creates new nonequal genes

model = Model("model")
rxn3 = Reaction('rxn3')
rxn3.gene_reaction_rule = "(A and B) or C"
rxn3_genes_orig = list(rxn3.genes)

model.add_reactions([rxn3])
In [8]: 
   ...: rxn3.genes == rxn3_genes_orig  # False!
Out[8]: False

In [9]: print_genes(rxn3_genes_orig)
   ...: print_genes(rxn3.genes)
[<Gene A at 0x7e3f1682e5a0>, <Gene B at 0x7e3f998235c0>, <Gene C at 0x7e3f1a7304a0>]
[<Gene A at 0x7e3f98e485f0>, <Gene B at 0x7e3f98e482f0>, <Gene C at 0x7e3f98e481d0>]

I have seen a situation, though I don't have a MWE where I ended up with duplicate genes in the same model.Genes with the same id.
Which breaks the DictList.
But I don't have a small reproducer yet.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions