Skip to content

Commit

Permalink
Fix a bug where get_loops interferes with loopless FVA, and where tra…
Browse files Browse the repository at this point in the history
…nsport reactions are not opened when preparing models for merging
  • Loading branch information
mpredl committed Jul 8, 2024
1 parent 1059c44 commit 1736372
Show file tree
Hide file tree
Showing 3 changed files with 909 additions and 415 deletions.
14 changes: 8 additions & 6 deletions src/pycomo/helper/multiprocess.py
Original file line number Diff line number Diff line change
Expand Up @@ -127,26 +127,28 @@ def _loopless_fva_step(rxn_id):
solution = _model.optimize("minimize")
if not solution.status == "infeasible":
if perform_loopless_on_rxn:
logger.debug("Starting loop correction")
logger.debug(f"{rxn.id} Starting loop correction")
with _model:
_add_loopless_constraints_and_objective(solution.fluxes)
logger.debug("Optimize for loopless flux")
logger.debug(f"{rxn.id} Optimize for loopless flux")
solution = _model.optimize()
logger.debug("Optimization for loopless flux finished")
logger.debug(f"{rxn.id} Optimization for loopless flux finished")
min_flux = solution.fluxes[rxn.id] if not solution.status == "infeasible" else 0.
else:
logger.debug(f"{rxn.id} min flux is infeasible")
min_flux = 0.
solution = _model.optimize("maximize")
if not solution.status == "infeasible":
if perform_loopless_on_rxn:
logger.debug("Starting loop correction")
logger.debug(f"{rxn.id} Starting loop correction")
with _model:
_add_loopless_constraints_and_objective(solution.fluxes)
logger.debug("Optimize for loopless flux")
logger.debug(f"{rxn.id} Optimize for loopless flux")
solution = _model.optimize()
logger.debug("Optimization for loopless flux finished")
logger.debug(f"{rxn.id} Optimization for loopless flux finished")
max_flux = solution.fluxes[rxn.id] if not solution.status == "infeasible" else 0.
else:
logger.debug(f"{rxn.id} max flux is infeasible")
max_flux = 0.
return rxn_id, max_flux, min_flux

Expand Down
30 changes: 19 additions & 11 deletions src/pycomo/pycomo_models.py
Original file line number Diff line number Diff line change
Expand Up @@ -304,7 +304,7 @@ def find_metabolites_without_exchange_rxn(self, model, exchg_comp=""):
mets_without_exchg = list((set(comp_mets) - set(exchg_mets)) - {self.biomass_met})
return mets_without_exchg

def convert_exchange_to_transport_reaction(self, model, old_comp, inplace=True):
def convert_exchange_to_transport_reaction(self, model, old_comp, inplace=True, max_flux=1000.):
"""
When preprocessing the model for merging into a community metabolic model, a new shared compartment is added.
This new compartment will be the compartment for boundary metabolites, containing all metabolites with exchange
Expand All @@ -314,6 +314,7 @@ def convert_exchange_to_transport_reaction(self, model, old_comp, inplace=True):
:param model: The model to be changed
:param old_comp: The ID of the old exchange compartment
:param inplace: If True, the input model will be changed, else a copy is made
:param max_flux: The maximum allowed flux in the model (used to open the constraints of the transport reactions)
:return: None if the inplace flag is set, otherwise the model with changed reactions
"""
if not inplace:
Expand All @@ -335,6 +336,7 @@ def convert_exchange_to_transport_reaction(self, model, old_comp, inplace=True):
rxn.add_metabolites(new_met_stoich)
if "Exchange" in rxn.name:
rxn.name = rxn.name.replace("Exchange", "Transport")
rxn.bounds = (-max_flux, max_flux)

model.repair()

Expand Down Expand Up @@ -423,7 +425,8 @@ def ensure_compartment_suffix(self, model, inplace=True):

return model

def add_exchange_compartment(self, model, exchg_comp_name=None, add_missing_transports=False, inplace=True):
def add_exchange_compartment(self, model, exchg_comp_name=None, add_missing_transports=False, inplace=True,
max_flux=1000.):
"""
Add a new exchange compartment. This will also copy all metabolites in the current external compartment to the
exchange compartment, establish transport reactions between the two compartments for all metabolites and add
Expand All @@ -434,6 +437,7 @@ def add_exchange_compartment(self, model, exchg_comp_name=None, add_missing_tran
:param add_missing_transports: If True, add exchange reactions and transport reactions for all metabolites in
the old external compartment that did not have any exchange reactions
:param inplace: If True, the input model will be changed, else a copy is made
:param max_flux: Maximum allowed flux
:return: Updated model
"""
if exchg_comp_name is None:
Expand All @@ -448,7 +452,7 @@ def add_exchange_compartment(self, model, exchg_comp_name=None, add_missing_tran
if add_missing_transports:
mets_without_exchg = self.find_metabolites_without_exchange_rxn(model)
self.add_exchange_reactions_to_metabolites(model, mets_without_exchg, inplace=True)
self.convert_exchange_to_transport_reaction(model, old_exc_comp, inplace=True)
self.convert_exchange_to_transport_reaction(model, old_exc_comp, inplace=True, max_flux=max_flux)
# Add exchange reactions
self.add_exchange_reactions_to_compartment(model, exchg_comp_name, inplace=True)

Expand Down Expand Up @@ -541,7 +545,7 @@ def remove_biomass_exchange_rxn(self, model, remove_all_consuming_rxns=True):
reaction.remove_from_model(remove_orphans=True)
return

def prepare_for_merging(self, shared_compartment_name=None):
def prepare_for_merging(self, shared_compartment_name=None, max_flux=1000.):
"""
Prepares the model for merging into a community metabolic model. The generated model is SBML conform and all
genes, reactions, compartments and metabolites contain the information that they belong to this model (which is
Expand All @@ -554,6 +558,7 @@ def prepare_for_merging(self, shared_compartment_name=None):
- Prefix all compartments, genes, reactions and metabolites with the name of the model
:param shared_compartment_name: The ID of the new shared exchange compartment
:param max_flux: Maximum allowed flux
:return: The model prepared for merging
"""
if shared_compartment_name is not None:
Expand All @@ -573,7 +578,7 @@ def prepare_for_merging(self, shared_compartment_name=None):
for comp in self.prepared_model.compartments:
rename[comp] = self.name + "_" + comp
self.rename_compartment(self.prepared_model, rename)
self.add_exchange_compartment(self.prepared_model, add_missing_transports=True)
self.add_exchange_compartment(self.prepared_model, add_missing_transports=True, max_flux=max_flux)
self.prefix_metabolite_names(self.prepared_model, self.name + "_",
exclude_compartment=self.shared_compartment_name)
self.prefix_reaction_names(self.prepared_model, self.name + "_",
Expand Down Expand Up @@ -897,7 +902,8 @@ def get_loops(self, processes=None):
no_medium = {}
self.model.medium = no_medium

solution_df = find_loops_in_model(self.convert_to_model_without_fraction_metabolites(), processes=processes)
with self.model:
solution_df = find_loops_in_model(self.convert_to_model_without_fraction_metabolites(), processes=processes)

self.medium = original_medium
self.apply_medium()
Expand Down Expand Up @@ -981,7 +987,8 @@ def generate_community_model(self):
for model in self.member_models:
idx += 1
if idx == 1:
merged_model = model.prepare_for_merging(shared_compartment_name=self.shared_compartment_name)
merged_model = model.prepare_for_merging(shared_compartment_name=self.shared_compartment_name,
max_flux=self.max_flux)
biomass_met_id = model.biomass_met.id
biomass_met = merged_model.metabolites.get_by_id(biomass_met_id)
biomass_mets[model.name] = biomass_met
Expand All @@ -992,7 +999,8 @@ def generate_community_model(self):
self.create_fraction_reaction(merged_model, member_name=model.name)

else:
extended_model = model.prepare_for_merging(shared_compartment_name=self.shared_compartment_name)
extended_model = model.prepare_for_merging(shared_compartment_name=self.shared_compartment_name,
max_flux=self.max_flux)

unbalanced_metabolites = check_mass_balance_of_metabolites_with_identical_id(extended_model,
merged_model)
Expand Down Expand Up @@ -1585,7 +1593,7 @@ def loopless_fva(self, reaction_ids, fraction_of_optimum=None, use_loop_reaction
use_loop_reactions_for_ko=use_loop_reactions_for_ko,
ko_candidate_ids=ko_candidate_ids,
verbose=verbose,
processes=None)
processes=processes)

def run_fva(self, fraction_of_optimum=0.9, composition_agnostic=False, loopless=False, fva_mu_c=None,
only_exchange_reactions=True, reactions=None, verbose=False, processes=None):
Expand Down Expand Up @@ -1659,7 +1667,7 @@ def run_fva(self, fraction_of_optimum=0.9, composition_agnostic=False, loopless=
self.change_reaction_bounds("community_biomass", lower_bound=0., upper_bound=0.)
for member_name, fraction_met in f_bio_mets.items():
biomass_rxn = model.reactions.get_by_id(f"{member_name}_to_community_biomass")
biomass_rxn.add_metabolites({fraction_met: -1}, combine=False)
biomass_rxn.add_metabolites({fraction_met: -1}, combine=True)

self.apply_fixed_growth_rate(mu_c)
else:
Expand Down Expand Up @@ -1698,7 +1706,7 @@ def run_fva(self, fraction_of_optimum=0.9, composition_agnostic=False, loopless=
self.change_reaction_bounds("community_biomass", lower_bound=0., upper_bound=0.)
for member_name, fraction_met in f_bio_mets.items():
biomass_rxn = model.reactions.get_by_id(f"{member_name}_to_community_biomass")
biomass_rxn.add_metabolites({fraction_met: -1}, combine=False)
biomass_rxn.add_metabolites({fraction_met: -1}, combine=True)

self.convert_to_fixed_abundance()
else:
Expand Down
Loading

0 comments on commit 1736372

Please sign in to comment.