Merge pull request #10 from cdanielmachado/experimental_constraints

Experimental constraints

carveme/reconstruction/carving.py

@@ -1,4 +1,5 @@
 import numpy as np
+import warnings
 import pandas as pd
 
 from framed.cobra.ensemble import EnsembleModel, save_ensemble
@@ -33,7 +34,8 @@ def inactive_reactions(model, solution):
     return inactive + inactive_ext
 
 
-def minmax_reduction(model, scores, min_growth=0.1, min_atpm=0.1, eps=1e-3, bigM=1e3, default_score=-1, uptake_score=0, solver=None):
+def minmax_reduction(model, scores, min_growth=0.1, min_atpm=0.1, eps=1e-3, bigM=1e3, default_score=-1, uptake_score=0,
+                     soft_score=1, soft_constraints=None, hard_constraints=None, solver=None, debug_output=None):
     """ Apply minmax reduction algorithm (MILP).
 
     Computes a binary reaction vector that optimizes the agreement with reaction scores (maximizes positive scores,
@@ -47,8 +49,11 @@ def minmax_reduction(model, scores, min_growth=0.1, min_atpm=0.1, eps=1e-3, bigM
         min_atpm (float): minimal maintenance ATP constraint
         eps (float): minimal flux required to consider leaving the reaction in the model
         bigM (float): maximal reaction flux
-        default_score (float): penalty score for reactions without an annotation score (default: -1).
-        uptake_score (float): penalty score for using uptake reactions (default: -1).
+        default_score (float): penalty score for reactions without an annotation score (default: -1.0).
+        uptake_score (float): penalty score for using uptake reactions (default: 0.0).
+        soft_score (float): score for soft constraints (default: 1.0)
+        soft_constraints (dict): dictionary from reaction id to expected flux direction (-1, 1, 0)
+        hard_constraints (dict): dictionary of flux bounds
         solver (Solver): solver instance (optional)
 
     Returns:
@@ -63,6 +68,14 @@ def minmax_reduction(model, scores, min_growth=0.1, min_atpm=0.1, eps=1e-3, bigM
     scores = scores.copy()
     reactions = scores.keys()
 
+    if not soft_constraints:
+        soft_constraints = {}
+
+    reactions += [r_id for r_id in soft_constraints if r_id not in reactions]
+
+    if hard_constraints:
+        solver.set_bounds(hard_constraints)
+
     # Add default score
     if default_score != 0:
         for r_id in model.reactions:
@@ -119,27 +132,48 @@ def minmax_reduction(model, scores, min_growth=0.1, min_atpm=0.1, eps=1e-3, bigM
 
     for r_id in reactions:
         y_r, y_f = 'yr_' + r_id, 'yf_' + r_id
+
+        if r_id in soft_constraints:
+            sign = soft_constraints[r_id]
+
+            if sign > 0:
+                w_f, w_r = soft_score, 0
+            elif sign < 0:
+                w_f, w_r = 0, soft_score
+            else:
+                w_f, w_r = -soft_score, -soft_score
+
         if y_f in solver.pos_vars:
-            objective[y_f] = scores[r_id]
+            if r_id in scores:
+                objective[y_f] = scores[r_id]
+            if r_id in soft_constraints:
+                objective[y_f] = w_f
+
         if y_r in solver.neg_vars:
-            objective[y_r] = scores[r_id]
+            if r_id in scores:
+                objective[y_r] = scores[r_id]
+            if r_id in soft_constraints:
+                objective[y_r] = w_r
 
     if uptake_score != 0:
         for r_id in model.reactions:
-            if r_id.startswith('R_EX'):
+            if r_id.startswith('R_EX') and r_id not in soft_constraints:
                 objective['y_' + r_id] = uptake_score
 
     solver.set_objective(linear=objective, minimize=False)
-    solution = solver.solve()
 
-#    for r_id in model.reactions:
-#        print r_id, solution.values[r_id], solution.values.get('yr_'+r_id, None), solution.values.get('yf_'+r_id, None)
+    if debug_output:
+        solver.write_to_file(debug_output + "_milp_problem.lp")
+
+    solution = solver.solve()
 
     return solution
 
 
 def carve_model(model, reaction_scores, outputfile=None, flavor=None, inplace=True,
-                default_score=-1.0, uptake_score=0.0, init_env=None, debug_output=None):
+                default_score=-1.0, uptake_score=0.0, soft_score=1.0,
+                soft_constraints=None, hard_constraints=None,
+                init_env=None, debug_output=None):
     """ Reconstruct a metabolic model using the CarveMe approach.
 
     Args:
@@ -150,6 +184,9 @@ def carve_model(model, reaction_scores, outputfile=None, flavor=None, inplace=Tr
         inplace (bool): Change model in place (default: True)
         default_score (float): penalty for non-annotated intracellular reactions (default: -1.0)
         uptake_score (float): penalty for utilization of extracellular compounds (default: 0.0)
+        soft_score (float): score for soft constraints (default: 1.0)
+        soft_constraints (dict): dictionary from reaction id to expected flux direction (-1, 1, 0)
+        hard_constraints (dict): dictionary of flux bounds
         init_env (Environment): initialize final model with given Environment (optional)
 
     Returns:
@@ -167,7 +204,24 @@ def carve_model(model, reaction_scores, outputfile=None, flavor=None, inplace=Tr
     if uptake_score is None:
         uptake_score = 0.0
 
-    sol = minmax_reduction(model, scores, default_score=default_score, uptake_score=uptake_score)
+    if soft_score is None:
+        soft_score = 1.0
+
+    if soft_constraints:
+        not_in_model = set(soft_constraints) - set(model.reactions)
+        if not_in_model:
+            soft_constraints = {r_id: val for r_id, val in soft_constraints.items() if r_id in model.reactions}
+            warnings.warn("Soft constraints contain reactions not in the model:\n" + "\n".join(not_in_model))
+
+    if hard_constraints:
+        not_in_model = set(hard_constraints) - set(model.reactions)
+        if not_in_model:
+            hard_constraints = {r_id: (lb, ub) for r_id, (lb, ub) in hard_constraints.items() if r_id in model.reactions}
+            warnings.warn("Hard constraints contain reactions not in the model:\n" + "\n".join(not_in_model))
+
+    sol = minmax_reduction(model, scores, default_score=default_score, uptake_score=uptake_score, soft_score=soft_score,
+                           soft_constraints=soft_constraints, hard_constraints=hard_constraints,
+                           debug_output=debug_output)
 
     if sol.status == Status.OPTIMAL:
         inactive = inactive_reactions(model, sol)

carveme/reconstruction/gapfilling.py

@@ -69,8 +69,8 @@ def gapFill(model, universe, constraints=None, min_growth=0.1, scores=None, inpl
         inactive = [r_id for r_id in new_reactions if abs(solution.values[r_id]) < abstol]
 
     else:
-        inactive = new_reactions
-        print 'Failed to gapfill model', ('for medium {}'.format(tag) if tag else '')
+#        inactive = new_reactions
+        raise RuntimeError('Failed to gapfill model for medium {}'.format(tag))
 
     model.remove_reactions(inactive)
     del_metabolites = disconnected_metabolites(model)
@@ -80,7 +80,8 @@ def gapFill(model, universe, constraints=None, min_growth=0.1, scores=None, inpl
         return model
 
 
-def multiGapFill(model, universe, media, media_db, min_growth=0.1, max_uptake=10, scores=None, inplace=True, bigM=1e3):
+def multiGapFill(model, universe, media, media_db, min_growth=0.1, max_uptake=10, scores=None, inplace=True, bigM=1e3,
+                 exchange_format=None):
     """ Gap Fill a metabolic model for multiple environmental conditions
 
     Args:
@@ -116,7 +117,8 @@ def multiGapFill(model, universe, media, media_db, min_growth=0.1, max_uptake=10
     for medium_name in media:
         if medium_name in media_db:
             compounds = media_db[medium_name]
-            constraints = medium_to_constraints(merged_model, compounds, max_uptake=max_uptake, inplace=False)
+            constraints = medium_to_constraints(merged_model, compounds, max_uptake=max_uptake, inplace=False,
+                                                exchange_format=exchange_format, verbose=False)
 
             gapFill(model, universe, constraints=constraints, min_growth=min_growth,
                     scores=scores, inplace=True, bigM=bigM, solver=solver, tag=medium_name)

carveme/reconstruction/utils.py

@@ -45,26 +45,6 @@ def add_maintenance_atp(model, lb=0, ub=1000):
     model.add_reaction_from_str(rxn_str)
 
 
-# def simulate_biomass_vs_medium(model, biomass_dict, media_db, constraints=None):
-#
-#     result = {}
-#
-#     for organism, equation in biomass_dict.items():
-#         result[organism] = {}
-#         r_id = add_biomass_equation(model, equation, organism)
-#
-#         for medium_id, compounds in media_db.items():
-#             constr = medium_to_constraints(model, compounds)
-#             if constraints:
-#                 constr.update(constraints)
-#             sol = FBA(model, constraints=constr, objective={r_id: 1})
-#             result[organism][medium_id] = sol.fobj
-#
-#         model.remove_reaction(r_id)
-#
-#     return pd.DataFrame(result).T
-
-
 def tab2fasta(inputfile, outputfile, filter_by_model=None):
 
     data = pd.read_csv(inputfile, sep='\t', header=0)
@@ -86,9 +66,12 @@ def load_media_db(filename, sep='\t', medium_col='medium', compound_col='compoun
     return media_db[compound_col].to_dict()
 
 
-def medium_to_constraints(model, compounds, max_uptake=10, inplace=False, verbose=False):
+def medium_to_constraints(model, compounds, max_uptake=10, inplace=False, verbose=False, exchange_format=None):
 
-    env = Environment.from_compounds(compounds, max_uptake=max_uptake)
+    if not exchange_format:
+        exchange_format = "'R_EX_{}_e'"
+
+    env = Environment.from_compounds(compounds, max_uptake=max_uptake, exchange_format=exchange_format)
     return env.apply(model, inplace=inplace, warning=verbose)
 
 
@@ -142,3 +125,14 @@ def add_biomass_equation(model, stoichiometry, label=None):
     name = 'Biomass reaction'
     reaction = CBReaction(r_id, name=name, reversible=False, stoichiometry=stoichiometry, objective=1.0)
     model.add_reaction(reaction)
+
+
+def load_soft_constraints(filename):
+    df = pd.read_csv(filename, sep='\t', header=None)
+    return dict(zip(df[0], df[1]))
+
+
+def load_hard_constraints(filename):
+    df = pd.read_csv(filename, sep='\t', header=None)
+    return dict(zip(df[0], zip(df[1], df[2])))
+