Feat equilibrator thermodata (#38)

* Dev: Incremental update (#29)

* MNT: Pypi details

* MNT: Adding veriable creation logs

* FIX: typo

* FIX: Issue #28, incorrect variable referencing

* TST: adding tests for some of the analysis functions

* FIX: more fixes from #28

* FIX: fixed transport reaction function

* TST: Better tests, with a smaller model for I/O

* TST/FIX: relative<- absolute paths for small_model resources

* FIX: typo in paths to small_model resources for tests

* VER: Bump v0.9.0-b2

* VER: bump to v0.9.1

* Update LICENSE.txt

* Fixing a couple of minor bugs and added new variable class (#35)

Thank you @remidhum 

* FIX: fixed the apply_directionality function

solution.raw dataframe index does not contain the fwd and bwd use variables.

* ENH: added new binary variable class

This class is usefull to deal with inactive reactions in a model.

* ENH: deal with models without an objective function

Objective function is set to Zero (symbol("0") does not work!) if there is no defined objective function

* OOPS: change proper function ...

* FIX: fixed the failing test

As suggested, I added a check to determine what object type is passed.

* MNT: solution object type testing improved

* feat: TOO SLOW eQuilibrator integration

* refactor: make use of new changes on equilibrator_api (https://gitlab.com/elad.noor/equilibrator-api/merge_requests/13)

* chore: equilibrator_api as extra dependency

* fix: make thermo_data optional

* chore: test equilibrator_api integration

* chore: set equilibrator_api as extra dependency

* refactor: back to previous init to impose thermo_data argument

* refactor: build thermo data from equilibrator as external function

* fix: install equilibrator as package

Co-Authored-By: Moritz E. Beber <midnighter@posteo.net>

* FIX: correct typo on error attribute

* refactor: load equilibrator on import of file

* REFACTOR: move print to logger warning

* FIX: set lower bound version on equilibrator_cache compatibility

* FIX: remove f-string syntax

* FIX: test equilibrator compatibility only for >= py3.6

* FIX: try to overcome py3.5 problems

* REFACTOR: ignore equilibrator test collection on py<3.6

* FIX: collect tests properly

* DOCS: add script as tutorial for equilibrator integration

* refactor: move equilibrator test to another module

* REFACTOR: encapsulate equilibrator testing

* "REFACTOR: REVERT encapsulate equilibrator testing, it doesn't work"

This reverts commit edc240d.

* "REFACTOR: REVERT move equilibrator test to another module"

This reverts commit 2f256de.

* Update test_equilibrator.py

* Update conftest.py

* fix: require Py>=3.6 when for equilibrator extra

* fix: import on function

* fix: remove conftest, proper coverage

* FIX: move pytfa import to function

* DOCS: warn about python version for equilibrator integration

* REFACTOR: remove GLPK solver enforcement

* fix: avoid side effects on tmodel

* FIX: change equilibrator-cache version

* TEST: change extra dependency management before_install

* FIX: typo

* FIX: add equilibrator_api to build

* fix: change travis pytest command to properly account for coverage

* FIX: typo

* FIX: separate coverage from pytest

* FIX: call codecov

* FIX: give coverage a try

* FIX: correct coverage install and calls

* Revert "FIX: correct coverage install and calls"

This reverts commit c1dfe77.

* Revert "FIX: give coverage a try"

This reverts commit eda6365.

* Revert "FIX: call codecov"

This reverts commit 38399e3.

* Revert "FIX: separate coverage from pytest"

This reverts commit 43ba69b.

* Revert "FIX: typo"

This reverts commit f10f1f8.

* Revert "fix: change travis pytest command to properly account for coverage"

This reverts commit 2b7b257.

* FIX: collect equilibrator tests at the end

* FIX: try with edit mode

* TEST: change file name

Co-authored-by: realLCSB <31034133+realLCSB@users.noreply.github.com>
Co-authored-by: Pierre Salvy <psalvy@users.noreply.github.com>
Co-authored-by: RémiDhum <rdhumeaux@gmail.com>
Co-authored-by: Moritz E. Beber <midnighter@posteo.net>

.travis.yml

@@ -1,16 +1,23 @@
 language: python
 dist: xenial
-python:
-  - "3.5"
-  - "3.6"
-  - "3.7"
+jobs:
+ include:
+ - python: 3.5
+   before_install:
+    - pip install pytest-cov
+ - python: 3.6
+   before_install:
+    - pip install pytest-cov
+    - pip install equilibrator_cache>=0.2.6 equilibrator_api
+ - python: 3.7
+   before_install:
+    - pip install pytest-cov
+    - pip install equilibrator_cache>=0.2.6 equilibrator_api
 branches:
  only:
  - master
  - dev
 
-before_install:
- - pip install pytest-cov
 install:
  - pip install -e .
 

LICENSE.txt

@@ -186,7 +186,7 @@
       same "printed page" as the copyright notice for easier
       identification within third-party archives.
 
-   Copyright {yyyy} {name of copyright owner}
+   Copyright 2019 Laboratory of Computational Systems Biotechnology, EPFL, Switzerland 
 
    Licensed under the Apache License, Version 2.0 (the "License");
    you may not use this file except in compliance with the License.

doc/conf.py

@@ -95,7 +95,7 @@ def __getattr__(cls, name):
 version = '0.9'
 # The full version, including alpha/beta/rc tags.
 
-release = '0.9.0-b2'
+release = '0.9.1'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

pytfa/analysis/manipulation.py

@@ -1,3 +1,6 @@
+from ..core.model import Solution
+import pandas as pd
+
 def apply_reaction_variability(tmodel, va, inplace = True):
     """
     Applies the VA results as bounds for the reactions of a cobra_model
@@ -46,7 +49,7 @@ def apply_generic_variability(tmodel,va, inplace = True):
     return _tmodel
 
 
-def apply_directionality(tmodel,solution, inplace = True):
+def apply_directionality(tmodel, solution, inplace = True):
     """
     Takes a flux solution and transfers its reaction directionality as
     constraints for the cobra_model
@@ -62,15 +65,22 @@ def apply_directionality(tmodel,solution, inplace = True):
     else:
         _tmodel = tmodel.copy()
 
+    if isinstance(solution, Solution):
+        sol = solution.raw
+    elif isinstance(solution, pd.Series) or isinstance(solution, pd.DataFrame):
+        sol = solution
+    else:
+        raise ArgumentError('solution object should be of class Solution or pandas.Series')
+
     for this_reaction in _tmodel.reactions:
 
         backward_use = _tmodel.backward_use_variable.get_by_id(this_reaction.id)
         forward_use = _tmodel.forward_use_variable.get_by_id(this_reaction.id)
 
-        backward_use.variable.lb = round(solution.raw[backward_use.name])
-        backward_use.variable.ub = round(solution.raw[backward_use.name])
+        backward_use.variable.lb = round(sol[backward_use.name])
+        backward_use.variable.ub = round(sol[backward_use.name])
 
-        forward_use.variable.lb  = round(solution.raw[forward_use.name])
-        forward_use.variable.ub  = round(solution.raw[forward_use.name])
+        forward_use.variable.lb  = round(sol[forward_use.name])
+        forward_use.variable.ub  = round(sol[forward_use.name])
 
-    return _tmodel
+    return _tmodel

pytfa/io/dict.py

@@ -22,6 +22,8 @@
 
 from copy import copy
 
+import sympy
+
 
 def get_all_subclasses(cls):
     all_subclasses = []
@@ -371,8 +373,11 @@ def model_from_dict(obj, solver=None, custom_hooks = None):
     return new
 
 def rebuild_obj_from_dict(new, objective_dict):
-    obj_expr = symbol_sum([v*new.variables.get(x) for x,v in objective_dict.items()])
-    new.objective = obj_expr
+    if objective_dict.__class__ is dict:
+        obj_expr = symbol_sum([v*new.variables.get(x) for x,v in objective_dict.items()])
+        new.objective = obj_expr
+    else:
+        new.objective = sympy.S.Zero
 
 def add_custom_classes(model, custom_hooks):
     """
@@ -446,5 +451,3 @@ def _rebuild_stoichiometry(new, stoich):
     return defaultdict(int,
                        {new.metabolites.get_by_id(k):v
                         for k,v in stoich.items()})
-
-

pytfa/optim/variables.py

@@ -454,7 +454,6 @@ def model(self):
 
     prefix = 'MV_'
 
-
 class ForwardUseVariable(ReactionVariable, BinaryVariable):
     """
     Class to represent a forward use variable, a type of binary variable used to
@@ -486,6 +485,24 @@ def __init__(self, reaction, **kwargs):
 
     prefix = 'BU_'
 
+class ForwardBackwardUseVariable(ReactionVariable, BinaryVariable):
+    """
+    Class to represent a type of binary variable used to tell whether the
+    reaction is active or not such that:
+        FU + BU + BFUSE = 1
+    """
+
+    def __init__(self, reaction, **kwargs):
+        if not 'lb' in kwargs:
+            kwargs['lb'] = 0
+        if not 'ub' in kwargs:
+            kwargs['ub'] = 1
+
+        ReactionVariable.__init__(self, reaction,
+                                  type=get_binary_type(),
+                                  **kwargs)
+
+    prefix = 'BFUSE_'
 
 class LogConcentration(MetaboliteVariable):
     """
@@ -494,23 +511,20 @@ class LogConcentration(MetaboliteVariable):
 
     prefix = 'LC_'
 
-
 class DeltaGErr(ReactionVariable):
     """
     Class to represent a DeltaGErr
     """
 
     prefix = 'DGE_'
 
-
 class DeltaG(ReactionVariable):
     """
     Class to represent a DeltaG
     """
 
     prefix = 'DG_'
 
-
 class DeltaGstd(ReactionVariable):
     """
     Class to represent a DeltaG^o (naught) - standard conditions
@@ -536,7 +550,6 @@ def __init__(self,reaction,**kwargs):
 
     prefix = 'PosSlack_'
 
-
 class NegSlackVariable(ReactionVariable):
     """
     Class to represent a negative slack variable for relaxation problems
@@ -554,4 +567,3 @@ class PosSlackLC(MetaboliteVariable):
 class NegSlackLC(MetaboliteVariable):
 
     prefix = 'NegSlackLC_'
-

pytfa/thermo/equilibrator.py

@@ -0,0 +1,99 @@
+# -*- coding: utf-8 -*-
+"""Thermodynamic information for metabolites from eQuilibrator.
+
+.. module:: pytfa
+   :platform: Unix, Windows
+   :synopsis: Thermodynamics-based Flux Analysis
+
+.. moduleauthor:: pyTFA team
+"""
+
+import equilibrator_cache.compatibility as compat
+
+# this file is not imported in the __init__.py. The explicit import of the file
+# triggers the loading of eQuilbrator
+from equilibrator_api import ComponentContribution, Q_
+from equilibrator_api.phased_reaction import PhasedReaction
+from equilibrator_cache import create_compound_cache_from_quilt
+from .std import TEMPERATURE_0
+from ..utils.numerics import BIGM_DG
+from ..utils.logger import get_bistream_logger
+
+# load compound cache just one to aliviate resource usage
+ccache = None
+logger = get_bistream_logger("eQuilibrator_formation")
+
+
+def build_thermo_from_equilibrator(model, T=TEMPERATURE_0):
+    """Build `thermo_data` structure from a cobra Model.
+
+    The structure of the returned dictionary is specified in the pyTFA
+    [documentation](https://pytfa.readthedocs.io/en/latest/thermoDB.html).
+
+    :param model: cobra.Model
+    :return thermo_data: dict
+        to be passed as argument to initialize a `ThermoModel`.
+
+    """
+    global ccache
+    if ccache is None:
+        ccache = create_compound_cache_from_quilt()
+        logger.debug("eQuilibrator compound cache is loaded.")
+
+    cc = ComponentContribution(temperature=Q_(str(T) + "K"))
+
+    thermo_data = {"name": "eQuilibrator", "units": "kJ/mol", "cues": {}}
+    met_to_comps = compat.map_cobra_metabolites(ccache, model.metabolites)
+    thermo_data["metabolites"] = [
+        compound_to_entry(met, cc) for met in met_to_comps
+    ]
+    return thermo_data
+
+
+def compute_dGf(compound, cc):
+    """Get ΔGf from equilibrator `compound`."""
+    dG0_prime, dG0_uncertainty = cc.dG0_prime(
+        PhasedReaction(sparse={compound: 1}, rid="tmp_"+compound.id)
+    )
+    return dG0_prime, dG0_uncertainty
+
+
+def compound_to_entry(compound, cc):
+    """Build thermo structure entry from a `equilibrator_cache.Compound`.
+
+    eQuilibrator works with Component Contribution instead of groups, so it is
+    not possible to generate cues from it.
+
+    :param compound: equilibrator_cache.Compound
+    :return: dict
+        with keys ['deltaGf_std', 'deltaGf_err', 'error', 'struct_cues',
+        'id', 'pKa', 'mass_std', 'charge_std', 'nH_std', 'name', 'formula',
+        'other_names']
+
+    """
+    deltaGf_std, deltaGf_err = (BIGM_DG, BIGM_DG)
+    nH_std = compound.atom_bag["H"] if "H" in compound.atom_bag else 0
+    try:
+        deltaGf_std, deltaGf_err = compute_dGf(compound, cc)
+        err = "Nil"
+    except Exception as e:
+        err = 1
+        logger.debug(
+            "{} : thermo data NOT created, error : {}".format(
+                compound.id, e
+            )
+        )
+    return dict(
+        deltaGf_std=deltaGf_std,
+        deltaGf_err=deltaGf_err,
+        error=err,
+        struct_cues=None,
+        id=compound.id,
+        pKa=compound.dissociation_constants,
+        mass_std=compound.mass,
+        charge_std=None,
+        nH_std=nH_std,
+        name=compound.id,
+        formula=compound.formula,
+        other_names=compound.identifiers,
+    )