spent medium option in gapfill

.gitignore

@@ -1,26 +1,10 @@
-
-.idea/inspectionProfiles/Project_Default.xml
-
-.idea/carveme.iml
-
-.idea/misc.xml
-
-.idea/modules.xml
-
-.idea/vcs.xml
-
-.idea/workspace.xml
-
+.idea/*
 build/*
-
 dist/*
-
 carveme.egg-info/*
-
 *.pyc
-
 .eggs/*
-
 docs/_build/*
 Dockerfile
 .idea/encodings.xml
+*.DS_Store

carveme/reconstruction/gapfilling.py

@@ -4,6 +4,7 @@
 from framed.model.transformation import disconnected_metabolites
 from framed.solvers import solver_instance
 from framed.solvers.solver import VarType, Status
+from framed import FBA
 
 
 def gapFill(model, universe, constraints=None, min_growth=0.1, scores=None, inplace=True, bigM=1e3, abstol=1e-9,
@@ -71,7 +72,6 @@ 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
         raise RuntimeError('Failed to gapfill model for medium {}'.format(tag))
 
     model.remove_reactions(inactive)
@@ -83,7 +83,7 @@ def gapFill(model, universe, constraints=None, min_growth=0.1, scores=None, inpl
 
 
 def multiGapFill(model, universe, media, media_db, min_growth=0.1, max_uptake=10, scores=None, inplace=True, bigM=1e3,
-                 exchange_format=None):
+                 exchange_format=None, spent_model=None):
     """ Gap Fill a metabolic model for multiple environmental conditions
 
     Args:
@@ -96,6 +96,8 @@ def multiGapFill(model, universe, media, media_db, min_growth=0.1, max_uptake=10
         scores (dict): reaction scores (optional, see *gapFill* for details)
         inplace (bool): modify given model in place (default: True)
         bigM (float): maximal reaction flux (default: 1000)
+        exchange_format (str): format string to convert compounds to exchange reactions (default: "'R_EX_{}_e'")
+        spent_model (CBModel): additional species to generate spent medium compounds
 
     Returns:
         CBModel: gap filled model (if inplace=False)
@@ -104,6 +106,8 @@ def multiGapFill(model, universe, media, media_db, min_growth=0.1, max_uptake=10
         *media_db* is a dict from medium name to the list of respective compounds.
     """
 
+    ABSTOL = 1e-6
+
     if not inplace:
         model = model.copy()
 
@@ -116,12 +120,26 @@ def multiGapFill(model, universe, media, media_db, min_growth=0.1, max_uptake=10
     merged_model = merge_models(model, universe, inplace=False)
     solver = solver_instance(merged_model)
 
+    if spent_model:
+        solver0 = solver_instance(spent_model)
+
     for medium_name in media:
         if medium_name in media_db:
-            compounds = media_db[medium_name]
+            compounds = set(media_db[medium_name])
+
             constraints = medium_to_constraints(merged_model, compounds, max_uptake=max_uptake, inplace=False,
                                                 exchange_format=exchange_format, verbose=False)
 
+            if spent_model:
+                constraints0 = medium_to_constraints(spent_model, compounds, max_uptake=max_uptake, inplace=False,
+                                                     exchange_format=exchange_format, verbose=False)
+                for r_id in spent_model.get_exchange_reactions():
+                    if r_id in constraints:
+                        sol = FBA(spent_model, objective={r_id: 1}, constraints=constraints0, solver=solver0, get_values=False)
+                        if sol.fobj > ABSTOL:
+                            constraints[r_id] = (-max_uptake, None)
+                            print("added", r_id[5:-2], "to", medium_name)
+
             gapFill(model, universe, constraints=constraints, min_growth=min_growth,
                     scores=scores, inplace=True, bigM=bigM, solver=solver, tag=medium_name)
         else:

scripts/gapfill

@@ -10,7 +10,7 @@ import os
 
 
 def main(inputfile, media, mediadb=None, universe=None, universe_file=None, outputfile=None, flavor=None,
-         exchange_format=None, verbose=False):
+         exchange_format=None, spent=None, verbose=False):
 
     if verbose:
         print('Loading model...')
@@ -20,6 +20,17 @@ def main(inputfile, media, mediadb=None, universe=None, universe_file=None, outp
     except IOError:
         raise IOError('Failed to load model:' + inputfile)
 
+    if spent:
+        if verbose:
+            print('Loading model for spent medium species...')
+
+        try:
+            spent_model = load_cbmodel(spent, flavor=flavor)
+        except IOError:
+            raise IOError('Failed to load model:' + spent)
+    else:
+        spent_model = None
+
     if verbose:
         print('Loading reaction universe...')
 
@@ -54,7 +65,7 @@ def main(inputfile, media, mediadb=None, universe=None, universe_file=None, outp
 
     max_uptake = config.getint('gapfill', 'max_uptake')
     multiGapFill(model, universe_model, media, media_db, max_uptake=max_uptake, inplace=True,
-                 exchange_format=exchange_format)
+                 exchange_format=exchange_format, spent_model=spent_model)
 
     if verbose:
         m2, n2 = len(model.metabolites), len(model.reactions)
@@ -84,6 +95,9 @@ if __name__ == '__main__':
     parser.add_argument('-m', '--media', dest='media', required=True, help="List of media (comma-separated)")
     parser.add_argument('--mediadb', help="Media database file")
 
+    parser.add_argument('--spent-medium', metavar='SPECIES', dest='spent',
+                        help="Add spent medium compounds generated from given species (SBML model).")
+
     univ = parser.add_mutually_exclusive_group()
     univ.add_argument('-u', '--universe', dest='universe', help="Pre-built universe model (default: bacteria)")
     univ.add_argument('--universe-file', dest='universe_file', help="Reaction universe file (SBML format)")
@@ -115,4 +129,5 @@ if __name__ == '__main__':
          outputfile=args.output,
          flavor=flavor,
          exchange_format=args.exchange_format,
+         spent=args.spent,
          verbose=args.verbose)