implemented pFBA and FVA

xbanalysis/model/xba_compartment.py

@@ -3,11 +3,15 @@
 Peter Schubert, HHU Duesseldorf, May 2022
 """
 
+
 class XbaCompartment:
 
     def __init__(self, s_compartment):
         self.id = s_compartment.name
-        self.name = s_compartment.get('name', '')
+        self.name = s_compartment.get('name', self.id)
+        if 'size' in s_compartment:
             self.size = s_compartment['size']
+        if 'units' in s_compartment:
             self.units = s_compartment['units']
+        if 'spatialDimension' in s_compartment:
             self.dimension = s_compartment['spatialDimension']
\ No newline at end of file

xbanalysis/model/xba_function.py

@@ -11,7 +11,7 @@ class XbaFunction:
 
     def __init__(self, s_function):
         self.id = s_function.name
-        self.name = s_function.get('name', '')
+        self.name = s_function.get('name', self.id)
         self.math = s_function['math']
 
         m = re.search(r'\blambda\((.*)\)', sbmlxdf.misc.mathml2numpy(s_function['math']))

xbanalysis/model/xba_model.py

@@ -36,6 +36,9 @@ class XbaModel:
         if 'funcDefs' in model_dict:
             self.functions = {fid: XbaFunction(row)
                               for fid, row in model_dict['funcDefs'].iterrows()}
+        else:
+            self.functions = {}
+
         self.species = {sid: XbaSpecies(row)
                         for sid, row in model_dict['species'].iterrows()}
         self.reactions = {rid: XbaReaction(row, self.functions, self.compartments)

xbanalysis/model/xba_parameter.py

@@ -8,8 +8,9 @@ class XbaParameter:
 
     def __init__(self, s_parameter):
         self.id = s_parameter.name
-        self.name = s_parameter.get('name', '')
-        self.sboterm = s_parameter.get('sboterm', '')
+        self.name = s_parameter.get('name', self.id)
         self.value = s_parameter['value']
         self.constant = s_parameter['constant']
         self.units = s_parameter['units']
+        if 'sboterm' in s_parameter:
+            self.sboterm = s_parameter['sboterm']
\ No newline at end of file

xbanalysis/model/xba_reaction.py

@@ -12,16 +12,22 @@ class XbaReaction:
 
     def __init__(self, s_reaction, functions, compartments):
         self.id = s_reaction.name
-        self.name = s_reaction.get('name', '')
+        self.name = s_reaction.get('name', self.id)
+        if 'sboterm' in s_reaction:
             self.sboterm = s_reaction['sboterm']
         self.reaction_string = s_reaction['reactionString']
         self.reversible = s_reaction['reversible']
         self.products = get_species_stoic(s_reaction['products'])
         self.reactants = get_species_stoic(s_reaction['reactants'])
+        if 'modifiers' in s_reaction:
             self.modifiers = get_species_stoic(s_reaction['modifiers'])
-        self.fbc_lb = s_reaction.get('fbcLb', float('nan'))
-        self.fbc_ub = s_reaction.get('fbcUb', float('nan'))
+        else:
+            self.modifiers = {}
+        if 'fbcLb' in s_reaction and 'fbcUb' in s_reaction:
+            self.fbc_lb = s_reaction['fbcLb']
+            self.fbc_ub = s_reaction['fbcUb']
         self.local_params = get_local_params(s_reaction)
+        if 'kineticLaw' in s_reaction:
             self.kinetic_law = s_reaction['kineticLaw'] if type(s_reaction['kineticLaw']) == str else ''
             # convert to numpy, replace local parameters with numerical values, inline lambda functions:
             expanded_kl = sbmlxdf.misc.mathml2numpy(self.kinetic_law)

xbanalysis/model/xba_species.py

@@ -14,13 +14,17 @@ class XbaSpecies:
 
     def __init__(self, s_species):
         self.id = s_species.name
-        self.name = s_species.get('name', '')
-        self.sboterm = s_species['sboterm']
+        self.name = s_species.get('name', self.id)
+        # set sboterm to 000247 to support Flux Balance Analysis
+        self.sboterm = s_species.get('sboterm', 'SBO:0000247')
         self.compartment = s_species['compartment']
-        self.initial_conc = s_species['initialConcentration']
         self.constant = s_species['constant']
         self.boundary = s_species['boundaryCondition']
+        if 'initialConcentration' in s_species:
+            self.initial_conc = s_species['initialConcentration']
+        if 'units' in s_species:
             self.units = s_species['substanceUnits']
+        if 'xmlAnnotation' in s_species:
             attrs = sbmlxdf.extract_xml_attrs(s_species['xmlAnnotation'], ns=xml_gba_ns, token=xml_token)
             self.mw = float(attrs['weight_Da'])
         self.ps_rid = None

xbanalysis/problems/fba_problem.py

@@ -4,6 +4,7 @@ Peter Schubert, HHU Duesseldorf, June 2022
 """
 
 import numpy as np
+import pandas as pd
 
 from xbanalysis.problems.lp_problem import LpProblem
 
@@ -39,6 +40,7 @@ class FbaProblem:
 
         self.metab_sids = [s.id for s in xba_model.species.values() if s.sboterm == 'SBO:0000247']
         metabs = set(self.metab_sids)
+        # identify reactions that do not use enzymes in reactants and products
         self.mrids = [rid for rid, rxn in xba_model.reactions.items()
                       if len((set(rxn.products) | set(rxn.reactants)) - metabs) == 0]
         self.mrid2idx = {mrid: idx for idx, mrid in enumerate(self.mrids)}
@@ -56,58 +58,125 @@ class FbaProblem:
                 self.set_obj_coefs(fbc_obj.coefficiants)
                 break
 
+    def create_fba_lp(self):
+        """create Linear Programming Problem for FBA with split reactions
+
+        Splitting reactions in forward/backward directions reduces high values during FBA
+        Splitting required for pFBA (minimize absolute flux values)
+
+        :return: LpProblem configured for FBA
+        :rtype: LpProblem
+        """
+        lp = LpProblem()
+        # split forward and backward reactions
+        fwd_bwd_s_mat = np.hstack((self.s_mat, -self.s_mat))
+        fwd_bwd_obj_coefs = np.hstack((self.obj_coefs, -self.obj_coefs))
+        fwd_bwd_bounds = np.hstack((np.fmax(self.flux_bounds, 0),
+                                    np.vstack((np.fmax(-self.flux_bounds[1], 0),
+                                               np.fmax(-self.flux_bounds[0], 0)))))
+        lp.configure(fwd_bwd_obj_coefs, fwd_bwd_bounds,
+                     fwd_bwd_s_mat, np.zeros((2, fwd_bwd_s_mat.shape[0])), self.obj_dir)
+        return lp
+
     def fba_optimize(self):
         """FBA optimization of current problem
 
-        Instantiate an LpProblem, configure it and delete it again.
-
         :return: results information from the optimization
         :rtype: dict
         """
-        lp = LpProblem()
-        lp.configure(self.obj_coefs, self.flux_bounds,
-                     self.s_mat, np.zeros((2, self.s_mat.shape[0])), self.obj_dir)
+        lp = self.create_fba_lp()
         res = lp.solve()
         del lp
+
+        # combine results for split fluxes
+        n_cols = len(self.mrids)
+        res['x'] = res['x'][:n_cols] - res['x'][n_cols:]
+        res['reduced_costs'] = res['reduced_costs'][:n_cols] - res['reduced_costs'][n_cols:]
         return res
 
     def pfba_optimize(self, fract_of_optimum=1.0):
         """ pfba analysis
 
-        As we are interested in minimizing the absolute value of fluxes,
-        we split fluxes in forward and reverse fluxes to make them all positive
+        Reactions need to be split in forward/backard (to minimize absolute flux value)
 
         :param fract_of_optimum: value between 0.0 and 1.0.
         :type fract_of_optimum: float (default 1.0)
         :return:
         """
-        # create and solve the FBA problem with the configured objective function
-        # for pFBA we need to split fluxes in forward and backward fluxes, with bounds set as such
-        #  that fluxes can only be positive, lb>=0, ub>=0
-
-        n_cols = len(self.mrids)
-        lp = LpProblem()
-
-        fwd_bwd_s_mat = np.hstack((self.s_mat, -self.s_mat))
-        fwd_bwd_obj_coefs = np.hstack((self.obj_coefs, -self.obj_coefs))
-        fwd_bwd_bounds = np.hstack((np.fmax(self.flux_bounds, 0),
-                                    np.vstack((np.fmax(-self.flux_bounds[1], 0),
-                                               np.fmax(-self.flux_bounds[0], 0)))))
-        lp.configure(fwd_bwd_obj_coefs, fwd_bwd_bounds,
-                     fwd_bwd_s_mat, np.zeros((2, fwd_bwd_s_mat.shape[0])), self.obj_dir)
+        lp = self.create_fba_lp()
         res = lp.solve(short_result=True)
-        # fix FBA objective as constraint with lower bound and minimize total fluxes
         if res['simplex_status'] == 0:
-            row_bds = np.array([fract_of_optimum * res['fun'], np.nan])
+            # fix FBA objective as constraint
+            fwd_bwd_obj_coefs = np.hstack((self.obj_coefs, -self.obj_coefs))
+            if self.obj_dir == 'maximize':
+                row_bds = np.array([res['fun'] * fract_of_optimum, np.nan])
+            else:
+                row_bds = np.array([np.nan, res['fun'] / fract_of_optimum])
             lp.add_constraint(fwd_bwd_obj_coefs, row_bds)
+            # new LP objective is to minimize sum of fluxes
+            n_cols = len(self.mrids)
             lp.set_obj_coefs(np.ones(2 * n_cols))
             lp.set_obj_dir('minimize')
             res = lp.solve()
-        del lp
+            # combine results for split fluxes
             res['x'] = res['x'][:n_cols] - res['x'][n_cols:]
             res['reduced_costs'] = res['reduced_costs'][:n_cols] - res['reduced_costs'][n_cols:]
+        del lp
         return res
 
+    def fva_optimize(self, rids=None, fract_of_optimum=1.0):
+        """FVA - flux variability analysis for all or selected reactions.
+
+        :param rids: selected reaction ids to run FVA (alternatively, all)
+        :type rids: list of strings (default: None - all reactions)
+        :param fract_of_optimum: value between 0.0 and 1.0.
+        :type fract_of_optimum: float (default 1.0)
+        :return: Flux ranges and objective values
+        :rtype: pandas DataFrame
+        """
+        lp = self.create_fba_lp()
+        res = lp.solve(short_result=True)
+
+        if rids is None:
+            rids = self.mrids
+        flux_ranges = np.zeros((len(rids), 4))
+        if res['simplex_status'] == 0:
+            # fix FBA objective as constraint
+            fwd_bwd_obj_coefs = np.hstack((self.obj_coefs, -self.obj_coefs))
+            if self.obj_dir == 'maximize':
+                row_bds = np.array([res['fun'] * fract_of_optimum, np.nan])
+            else:
+                row_bds = np.array([np.nan, res['fun'] / fract_of_optimum])
+            fba_constr_row = lp.add_constraint(fwd_bwd_obj_coefs, row_bds)
+
+            n_cols = len(self.mrids)
+            for idx, rid in enumerate(rids):
+                # select reaction maximize/minimize flux
+                obj_coefs = np.zeros(2 * n_cols)
+                obj_coefs[self.mrid2idx[rid]] = 1
+                obj_coefs[self.mrid2idx[rid] + n_cols] = -1
+                lp.set_obj_coefs(obj_coefs)
+
+                lp.set_obj_dir('minimize')
+                res = lp.solve(short_result=True)
+                if res['simplex_status'] == 0:
+                    flux_ranges[idx, 0] = res['fun']
+                    flux_ranges[idx, 2] = lp.get_row_value(fba_constr_row)
+                else:
+                    print(f"FVA min failed for {rid}, {res['simplex_message']}")
+                    flux_ranges[idx, 0] = np.nan
+
+                lp.set_obj_dir('maximize')
+                res = lp.solve(short_result=True)
+                if res['simplex_status'] == 0:
+                    flux_ranges[idx, 1] = res['fun']
+                    flux_ranges[idx, 3] = lp.get_row_value(fba_constr_row)
+                else:
+                    print(f"FVA max failed for {rid}, {res['simplex_message']}")
+                    flux_ranges[idx, 1] = np.nan
+        del lp
+        return pd.DataFrame(flux_ranges, index=rids, columns=['min', 'max', 'obj_min', 'obj_max'])
+
     @property
     def obj_dir(self):
         return self.__obj_dir

xbanalysis/problems/lp_problem.py

@@ -255,6 +255,13 @@ class LpProblem:
         return self.results(short_result)
 
     def get_row_value(self, row):
+        """Retrieve right hand side of a constraint (auxiliary variable)
+
+        :param row: index into constraints
+        :type row: int
+        :return: value of auxiliary varibale (value of rhs)
+        :rtype: float
+        """
         assert row <= glpk.glp_get_num_rows(self.lp)
         return glpk.glp_get_row_prim(self.lp, row)
 

xbanalysis/utils/utils.py

@@ -93,7 +93,6 @@ def get_local_params(reaction):
     return local_params
 
 
-# added
 def get_species_stoic(srefs):
     """Retrieve species and stoichiometry from reaction srefs.