Merge pull request #11 from MetExplore/development

Version 0.9.0-beta.0

docs/examples/fastcore.md

@@ -10,13 +10,11 @@ Fastcore defines a set of core reactions that is forced to be active (carry a no
 The Fastcore algorithm is greedy approximation of the exact problem which can be modelled as a MIP problem. With MIOM, the exact problem can be defined in a few lines, using the `opt_tol` parameter to control the level of optimality required:
 
 ```python
-from miom import miom, Solvers
-from miom.mio import load_gem
-from miom.miom import Comparator, ExtractionMode
+import miom
 import numpy as np
 
 # Use the flux-consistent subnetwork (fcm) of the Human1 GEM model 
-m = load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/homo_sapiens_human1_fcm.miom')
+m = miom.mio.load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/homo_sapiens_human1_fcm.miom')
 # Select reactions from the cholesterol metabolism as the core reactions to keep
 core_rxn = m.find_reactions_from_pathway("Cholesterol metabolism")
 print(sum(core_rxn))
@@ -26,18 +24,19 @@ weights = -1 * np.ones(m.num_reactions)
 weights[core_rxn == 1] = 1
 
 # Exact-Fastcore
-fmc = (miom(m, solver=Solvers.GUROBI_PYMIP)
-          .setup(opt_tol=0.01)
-          .steady_state()
-          .subset_selection(weights)
-          .keep(np.where(core_rxn == 1)[0])
-          .solve(verbosity=1)
-          .select_subnetwork(
-               mode=ExtractionMode.ABSOLUTE_FLUX_VALUE,
-               comparator=Comparator.GREATER_OR_EQUAL,
-               value=1e-8
-          )
-          .network
+fmc = (miom
+        .load(m, solver=miom.Solvers.GUROBI_PYMIP)
+        .setup(opt_tol=0.01)
+        .steady_state()
+        .subset_selection(weights)
+        .keep(core_rxn == 1)
+        .solve(verbosity=1)
+        .select_subnetwork(
+            mode=miom.ExtractionMode.ABSOLUTE_FLUX_VALUE,
+            comparator=miom.Comparator.GREATER_OR_EQUAL,
+            value=1e-8
+        )
+        .network
 )
 print(fmc.num_reactions)
 ```

docs/examples/fba.md

@@ -3,10 +3,11 @@
 A simple Flux Balance Analysis can be easily defined and solved with any of the commercial or open-source solvers available:
 
 ```python
-from miom import miom, load_gem
-network = load_gem("https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom")
+import miom
+network = miom.mio.load_gem("https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom")
 target = "BIOMASS_reaction"
-flux = (miom(network)
+flux = (miom
+        .load(network)
         .steady_state()
         .set_rxn_objective(target)
         .solve(verbosity=1)

docs/examples/imat.md

@@ -6,12 +6,10 @@ Example implementation of iMAT using MIOM. Note that this implementation support
     Shlomi, T., Cabili, M. N., Herrgård, M. J., Palsson, B. Ø., & Ruppin, E. (2008). Network-based prediction of human tissue-specific metabolism. [Nature biotechnology, 26(9), 1003-1010](https://www.nature.com/articles/nbt.1487).
 
 ```python
-from miom import miom, Solvers
-from miom.mio import load_gem
-from miom.miom import Comparator, ExtractionMode
+import miom
 
 # Use the iHuman-GEM model
-m = load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/homo_sapiens_human1.miom')
+m = miom.mio.load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/homo_sapiens_human1.miom')
 
 # Add all the reactions from the Cholesterol pathway to the highly expressed set
 RH = m.find_reactions_from_pathway("Cholesterol metabolism")
@@ -20,14 +18,15 @@ RL = -1 * m.find_reactions_from_pathway("Pyruvate metabolism")
 w = RH + RL
 print("RH:", sum(RH), "RL:", sum(abs(RL)))
 
-m = (miom(m, solver=Solvers.GUROBI)
+m = (miom
+     .load(m, solver=miom.Solvers.GUROBI)
      .setup(int_tol=1e-8, opt_tol=0.01, verbosity=1)
      .steady_state()
      .subset_selection(w)
      .solve(max_seconds=30)
      .select_subnetwork(
-          mode=ExtractionMode.ABSOLUTE_FLUX_VALUE,
-          comparator=Comparator.GREATER_OR_EQUAL,
+          mode=miom.ExtractionMode.ABSOLUTE_FLUX_VALUE,
+          comparator=miom.Comparator.GREATER_OR_EQUAL,
           value=1e-8
      )
      .network)

docs/examples/sparse_fba.md

@@ -7,18 +7,19 @@ Toolbox includes different LP heuristics to minimize different approximations of
 Here is the implementation of sparse-FBA with MIOM:
 
 ```python
-from miom import miom, load_gem, Solvers
+import miom
 
 # Load a genome-scale metabolic network. You can load SMBL or Matlab metabolic networks
 # as well using the same method, but it requires to have the cobratoolbox python library
 # installed.
-network = load_gem("https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom")
+network = miom.mio.load_gem("https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom")
 
 # Create the sparse FBA problem to get a solution that maximizes
 # the optimal flux through the BIOMASS_reaction minimizing the
 # number of active reactions. The solution should be not more than
 # 5% of the optimal solution (opt_tol = 0.05).
-V, X = (miom(network, solver=Solvers.GUROBI_PYMIP)
+V, X = (miom
+        .load(network, solver=miom.Solvers.GUROBI_PYMIP)
         # Set-up the solver options
         .setup(int_tol=1e-8, opt_tol=0.05, verbosity=1)
         # Add the steady-state constraints (S*V = 0)
@@ -51,7 +52,7 @@ V, X = (miom(network, solver=Solvers.GUROBI_PYMIP)
         .get_values())
 
 # Show reactions with a flux > 1e-7
-print("Number of reactions with flux above +/- 1e-7:", sum(abs(V)>1e-7))
+print("Number of reactions with flux above +/- 1e-8:", sum(abs(V) > 1e-8))
 
 # Count reactions with an indicator value of 0 (active). Note that since
 # the weights of the reactions are negative (for all rxns), an indicator

docs/index.md

@@ -5,12 +5,10 @@
 
 
 
-__MIOM__ (Mixed Integer Optimization for Metabolism) is a python library for creating and solving complex optimization problems using genome-scale metabolic networks, in just a few lines. 
+__MIOM__ (Mixed Integer Optimization for Metabolism) is a python library for creating and solving complex constraint-based optimization problems for metabolism, in just a few lines. 
 
 By leveraging the power of modern Mixed Integer Optimization (MIO) solvers, it can transform any simple constraint-based problem, such as Flux Balance Analysis (FBA) into a more complex optimization problem, such as sparse FBA or context-specific reconstruction problems very easily, and solve it with the __required level of optimality__.
 
-But what is even better is that most of the time, algorithms formulated as Mixed Integer Optimization problems with MIOM can be solved faster and with better quality than currently existing alternatives that are approximations of the original problem. By using the MIO formulation, you can get also an estimation of how close to optimality a solution is, so you don't need to waste more time than needed.
-
 MIOM uses the [PICOS](https://picos-api.gitlab.io/picos/) and the [Python-MIP](https://www.python-mip.com/) libraries to build and solve the optimization problems using many commercial, academic and free solvers.
 
 !!! warning
@@ -37,17 +35,18 @@ CPLEX is also supported, but requires a license. To install MIOM with CPLEX supp
 Here is an example of how to load a metabolic network and maximize the flux through a target reaction using FBA, and then how to modify the original problem to implement the sparse FBA problem adding only a few lines to the original problem:
 
 ```python
-from miom import miom, load_gem, Solvers
+import miom
 
 # Load a genome-scale metabolic network using the miom format. 
 # You can load SMBL or Matlab metabolic networks as well using 
 # the same method, but it requires to have the cobratoolbox python library
 # installed (and scipy for mat files). To install these dependencies, run:
 # $ pip install cobra scipy
-network = load_gem("https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom")
+network = miom.mio.load_gem("https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom")
 target_rxn = "BIOMASS_reaction"
 # Create the optimization problem with miom and solve
-model = (miom(network)
+model = (miom
+        .load(network)
         .steady_state()
         .set_rxn_objective(target_rxn)
         .solve(verbosity=1))
@@ -109,7 +108,8 @@ approximate solution, controlled by the `opt_tol` parameter.
 To use other solvers, you only need to provide the specific solver to the `miom` method, for example:
 
 ```python
-model = (miom(network, solver=Solvers.GLPK)
+model = (miom
+        .load(network, solver=Solvers.GLPK)
         .steady_state()
         .set_rxn_objective(target_rxn)
         .solve(verbosity=1))

examples/exact_fastcore_example.py

@@ -1,27 +1,27 @@
-from miom import miom, Solvers
-from miom.mio import load_gem
-from miom.miom import Comparator, ExtractionMode
+import miom
 import numpy as np
 
 # Load a model
-m = load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/homo_sapiens_human1_fcm.miom')
+m = miom.mio.load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/homo_sapiens_human1_fcm.miom')
 print("Reactions in the consistent iHuman GEM:", m.num_reactions)
 core_rxn = m.find_reactions_from_pathway("Cholesterol metabolism")
 print("Num. of core reactions:", sum(core_rxn))
-# Assign a negative weight for reactions not in the core
+
+# Assign a negative weight for reactions not in the core and a positive value
+# for reactions in the core.
 weights = -1 * np.ones(m.num_reactions)
 weights[core_rxn == 1] = 1
 
-# Fastcore
-fmc = (miom(m, solver=Solvers.GUROBI_PYMIP)
+# Weighted Exact Fastcore algorithm with MIOM:
+fmc = (miom.load(m, solver=miom.Solvers.COIN_OR_CBC)
           .setup(opt_tol=0.05)
           .steady_state()
           .subset_selection(weights)
-          .keep(np.where(core_rxn == 1)[0])
+          .keep(core_rxn == 1)
           .solve(verbosity=1)
           .select_subnetwork(
-               mode=ExtractionMode.ABSOLUTE_FLUX_VALUE,
-               comparator=Comparator.GREATER_OR_EQUAL,
+               mode=miom.ExtractionMode.ABSOLUTE_FLUX_VALUE,
+               comparator=miom.Comparator.GREATER_OR_EQUAL,
                value=1e-8
           )
           .network

examples/fba_example.py

@@ -1,11 +1,10 @@
-from miom import load_gem, miom, Solvers
-import numpy as np
+import miom
 
 # Load a model
-m = load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom')
+m = miom.mio.load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom')
 
 target = 'BIOMASS_reaction'
-opt_flux = (miom(network=m, solver=Solvers.GLPK)
+opt_flux = (miom.load(network=m, solver=miom.Solvers.GLPK)
                 .steady_state()
                 .set_rxn_objective(target)
                 .solve(verbosity=1)

examples/imat_example.py

@@ -18,7 +18,7 @@
 w = RH + RL
 print("RH:", sum(RH), "RL:", sum(abs(RL)))
 
-m = (miom(m, solver=Solvers.GUROBI)
+m = (miom(m, solver=Solvers.COIN_OR_CBC)
      .setup(int_tol=1e-8, opt_tol=0.01, verbosity=1)
      .steady_state()
      .subset_selection(w)
@@ -29,5 +29,5 @@
           value=1e-8
      )
      .network)
-print(m.num_reactions)
+print("Number of reactions in the subnetwork with active flux:", m.num_reactions)
 

examples/introduction_example.py

@@ -1,4 +1,5 @@
-from miom import miom, load_gem, Solvers
+from miom import miom, Solvers
+from miom.mio import load_gem
 
 # Load a genome-scale metabolic network. You can load SMBL or Matlab metabolic networks
 # as well using the same method, but it requires to have the cobratoolbox python library
@@ -11,8 +12,9 @@
         .set_rxn_objective(target_rxn)
         .solve(verbosity=1))
 print("Optimal flux:", model.get_fluxes(target_rxn), "mmol/(h·gDW)")
+
 # Show reactions with non-zero flux
-V, _ = model.get_values()
+V = model.get_fluxes()
 print("Number of reactions with flux above +/- 1e-8:", sum(abs(V) > 1e-8))
 
 V, X = (model
@@ -31,6 +33,6 @@
         .get_values())
 
 print("Number of reactions with an absolute flux value above 1e-8:", sum(abs(V) > 1e-8))
-print("Active reactions:", sum(1 if activity == 1 else 0 for activity in model.variables.indicator_rxn_activity))
-print("Inconsistencies:", sum(1 if activity != activity else 0 for activity in model.variables.indicator_rxn_activity))
+print("Active reactions:", sum(1 if activity == 1 else 0 for activity in model.variables.reaction_activity))
+print("Inconsistencies:", sum(1 if activity != activity else 0 for activity in model.variables.reaction_activity))
 print("Solver status:", model.get_solver_status())

examples/sparse_fba_example.py

@@ -1,11 +1,12 @@
-from miom import miom, load_gem, Solvers
+import miom
 
 # Load a model
-m = load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom')
-print(m.num_reactions)
+m = miom.mio.load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom')
+print("Num. of reactions in the network:", m.num_reactions)
 
 # Solve with Gurobi, CPLEX or CBC (other MIP solvers struggle with mediudm/large networks)
-V, X = (miom(network=m, solver=Solvers.GUROBI_PYMIP)
+V, X = (miom
+        .load(network=m, solver=miom.Solvers.COIN_OR_CBC)
         # Config the solver (e.g. set the optimality tolerance for MIP problems)
         .setup(int_tol=1e-8, opt_tol=0.05)
         # Add steady-state constraints to the model (S * V = 0)
@@ -30,7 +31,7 @@
         .get_values())
 
 # Show reactions with a flux > 1e-7
-print("Number of reactions with flux above +/- 1e-7:", sum(abs(V)>1e-7))
+print("Number of reactions with flux above +/- 1e-8:", sum(abs(V)>1e-8))
 
 # Count reactions with an indicator value of 0 (active). Note that since
 # the weights of the reactions are negative (for all rxns), an indicator

miom/__init__.py

@@ -1,11 +1,9 @@
-from miom.miom import (
-    miom,
-    PicosModel, 
-    PythonMipModel, 
-    Solvers
+from miom.miom import ( 
+    load,
+    Solvers,
+    Comparator,
+    ExtractionMode
 )
 
-from miom.mio import (
-    load_gem, 
-    export_gem
-)
+__all__ = ["load", "Solvers", "Comparator", "ExtractionMode"]
+__version__ = "0.9.0-beta.0"

miom/__main__.py

@@ -0,0 +1,35 @@
+import argparse
+import miom
+
+def convert_gem(args):
+    print(f"MIOM v{miom.__version__}")
+    m = miom.mio.load_gem(args.input)
+    print(f"Loaded network with {m.num_reactions} reactions")
+    print(f"Exporting to {args.output}...")
+    miom.mio.export_gem(m, args.output)
+    print("Done.")
+
+def get_args():
+    parser = argparse.ArgumentParser(description=f"MIOM: Mixed Integer Optimization for Metabolism")
+    parser.add_argument(
+        "--version",
+        action="version", 
+        version=f"MIOM v{miom.__version__}"
+    )
+    subparsers = parser.add_subparsers(help="sub-command help")
+    convert = subparsers.add_parser("convert", help="Convert a model to miom format")
+    convert.add_argument(
+        "input",
+        help="Input model file (if cobra is installed, any format supported by cobra is allowed)"
+    )
+    convert.add_argument(
+        "output",
+        help="Output GEM file in MIOM format"
+    )
+    convert.set_defaults(func=convert_gem)
+    return parser.parse_args()
+
+if __name__ == '__main__':
+    args = get_args()
+    if args.func:
+        args.func(args)

miom/miom.py

@@ -148,7 +148,7 @@ class Comparator(str, Enum):
     ])
 
 
-def miom(network, solver=Solvers.COIN_OR_CBC):
+def load(network, solver=Solvers.COIN_OR_CBC):
     """
     Create a MIOM optimization model for a given solver.
     If the solver is Coin-OR CBC, an instance of PythonMipModel is used
@@ -157,12 +157,13 @@ def miom(network, solver=Solvers.COIN_OR_CBC):
 
     Example:
         Example of how to perform FBA to maximize flux through the
-        BIOMASS_reaction in the iMM1865 model:
+        `BIOMASS_reaction` in the iMM1865 model:
 
         ```python
-        >>> from miom import miom, load_gem
-        >>> network = load_gem("https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom")
-        >>> V, _ = miom(network)
+        >>> import miom
+        >>> network = miom.mio.load_gem("https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom")
+        >>> V, X = (miom
+                    .load(network)
                     .steady_state()
                     .set_rxn_objective("BIOMASS_reaction")
                     .solve(verbosity=1)
@@ -480,6 +481,9 @@ def keep(self, reactions):
             return False
         if isinstance(reactions, str):
             reactions = [reactions]
+        # Check if it's a binary vector
+        if len(reactions) == self.network.num_reactions and max(reactions)==1:
+            reactions = np.where(reactions==1)[0]
         # Check if the reactions have an indicator variable
         reactions = set(self.network.find_reaction(rxn)[0] for rxn in reactions)
         available = set(rxn.index for rxn in self.variables.assigned_reactions if rxn.cost > 0)

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "miom"
-version = "0.9.0-alpha.4"
+version = "0.9.0-beta.0"
 description = "Mixed Integer Optimization for Metabolism"
 authors = ["Pablo R. Mier <[email protected]>"]
 keywords = ["optimization", "LP", "MIP", "metabolism", "metabolic-networks"]