GPR to Sympy

[akaviaLab]

This is a follow-up of the discussion on #1013

To summarize

1. I will use sympy, but make it possible to use optlang when some bugs in symengine are fixed.
2. Gene will inherit from both Species and sympy.Symbol

One problem I've run into is the fact that sympy uses name to identify Symbols (like cobrapy uses id).
I can't stop sympy or optlang from using name without rewriting most of Sympy, which is not an option. Sympy will not recognize id, since many functions look to *.name
Currently, the Gene.name is the longer string that is potentially human readable, but with the sympy changes Gene.name will be equivalent to id.

I'm thinking of adding a field called "alias" to Gene, which will have the role of the current name (a longer/different str). I will then refactor all the internal functions that use Gene.name to use Gene.alias.

Would that work fine?
If it does, should Gene be created with name as a parameter, even if it will be internally renamed alias?
Something like this to keep consistency with previous usages of Gene
class Gene(Symbol, Species):
def init(self, id=None, name="", functional=True):
super().init(name=id)
super(Species, self)
self.alias = name
self.id = self.name

Or should it become def init(self, name=None, alias="", functional=True)

What about functions like (in Reaction)
@Property
def gene_name_reaction_rule(self):
"""Display gene_reaction_rule with names intead.

    Do NOT use this string for computation. It is intended to give a
    representation of the rule using more familiar gene names instead of
    the often cryptic ids.

I think one option is to keep it, have it call gene_alias_reaction(self) and maybe print out a warning it will be deprecated. Or I can just keep it as is, and refer internally to alias.

If not, do you guys have suggestions for what would work?

One additional topic to think about is refactoring all of cobrapy to use alias (or whatever we decide). I think that might be overkill, but I'm opening to hearing what you say.

Thank you

Uri David

[matthiaskoenig]

Can't you just use a special object which is used for the GPR calculation which is based on a Gene?
Something along the line of

class GPRGene:
    def __init__(self, gene: Gene):
        self.name = gene.id
        self.gene = gene

You would use this for all the GPR parsing. I don't think it is a good idea to use aliases. People expect the name to be a certain thing. To change the complete code base because some library expects something is a bad idea. You have to create an abstraction layer as suggested above.

[akaviaLab]

I'm a bit confused.
This abstraction looks like it undercuts the idea of Gene being a Sympy/Symbol capable object.
If you want me to rewrite GPRs, so that GPRs are SYMPY expressions and Genes stay as they are, that's pretty similar to what I did before.
I could modify it that Genes stay as they are and GPRs are parsed as SYMPY (when reading), kept as SYMPY/SYMBOL and printed as str. This might use GPRGene, but I'm not sure this abstraction class would even be necessary in that case.

Am I missing/misunderstanding something here?

I could set up GPR to inherit from Gene and Symbol (or Sympy Expression), but in that case, what extra features would you like the GPR to have?

[akaviaLab]

Specifically, I'm confused by the comment made here that
#1013 (comment)
"Sorry, that was a typo. Was supposed to read "genes be symbols/entities". This way you can get them from a sympy expression using expression.atoms. Basically the Gene class would now derive from optlang.symbolics.Symbol."

I can set up GPRExpression to have a field called genes, which will use expression.atoms and a bit of processing to get you the Genes (which are NOT Sympy objects).

[akaviaLab]

One more thing. It looks like symengine has logical statements, but I can't figure out how to implement it in optlang. I've tried looking at the code and tutorials, but I'm not good with C++.

Could someone look at this symengine/symengine#154 and add the relevant code to optlang, please? I'll then write whatever extensions to use optlang.Symbol and optlang.Or and so forth.

Thank you

[cdiener]

Hmmm, that is unfortunate but not much we can do. What Matthias suggested is in line with what the rest of cobrapy is doing in with optlang variables already (which are sympy/symengine symbols as well). So that makes sense. I would call it GeneSymbol but apart from that that seems like the way to go. Each gene would have its own symbol at gene.symbol and GPRs are now always expressions of gene symbols and never strings. The reaction.functional system would also be switched to evaluate the expression.

Optlang symbols can already use boolean expressions, but sympify is still not working and this would be one of the parsers so that would be an issue.

[akaviaLab]

Okay, I think I got it.
GeneSymbol or GeneGPR, not sure which one is better. I have it GeneGPR for now, and I can change it later.
I would probably also need to create something like a GeneGPRExpression which would be derived from sympy.Expression.

When would sympify be used?
My code constructs an AST tree, and then turns it into a Sympy expression. This is complicated, but unfortunately RECON3D (and possibly others) have very long gprs that break sympify. My solution is an AST tree, which is processed. It is as fast as sympify, but it works for long gprs.

I think the same logic could be applied to symengine when visiting the AST tree, just calling optlang.or instead of sympy.or. I could use optlang if optlang were to call logical_or, logical_and. It looks like symengine refers to sympy.core.S, so either using sympy or symengine, I can keep using that. I'm just not sure how to modify optlang.

I would have to rewrite the reaction functional system, which is doable.
My feeling is that in this case, Gene.remove_from_model() should be removed and replaced with a simple call to cobra.manipulation.remove_genes instead. That way I/we don't have to maintain two versions.

Does this make sense?

[cdiener]

Found some time to run a simple benchmark and it looks pretty bad 😢

Python 3.9.8 (main, Nov  7 2021, 15:47:09) 
Type 'copyright', 'credits' or 'license' for more information
IPython 7.29.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: from sympy import Symbol, And

In [2]: syms = [Symbol("s%d" % i) for i in range(100)]

In [3]: sympy_ex = And(*syms)

In [4]: from cobra.core.gene import parse_gpr, eval_gpr

In [6]: str_ex = " and ".join(["s%d" % i for i in range(100)])

In [7]: ast, genes = parse_gpr(str_ex)

In [8]: subs = {s: True for s in syms}

In [9]: str_subs = {s: True for s in genes}

In [10]: # Build

In [11]: %timeit And(*syms)
4.02 ms ± 364 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [12]: %timeit ast, genes = parse_gpr(str_ex)
391 µs ± 20.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [13]: # Eval

In [14]: %timeit res = sympy_ex.xreplace(subs)
1.33 ms ± 25 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [15]: %timeit res = eval_gpr(ast, str_subs)
1.24 µs ± 42.5 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

Basically building the expression (even without parsing) is 10x slower with sympy and evaluating is about 1000x (!) slower than what we have now. This would add quite some time to parsing and a massive chunk of time to gene deletions. So for now I'm not sure how much sense it makes to have this especially since there seem to be no plans to add it to the symengine Python interface. So this is probably not in a usable state now.

Maybe as an alternative we could start to save the AST as the GPR and not the string representation. ast2str is pretty fast (10x faster than parse_gpr), so it would be feasible for printing. This way we would have a GPR class that can be build from a string and can convert to other formats (str, SBML, etc).

[Midnighter]

I think we need to distinguish two cases:

1. In modern SBML, the GPR is already an expression tree in the XML.
2. In older models, where the GPR is written in the notes field, this is text form.

Some alternatives for building an expression where the last one is equivalent to what optlang.symbols.add does (and an order of magnitude faster than parse_gpr):

In [1]: from sympy import Symbol, And

In [2]: syms = [Symbol(f"s{i}") for i in range(100)]

In [3]: %timeit And(*syms)
4.22 ms ± 36.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [4]: %timeit And(*syms, evaluate=False)
4.22 ms ± 1.38 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [5]: %timeit And.fromiter(args=syms)
4.2 ms ± 2.77 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [6]: %timeit And._from_args(syms)
57.5 µs ± 154 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

I didn't manage to speed up the evaluation part yet. I would like to try expr.subs(subs, simultaneous=True) but I get an error.

[cdiener]

I think there are more cases like JSON where it is also a string. So we need efficient parsers for each. Right now 1 and 2 are actually the same as the SBML reader builds up a string rep as well right now. Building up the AST from scratch is equally fast (that could be mapped directly from SBML):

In [38]: %timeit ex = Expression(BoolOp(And(), [Name(id=x, ctx=Load()) for x in s]))
46.5 µs ± 1.47 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

The question is whether there is a _from_args equivalent for a random order of boolean operations. Though that difference was smaller than the evaluations. Efficient evaluation is probably harder to get. At least from a cursory look neither .subs() nor .xreplace() is fast in any setting.

Symengine does not support booleans in the Python interface. But using the normal Add and Integers gives about 400us for sympify (which is good) and about 150us for xreplace (which is still 100x slower).

[cdiener]

So maybe better to wrap the AST in a nice encapsulated class? This would allow things like GPR.from_sbml, gpr.remove(gene1), gpr.to_symbolic, and gpr & gene1 which could behave somewhat to what we actually wanted without sacrificing too much speed.

[akaviaLab]

Some comments

1. I don't understand why OR isn't in the symengine interface
   Looking at https://github.com/symengine/symengine.py/blob/4978e0be5084e09bba54b3a0793de6da33a79e36/symengine/lib/symengine_wrapper.pyx#L428
   It seems that it parses sympy.Or to symengine logical_or, so it would be able to do so
   Am I missing something?
2. Does the gene keep track of other genes?
   Gene has the field functional.
   If I have
   rxn1.gene_reaction_rule = 'A & B'
   and then set model.get_genes_by_id('A').functional=false

What would/should eval_gpr(rxn1, list()) give me?
True or False

In general, I think that sympy might be a bit slower, but if we only do the conversion once, it should not be very slow. Right now, we convert from str to AST and eval the AST multiple times. I'll time it when I'm ready.

[cdiener]

1. It isn't supported in the Python library yet. In particular you can't build logical expressions from scratch with symengine which is something we would need:

In [1]: from symengine import Symbol

In [2]: x, y = Symbol("x"), Symbol("y")

In [3]: x & y
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-3-c19a6a941ba4> in <module>
----> 1 x & y

TypeError: unsupported operand type(s) for &: 'symengine.lib.symengine_wrapper.Symbol' and 'symengine.lib.symengine_wrapper.Symbol'

In [4]: x | y
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-4-9ffff6900e54> in <module>
----> 1 x | y

TypeError: unsupported operand type(s) for |: 'symengine.lib.symengine_wrapper.Symbol' and 'symengine.lib.symengine_wrapper.Symbol'

In [5]: from symengine import sympify

In [6]: sympify("x & y")
[1]    1415864 segmentation fault  ipython

2. You wouldn't directly set functional but use gene.knock_out which sets functional and checks for all associated reactions if they need to be knocked out. That just requires the gpr evaluation to be fast.

In general, I think that sympy might be a bit slower, but if we only do the conversion once, it should not be very slow. Right now, we convert from str to AST and eval the AST multiple times. I'll time it when I'm ready.

This is what I tested in the benchmarks above. Even with a parsed expression, evaluation is about 3-4x slower than what we have now. The parsing is also slower than what we have now, so that is a best case estimate. That would make gene deletions at least 4x slower which would be problematic. As a reference, when we converted to the optlang interface we allowed a maximal slow-down of 20% for core operations (so 1.2x slower).
It is true that you could make it much faster by using the current parser and saving the AST tree instead of the string as the GPR expression. This is also benchmarked above. This would be a >100x speed-up compared to what is used now, but would not use a symbolic interface.

[Midnighter]

Thinking about it some more, one only really needs symbolic expressions when transforming them, for example, simplifying the expression, conversion to conjunctive or disjunctive normal form. It seems to me that these are rare cases. Simplifying an expression would only be helpful for storing it that way for quicker evaluation and I don't think we'd gain much. Otherwise the expressions have a clear biological meaning and should not be simplified. Conversion to disjunctive normal form has its uses but probably is very rarely desired.

The only remaining part is building up an expression. Dropping symbolic expressions and only using AST would mean an entire string has to be supplied for parsing. Is that an issue? I suspect that in the majority of cases models with complete GPR are read and only applied, thus storing a parsed AST ready for evaluation seems like a good way forward.

[cdiener]

Yeah this would be my pref as well and we could still add functions to convert to sympy if required (for CNF for instance). You can build an AST without a string like in my example above or maybe I misunderstood what you mean by "an entire string has to be supplied". So there is a lot of potential speed up in SBML parsing of GPRs and eval.

Disadvantages:

• can't trivially check whether a gene is in a GPR (but you can walk the AST)
• transformations like CNF, simplification not simple
• equality checking not simple

Advantages:

• wicked fast (10x faster parsing 100x faster eval)
• already have a fast parser
• no modifications to Gene are necessary

[akaviaLab]

Okay.
So the way I understand this discussion, and what I would like to do is as follows. Can you please tell me if I've got this right?

1. Keep using AST, but shift GPR to be an AST, where
   Reaction
   @Property
   def gene_reaction_rule(self):
   return self._gene_reaction_rule
   Becomes something like
   return ast2str(self._gene_reaction_rule)

Shift all other modules to use this property. I think delete compiles the string AST, which would not be necessary for example. SBML io parses AST to string and back again, which could probably be simplified.

2. Add Sympy as an additional function, which will allow CNF, DNF and logical comparisons of GPRs
   One reason I started using Sympy is that AST can't understand that 'A & B' is identical to 'B & A', let alone more complicated expressions.

3. Gene expression question
   Does cobrapy have something that would calculate the expression of a reaction by parsing the GPR and gene expression?
   rxn1.gene_rule_set = 'A & B'
   Expression(A) = 5
   Expression(B) = 8
   Expression(rxn1) = min(Expression(A), Expression(B)) = 5

If cobrapy does have this, where is it defined? I've looked at the code and haven't found it
If it doesn't have this option, would you like me to add this feature?

[cdiener]

1. Yep pretty much. Might be nice to make gene reaction rule a class since it gives a much nicer API. Something like:

import ast

# ...

class GPR(ast.AST):

    @staticmethod
    def from_string(string_gpr):
        # ...

    @staticmethod
    def from_sbml(sbml_exp):
        # ...
    
    def eval(self):
        # ...
    
    def __str__(self):
        # ...
    
    def __repr_html__(self):
        # ...

    def as_symbolic(self):
        # ...

2. Yep that would be great. This should be fairly efficient from a parsed AST. You basically just have to map the AST operations and names to sympy operators and symbols. The eval_gpr code or ast2str code is a good example how to do that.

3. The GPR does not use anything else than binary expression values. Enzyme abundance and fluxes are usually not correlated so most people use some discrete scheme for expression here. An example for expression confidence scores (like used by CORDA, GIMME and Fastcore) can be found here. A lot of that was supposed to be implemented in driven but there is nobody that wanted to take the lead on that.

[akaviaLab]

1. Thank you. working on that.
2. Will do
3. I might add the bit linked from CORDA as a sub-function of GPR class if someone wants it, and I'll leave everything else in driven for someone who needs it.

I've set it up so reaction.gene_reaction_rule is basically the string version of GPR, and GPR is based on AST. reaction.gene_reaction_rule is no longer an entity by itself
Some questions
If reaction_gene_rule is deriving from AST, and no longer a string, what should happen with malformed gene rules?
https://github.com/akaviaLab/cobrapy/blob/0bbad14ee36406ea14c06f44e9cdd216182b945a/src/cobra/test/test_core/test_core_reaction.py#L16
contains defective GPRs. If I feed them to my class, it complains, and the gene_reaction_rule remains empty string.
What do you think should happen in these cases? To the gene_reaction_rule? To the reaction.genes?

[cdiener]

reaction.gene_reaction_rule is no longer an entity by itself

You mean that Reaction._gene_reaction_rule now points to a GPR instance, right?

If reaction_gene_rule is deriving from AST, and no longer a string, what should happen with malformed gene rules?

I think it should throw an error with a clear statement that it is not a valid GPR. Since the GPR is now an object its default value should switch to None. Just need to keep this on our radar since it will require a minor version bump.

Feel free to split up PRs into smaller parts. For instance, by starting with a smaller implementation that only gets to feature parity with the current implementation and another one that adds the new features. This makes it easier to review.

[akaviaLab]

Reaction._gene_reaction_rule is effectively
Reaction._gpr.gene_reaction_rule

which is
str(Reaction._gpr), which is ast2str(Reaction._gpr)

I would have like to split the PR to parts, where the first one adds GPR to cobra.core.gene and doesn't change almost anything (although some test might fail) and the second part removes all irrelevant functions so ast2str, parse_gpr, eval_gpr become internal functions.

However, I already pushed some of the extra commits to Github. How do I fix that, so the pull request won't be ALL the commits?

[cdiener]

I don't think you need Reaction._gpr.gene_reaction_rule. The GPR is the gene_reaction_rule so it does not need this entity. Just having Reaction._gene_reaction_rule point to a GPR object is sufficient. Keeping around stringified versions just gets rid of all the speed gains again. So just having something like self._gene_reaction_rule = GPR() in the init code should be enough. Then the user can assign with rxn.gene_reaction_rule = "gene1 and gene2" to set it where the setter translates this to rxn._gene_reaction_rule = GRP.from_string("gene1 and gene2"). You can also overwrite GPR.__repr__ if you want it to look nice when writing rxn.gene_reaction_rule but I think it's enough to overwrite GPR.__str__ and then call print(rxn.gene_reaction_rule) to actually see the rule. This is also super consistent with Reaction.

I would have like to split the PR to parts, where the first one adds GPR to cobra.core.gene and doesn't change almost anything (although some test might fail) and the second part removes all irrelevant functions so ast2str, parse_gpr, eval_gpr become internal functions.

I would structure the first PR so that it does pass all existing tests and moves the functions into the GPR class. So basically the same API but using the new GPR class in the background. If you have more stuff added and not pushed to your github repo you can just revert the commit with --soft so it keeps the changes. Then add only he parts you need to a new branch. Or copy the code to a new directory, burn everything down with git reset --hard and copy in the parts you want in a new branch.

EDIT: sorry, misread a part :)

[Midnighter]

You can also overwrite GPR.__repr__ if you want it to look nice when writing rxn.gene_reaction_rule but I think it's enough to overwrite GPR.__str__ and then call print(rxn.gene_reaction_rule) to actually see the rule. This is also super consistent with Reaction.

While I agree that this would be the right way to go, I think we have to maintain backwards compatibility here and return the full string. Something like:

class Reaction:

    ...

    @property
    def gene_reaction_rule(self):
        return str(self._gpr)

[akaviaLab]

I made GPR a class, and gave it its own pull request. Please review
#1134

This is exactly like Midnighter says, to keep compatibility.

I'm going to keep working on it while it is reviewed to shift all the code to use the new class.

[akaviaLab]

One general question.
I'd like to move ast2str and eval_gpr into the GPR() class. However, since it is recursive, I'm not sure how to do it.

I've tried to do it, but it throws up warnings and errors about unexpected inputs, like when trying to run self.eval_gpr(self.body, knockouts). Do you have advice or an example?

[cdiener]

Yeah that is a bit problematic. You could do a private method plus a public one like.

class GPR(Module):
    def __eval_gpr(self, tree, knockouts):
       # just always call the recursions as self.__eval_gpr(a, b) 
       ...
    
    def eval(self, knockouts):
        return self.__eval_gpr(self.body, knockouts)

or with a static method:

class GPR(Module):
    @staticmethod
    def eval_gpr(tree, knockouts):
        ...
    
    def eval(self, knockouts):
        return GPR.eval_gpr(self.body, knockouts)

Not sure what looks cleaner. Maybe (1) because it does not expose an extra method on the class?

Or just leave them outside of the class. Maybe the best for backwards compatibility. In that case maybe factor all the GPR code into its own source file like core/gpr.py.

[akaviaLab]

Okay. The PR should be updated with all comments and commits. I've removed SYMPY from this branch for now.

I am also wondering if gene.remove_from_model() will be removed, since I tried to get it to work with the new GPR, but had no luck

[Midnighter]

I'm closing this discussion due to the work by @akaviaLab that has been merged #1134.