Update expand.py #62
Update expand.py #62
Conversation
The new medusa with continuous gap-filling based on metabolic probability and ensembles of gap-filled metabolic models generating with some noise to the penality.
There was a problem hiding this comment.
Hi Basazin, I've highlighted some changes that are needed in this code review. The biggest change we need is to transfer your new functionality away from the iterative gapfilling function, which is only intended for scenarios where the user has multiple media conditions they are using for gapfilling (such as growth data from Biolog phenotype arrays). Your functionality should be in gapfill_to_ensemble. Please try to transfer over your code and leave the iterative gapfilling functions as they were originally. After that we can fine tune your implementation!
| """ | ||
| Performs gapfilling on model, pulling reactions from universal. | ||
| Any existing constraints on base_model are maintained during gapfilling, so | ||
| these should be set before calling gapfill_to_ensemble (e.g. secretion of | ||
| metabolites, choice of objective function etc.). | ||
|
|
||
| Currently, only iterative solutions are supported with accumulating | ||
| penalties (i.e. after each iteration, the penalty for each reaction | ||
| doubles). | ||
|
|
||
| Parameters | ||
| ---------- | ||
| model : cobra.Model | ||
| The model to perform gap filling on. | ||
| universal : cobra.Model | ||
| A universal model with reactions that can be used to complete the | ||
| model. | ||
| lower_bound : float, 0.05 | ||
| The minimally accepted flux for the objective in the filled model. | ||
| penalties : dict, None | ||
| A dictionary with keys being 'universal' (all reactions included in | ||
| the universal model), 'exchange' and 'demand' (all additionally | ||
| added exchange and demand reactions) for the three reaction types. | ||
| Can also have reaction identifiers for reaction specific costs. | ||
| Defaults are 1, 100 and 1 respectively. | ||
| integer_threshold : float, 1e-6 | ||
| The threshold at which a value is considered non-zero (aka | ||
| integrality threshold). If gapfilled models fail to validate, | ||
| you may want to lower this value. However, picking a threshold that is | ||
| too low may also result in reactions being added that are not essential | ||
| to meet the imposed constraints. | ||
| exchange_reactions : bool, False | ||
| Consider adding exchange (uptake) reactions for all metabolites | ||
| in the model. | ||
| demand_reactions : bool, False | ||
| Consider adding demand reactions for all metabolites. | ||
|
|
||
| Returns | ||
| ------- | ||
| ensemble : medusa.core.Ensemble | ||
| The ensemble object created from the gapfill solutions. | ||
| """ |
There was a problem hiding this comment.
Make sure you don't delete these docstrings for any functions--they are automatically parsed to build the documentation at https://medusa.readthedocs.io/en/latest/
| return ensemble | ||
|
|
||
| def iterative_gapfill_from_binary_phenotypes(model,universal,phenotype_dict, |
There was a problem hiding this comment.
We still need the iterative_gapfill_from_binary_phenotypes function--please add this back in your pull request
| return ensemble | ||
|
|
||
| def _continuous_iterative_binary_gapfill(model,phenotype_dict,cycle_order, | ||
| def _continuous_iterative_binary_gapfill(model = None, |
There was a problem hiding this comment.
The probabilistic gap filling functionality should be added to the gapfill_to_ensemble function instead of _continuous_iterative_binary_gapfill (which is for gapfilling when multiple media conditions are available for gap filling, such as in Biolog growth assays).
medusa/reconstruct/expand.py
Outdated
|
|
||
| ## set a constraint on flux through the original objective | ||
| for reaction in original_objective.keys(): | ||
| print("Constraining lower bound for " + reaction) | ||
| gapfiller.reactions.get_by_id(reaction).lower_bound = lower_bound | ||
|
|
||
| exchange_reactions = [rxn for rxn in gapfiller.reactions if\ | ||
| rxn.id.startswith(exchange_prefix)] | ||
| for rxn in exchange_reactions: | ||
| rxn.lower_bound = 0 | ||
|
|
||
| for cycle_num in range(0,output_ensemble_size): | ||
| print("starting cycle number " + str(cycle_num)) |
There was a problem hiding this comment.
We still need this for the original function to work--make sure you add this back in when you transfer your new code over to gapfill_to_ensemble
medusa/reconstruct/expand.py
Outdated
| original_coefficients = \ | ||
| gapfiller.objective.get_linear_coefficients(gapfiller.variables) | ||
|
|
||
| for condition in cycle_order[cycle_num]: | ||
| # set the medium for this condition. | ||
| for ex_rxn in phenotype_dict[condition].keys(): | ||
| gapfiller.reactions.get_by_id(ex_rxn).lower_bound = \ | ||
| -1.0*phenotype_dict[condition][ex_rxn] | ||
| gapfiller.reactions.get_by_id(ex_rxn).upper_bound = \ | ||
| 1.0*phenotype_dict[condition][ex_rxn] |
There was a problem hiding this comment.
We still need this for the original function to work--make sure you add this back in when you transfer your new code over to gapfill_to_ensemble
medusa/reconstruct/expand.py
Outdated
|
|
||
| # Generate the coefficients list and Randomly varying the coefficients | ||
| iou = [] | ||
| for y in coefficients.values(): | ||
| iou.append(y) | ||
| solutions = [] | ||
| t = 0.05 | ||
| for n in range(100): | ||
| new_coeff = iou + np.random.normal(0, t, len(iou)) | ||
| coefficients = dict(zip(coefficients.keys(), abs(new_coeff))) | ||
|
|
||
| gapfiller.objective.set_linear_coefficients(coefficients) | ||
| for reaction in original_objective.keys(): | ||
| print("Constraining lower bound for " + reaction) | ||
| gapfiller.reactions.get_by_id(reaction).lower_bound = lower_bound |
There was a problem hiding this comment.
We'll make some adjustments to this code later, but get it transferred over to gapfill_to_ensemble and we will fine-tune things then
medusa/reconstruct/expand.py
Outdated
| # save_directory = ("/Users/basazinbelhu/probannopy/iterative_pfba/data/medusa_gapfilled_model_pickle/" +str(ensemble.name)+".pkl") | ||
| # ensemble.to_pickle(save_directory) |
There was a problem hiding this comment.
We shouldn't save the ensemble, just return it. If a user needs to save the ensemble, they can do so using the returned object in their own script.
medusa/reconstruct/expand.py
Outdated
| def _build_ensemble_from_gapfill_solutions(model,solutions,universal=None): | ||
|
|
||
|
|
||
| model.reactions.EX_cpd00007_e.lower_bound = 0 # out of the function on medusa |
There was a problem hiding this comment.
Make sure you aren't setting anaerobic conditions in the Medusa code--this should be up to the user to do outside of Medusa.
Updated Medusa based on Greg's feedback!! The probabilistic gap-filling moved to the gapfill_to_ensemble function. Some code added back to the _continuous_iterative_binary_gapfill function
This is the code that added to generate a metabolite probability-based gap-filling and generate ensembles of models based on some noise to the penalty
The new medusa with continuous gap-filling based on metabolic probability and ensembles of gap-filled metabolic models generating with some noise to the penalty.