Commit dafaae94 authored by Peter Schubert's avatar Peter Schubert
Browse files

added support function get_obj_values()

parent 95dfe2c2
"""
Definition of version string.
"""
__version__ = "0.1.0"
__version__ = "0.1.2"
program_name = 'networkred'
......@@ -354,6 +354,35 @@ class ReduceModel:
self._active_constraints = None
logging.info('constraints removed')
def get_opt_flux_values(self):
"""Get optimal flux values for constraints in protected functions.
Useful for determining upper/lower bounds in protected functions.
:returns: copy of protected functions dataFrame with optimal
fluxes added
:rtype: pandas DataFrame
"""
if self._protect_funcs is not None:
df_fluxes = self._protect_funcs.copy()
df_fluxes.reset_index(inplace=True)
df_fluxes['opt_flux'] = np.nan
for k in set(self._protect_funcs.index):
self._add_constraints(k)
solution = self._red_model.optimize()
self._remove_constraints()
for idx, row in df_fluxes[df_fluxes['k'] == k].iterrows():
flux = 0.0
for item in get_items(row['expression']):
params = extract_params(item)
flux += solution.fluxes[params['rid']] * float(params['coeff'])
if solution.status == 'optimal':
df_fluxes.at[idx, 'opt_flux'] = flux
return df_fluxes
else:
return None
def check_feasibility(self, feas=True):
""" Check feasibility of protected functions, reactions and metabolites.
......@@ -446,43 +475,6 @@ class ReduceModel:
'essential_r': len(self._essential_rids),
}
def _remove_reactions(self, remove_rids):
"""Remove reactions from the model, including orphans.
also delete the corresponding group associations
:param remove_rids: reaction IDs to be removed
:type remove_rids: set or list of str
"""
# first, remove any group association
for rid in remove_rids:
for group in self._red_model.get_associated_groups(rid):
group.remove_members([rid])
self._red_model.remove_reactions(list(remove_rids), remove_orphans=True)
def remove_blocked(self):
"""Remove blocked reactions from the model
:returns: status information in a dictionary
:rtype: dict
"""
# identify blocked reactions:
blocked_rids = set(find_blocked_reactions(self._red_model, processes=self._cpus))
logging.info('number of blocked reactions identified: %d', len(blocked_rids))
protected_blocked_rids = blocked_rids.intersection(self._protect_rids)
if len(protected_blocked_rids) > 0:
logging.info('blocked reaction being protected: %s', ', '.join(protected_blocked_rids))
# remove blocked reactions:
remove_rids = blocked_rids.difference(protected_blocked_rids)
if len(remove_rids) > 0:
self._remove_reactions(remove_rids)
logging.info('blocked reactions removed: %s', ', '.join(remove_rids))
remain_mids = set([m.id for m in self._red_model.metabolites])
protected_removed_mids = self._protect_mids.difference(remain_mids)
if len(protected_removed_mids) > 0:
logging.warning('protected metabolites removed: %s', ', '.join(protected_removed_mids))
def set_solver_params(self, params):
"""Set Solver type and Solver Parameters
......@@ -520,7 +512,44 @@ class ReduceModel:
'timeout': self._red_model.solver.configuration.timeout
}
def remove_parallel_rids(self):
def _remove_reactions(self, remove_rids):
"""Remove reactions from the model, including orphans.
also delete the corresponding group associations
:param remove_rids: reaction IDs to be removed
:type remove_rids: set or list of str
"""
# first, remove any group association
for rid in remove_rids:
for group in self._red_model.get_associated_groups(rid):
group.remove_members([rid])
self._red_model.remove_reactions(list(remove_rids), remove_orphans=True)
def remove_blocked_reactions(self):
"""Remove blocked reactions from the model
:returns: status information in a dictionary
:rtype: dict
"""
# identify blocked reactions:
blocked_rids = set(find_blocked_reactions(self._red_model, processes=self._cpus))
logging.info('number of blocked reactions identified: %d', len(blocked_rids))
protected_blocked_rids = blocked_rids.intersection(self._protect_rids)
if len(protected_blocked_rids) > 0:
logging.info('blocked reaction being protected: %s', ', '.join(protected_blocked_rids))
# remove blocked reactions:
remove_rids = blocked_rids.difference(protected_blocked_rids)
if len(remove_rids) > 0:
self._remove_reactions(remove_rids)
logging.info('blocked reactions removed: %s', ', '.join(remove_rids))
remain_mids = set([m.id for m in self._red_model.metabolites])
protected_removed_mids = self._protect_mids.difference(remain_mids)
if len(protected_removed_mids) > 0:
logging.warning('protected metabolites removed: %s', ', '.join(protected_removed_mids))
def remove_parallel_reactions(self):
"""Identify and remove parallel reactions.
"""
......@@ -631,7 +660,6 @@ class ReduceModel:
def print_status(self):
""" Print status along reduction process
"""
self.remove_blocked()
stat = self.get_reduction_status()
print('reactions removed/remaining/protected: %4d/%4d/%4d; dof: %3d' %
(stat['delta_r'], stat['remain_r'], stat['essential_r'], stat['dof']))
......@@ -654,11 +682,11 @@ class ReduceModel:
self.print_status()
print('remove blocked reactions')
self.remove_blocked()
self.remove_blocked_reactions()
self.print_status()
print('remove parallel reactions')
self.remove_parallel_rids()
self.remove_parallel_reactions()
self.print_status()
if save_pts > 0:
......@@ -737,7 +765,7 @@ class ReduceModel:
os.remove(snapshot)
logging.info('removed snapshot')
def _add_orig_data(self, sbml_fname):
def _add_original_data(self, sbml_fname):
"""Adds units from original SBML file to reduced model"
NOTE: at the moment this only works properly for
......@@ -803,7 +831,7 @@ class ReduceModel:
cobra.io.write_sbml_model(self._red_model, sbml_fname)
if type(self._orig_model_ref) == str:
self._add_orig_data(sbml_fname)
self._add_original_data(sbml_fname)
pass
def get_cobra_model(self):
......
<?xml version="1.0" encoding="UTF-8"?>
<!-- Created by sbmlxdf version 0.2.5 on 2021-11-05 09:00 with libSBML version 5.19.0. -->
<!-- Created by sbmlxdf version 0.2.5 on 2021-11-05 11:24 with libSBML version 5.19.0. -->
<sbml xmlns="http://www.sbml.org/sbml/level3/version2/core" xmlns:fbc="http://www.sbml.org/sbml/level3/version1/fbc/version2" level="3" version="2" fbc:required="false">
<model metaid="Deniz_model_fba" id="Deniz_model_fba_reduced" name="Deniz_model_fba" substanceUnits="substance" timeUnits="time" extentUnits="substance" fbc:strict="true">
<notes>
......
......@@ -31,7 +31,7 @@ def reduce_model():
sbml_dir = 'sample_data/SBML_models'
data_dir = 'sample_data/data'
model_name = 'Deniz_model_fba'
protected_parts_file = 'Deniz_model_fba_pf.xlsx'
protected_parts_file = 'Deniz_model_fba_nrp.xlsx'
# load the original model
sbml_file = os.path.join(sbml_dir, model_name) + '.xml'
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment