This is an automated email from the git hooks/post-receive script. afif pushed a commit to branch master in repository python-cobra.
commit 544431853b10e6b06b098523786ac3f681a9df9d Author: Afif Elghraoui <[email protected]> Date: Tue Oct 17 22:11:35 2017 -0400 New upstream version 0.9.0 --- .travis.yml | 2 +- appveyor.yml | 2 +- cobra/__init__.py | 2 +- cobra/core/gene.py | 4 +- cobra/core/metabolite.py | 7 +- cobra/core/model.py | 65 ++++-- cobra/core/reaction.py | 43 ++-- cobra/core/solution.py | 11 +- cobra/flux_analysis/phenotype_phase_plane.py | 305 +++++++++++++++------------ cobra/flux_analysis/summary.py | 10 +- cobra/manipulation/validate.py | 4 - cobra/test/data/iJO1366.pickle | Bin 1810497 -> 1810487 bytes cobra/test/data/mini.pickle | Bin 32825 -> 33312 bytes cobra/test/data/raven.pickle | Bin 13758 -> 13310 bytes cobra/test/data/salmonella.pickle | Bin 2151196 -> 2151186 bytes cobra/test/data/textbook_solution.pickle | Bin 143584 -> 7330 bytes cobra/test/test_flux_analysis.py | 22 +- cobra/test/test_manipulation.py | 3 +- cobra/test/test_model.py | 27 +++ cobra/util/util.py | 7 + release-notes/0.8.2.md | 13 +- release-notes/0.9.0.md | 30 +++ release-notes/next-release.md | 2 + scripts/compare-benchmark.py | 16 +- setup.cfg | 2 +- setup.py | 2 +- 26 files changed, 365 insertions(+), 214 deletions(-) diff --git a/.travis.yml b/.travis.yml index db8a64a..5bb203c 100644 --- a/.travis.yml +++ b/.travis.yml @@ -68,7 +68,7 @@ matrix: before_install: - if [[ -n "${MB_PYTHON_VERSION}" ]]; then - (travis_retry git clone https://github.com/matthew-brett/multibuild.git && cd multibuild && git checkout edf5b691d0d565b4e65e655b983c11c883acbeca); + (travis_retry git clone https://github.com/matthew-brett/multibuild.git && cd multibuild && git checkout 37040e31b1044468027bd86991c805006a92bf09); TEST_DEPENDS="swiglpk optlang sympy decorator cython codecov coverage numpy scipy jsonschema six pytest pytest-cov pytest-benchmark tabulate"; BUILD_DEPENDS="swiglpk optlang sympy cython numpy scipy"; source multibuild/common_utils.sh; diff --git a/appveyor.yml b/appveyor.yml index 2571dc9..017a3a4 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -84,7 +84,7 @@ install: python -m pip install pip setuptools>=24.0 wheel --upgrade } - if not defined CONDA %WITH_COMPILER% python -m pip install --upgrade pytest - if defined CONDA conda install -q setuptools pytest numpy scipy - - python -m pip install pytest-cov pytest-benchmark pandas swiglpk optlang python-libsbml decorator Cython jsonschema twine pypandoc==1.1.3 + - python -m pip install pytest-cov pytest-benchmark pandas swiglpk optlang python-libsbml decorator Cython jsonschema twine pypandoc - "%WITH_COMPILER% python appveyor/build_glpk.py" build: off diff --git a/cobra/__init__.py b/cobra/__init__.py index 21666c1..ce83d7a 100644 --- a/cobra/__init__.py +++ b/cobra/__init__.py @@ -13,7 +13,7 @@ from cobra.core import ( DictList, Gene, Metabolite, Model, Object, Reaction, Species) from cobra.util.version_info import show_versions -__version__ = "0.8.2" +__version__ = "0.9.0" # set the warning format to be prettier and fit on one line _cobra_path = _dirname(_abspath(__file__)) diff --git a/cobra/core/gene.py b/cobra/core/gene.py index 4ea1f1b..b869b06 100644 --- a/cobra/core/gene.py +++ b/cobra/core/gene.py @@ -11,6 +11,7 @@ from warnings import warn from cobra.core.species import Species from cobra.util import resettable +from cobra.util.util import format_long_string keywords = list(kwlist) keywords.remove("and") @@ -289,4 +290,5 @@ class Gene(Species): functional=self.functional, address='0x0%x' % id(self), n_reactions=len(self.reactions), - reactions=', '.join(r.id for r in self.reactions)) + reactions=format_long_string( + ', '.join(r.id for r in self.reactions), 200)) diff --git a/cobra/core/metabolite.py b/cobra/core/metabolite.py index 8f78e39..c3cd0f7 100644 --- a/cobra/core/metabolite.py +++ b/cobra/core/metabolite.py @@ -12,6 +12,7 @@ from cobra.exceptions import OptimizationError from cobra.core.formula import elements_and_molecular_weights from cobra.core.species import Species from cobra.util.solver import check_solver_status +from cobra.util.util import format_long_string # Numbers are not required because of the |(?=[A-Z])? block. See the @@ -257,8 +258,10 @@ class Metabolite(Species): <td><strong>In {n_reactions} reaction(s)</strong></td><td> {reactions}</td> </tr> - </table>""".format(id=self.id, name=self.name, formula=self.formula, + </table>""".format(id=self.id, name=format_long_string(self.name), + formula=self.formula, address='0x0%x' % id(self), compartment=self.compartment, n_reactions=len(self.reactions), - reactions=', '.join(r.id for r in self.reactions)) + reactions=format_long_string( + ', '.join(r.id for r in self.reactions), 200)) diff --git a/cobra/core/model.py b/cobra/core/model.py index fb0be9a..7eabf78 100644 --- a/cobra/core/model.py +++ b/cobra/core/model.py @@ -24,7 +24,7 @@ from cobra.util.solver import ( SolverNotFound, get_solver_name, interface_to_str, set_objective, solvers, add_cons_vars_to_problem, remove_cons_vars_from_problem, choose_solver, check_solver_status, assert_optimal) -from cobra.util.util import AutoVivification +from cobra.util.util import AutoVivification, format_long_string LOGGER = logging.getLogger(__name__) @@ -63,8 +63,6 @@ class Model(Object): for y in ['reactions', 'genes', 'metabolites']: for x in getattr(self, y): x._model = self - if y == 'reactions': - x._reset_var_cache() if not hasattr(self, "name"): self.name = None @@ -94,7 +92,7 @@ class Model(Object): self.reactions = DictList() # A list of cobra.Reactions self.metabolites = DictList() # A list of cobra.Metabolites # genes based on their ids {Gene.id: Gene} - self.compartments = dict() + self._compartments = dict() self._contexts = [] # from cameo ... @@ -148,9 +146,6 @@ class Model(Object): # Do nothing if the solver did not change if self.problem == interface: return - - for reaction in self.reactions: - reaction._reset_var_cache() self._solver = interface.Model.clone(self._solver) @property @@ -165,10 +160,37 @@ class Model(Object): def get_metabolite_compartments(self): """Return all metabolites' compartments.""" + warn('use Model.compartments instead', DeprecationWarning) return {met.compartment for met in self.metabolites if met.compartment is not None} @property + def compartments(self): + return {met.compartment: self._compartments.get(met.compartment, '') + for met in self.metabolites if met.compartment is not None} + + @compartments.setter + def compartments(self, value): + """Get or set the dictionary of current compartment descriptions. + + Assigning a dictionary to this property updates the model's + dictionary of compartment descriptions with the new values. + + Parameters + ---------- + value : dict + Dictionary mapping compartments abbreviations to full names. + + Examples + -------- + >>> import cobra.test + >>> model = cobra.test.create_test_model("textbook") + >>> model.compartments = {'c': 'the cytosol'} + {'c': 'the cytosol', 'e': 'extracellular'} + """ + self._compartments.update(value) + + @property def medium(self): def is_active(reaction): @@ -300,9 +322,6 @@ class Model(Object): new_gene = new.genes.get_by_id(gene.id) new_reaction._genes.add(new_gene) new_gene._reaction.add(new_reaction) - - for reaction in new.reactions: - reaction._reset_var_cache() try: new._solver = deepcopy(self.solver) # Cplex has an issue with deep copies @@ -521,7 +540,6 @@ class Model(Object): # Add reactions. Also take care of genes and metabolites in the loop for reaction in reaction_list: - reaction._reset_var_cache() reaction._model = self # the reaction now points to the model # keys() is necessary because the dict will be modified during # the loop @@ -604,7 +622,6 @@ class Model(Object): reaction._model = None if context: - context(reaction._reset_var_cache) context(partial(setattr, reaction, '_model', self)) context(partial(self.reactions.add, reaction)) @@ -940,6 +957,28 @@ class Model(Object): value = {rxn: 1 for rxn in reactions} set_objective(self, value, additive=False) + @property + def objective_direction(self): + """ + Get or set the objective direction. + + When using a `HistoryManager` context, this attribute can be set + temporarily, reversed when exiting the context. + + """ + return self.solver.objective.direction + + @objective_direction.setter + @resettable + def objective_direction(self, value): + value = value.lower() + if value.startswith("max"): + self.solver.objective.direction = "max" + elif value.startswith("min"): + self.solver.objective.direction = "min" + else: + raise ValueError("Unknown objective direction '{}'.".format(value)) + def summary(self, solution=None, threshold=1E-8, fva=None, floatfmt='.3g'): """Print a summary of the input and output fluxes of the model. This method requires the model to have been previously solved. @@ -1064,7 +1103,7 @@ class Model(Object): address='0x0%x' % id(self), num_metabolites=len(self.metabolites), num_reactions=len(self.reactions), - objective=str(self.objective.expression), + objective=format_long_string(str(self.objective.expression), 100), compartments=", ".join( v if v else k for k, v in iteritems(self.compartments) )) diff --git a/cobra/core/reaction.py b/cobra/core/reaction.py index e2d424f..e5c3c42 100644 --- a/cobra/core/reaction.py +++ b/cobra/core/reaction.py @@ -20,6 +20,7 @@ from cobra.core.object import Object from cobra.util.context import resettable, get_context from cobra.util.solver import ( linear_reaction_coefficients, set_objective, check_solver_status) +from cobra.util.util import format_long_string # precompiled regular expressions # Matches and/or in a gene reaction rule @@ -137,11 +138,7 @@ class Reaction(Object): not associated with a model. """ if self.model is not None: - if self._forward_variable is None: - self._forward_variable = self.model.variables[ - self.id] - assert self._forward_variable.problem is self.model.solver - return self._forward_variable + return self.model.variables[self.id] else: return None @@ -155,20 +152,12 @@ class Reaction(Object): An optlang variable for the reverse flux or None if reaction is not associated with a model. """ - model = self.model - if model is not None: - if self._reverse_variable is None: - self._reverse_variable = model.variables[ - self.reverse_id] - assert self._reverse_variable.problem is self.model.solver - return self._reverse_variable + + if self.model is not None: + return self.model.variables[self.reverse_id] else: return None - def _reset_var_cache(self): - self._forward_variable = None - self._reverse_variable = None - @property def objective_coefficient(self): """ Get the coefficient for this reaction in a linear @@ -190,12 +179,10 @@ class Reaction(Object): def __copy__(self): cop = copy(super(Reaction, self)) - cop._reset_var_cache() return cop def __deepcopy__(self, memo): cop = deepcopy(super(Reaction, self), memo) - cop._reset_var_cache() return cop @property @@ -621,9 +608,15 @@ class Reaction(Object): """ new_reaction = self.copy() - new_reaction += other + if other == 0: + return new_reaction + else: + new_reaction += other + return new_reaction + __radd__ = __add__ + def __iadd__(self, other): self.add_metabolites(other._metabolites, combine=True) gpr1 = self.gene_reaction_rule.strip() @@ -911,6 +904,7 @@ class Reaction(Object): def get_compartments(self): """lists compartments the metabolites are in""" + warn('use Reaction.compartments instead', DeprecationWarning) return list(self.compartments) def _associate_gene(self, cobra_gene): @@ -1064,11 +1058,14 @@ class Reaction(Object): <td><strong>Upper bound</strong></td><td>{ub}</td> </tr> </table> - """.format(id=self.id, name=self.name, + """.format(id=format_long_string(self.id, 100), + name=format_long_string(self.name, 100), address='0x0%x' % id(self), - stoich_id=self.build_reaction_string(), - stoich_name=self.build_reaction_string(True), - gpr=self.gene_reaction_rule, + stoich_id=format_long_string( + self.build_reaction_string(), 200), + stoich_name=format_long_string( + self.build_reaction_string(True), 200), + gpr=format_long_string(self.gene_reaction_rule, 100), lb=self.lower_bound, ub=self.upper_bound) diff --git a/cobra/core/solution.py b/cobra/core/solution.py index e730565..05bba87 100644 --- a/cobra/core/solution.py +++ b/cobra/core/solution.py @@ -10,7 +10,7 @@ from warnings import warn from numpy import empty, nan from optlang.interface import OPTIMAL -from pandas import Series, DataFrame +from pandas import Series, DataFrame, option_context from cobra.util.solver import check_solver_status @@ -89,10 +89,11 @@ class Solution(object): def _repr_html_(self): if self.status == OPTIMAL: - html = ('<strong><em>Optimal</em> solution with objective value ' - '{:.3f}</strong><br>{}' - .format(self.objective_value, - self.to_frame()._repr_html_())) + with option_context('display.max_rows', 10): + html = ('<strong><em>Optimal</em> solution with objective ' + 'value {:.3f}</strong><br>{}' + .format(self.objective_value, + self.to_frame()._repr_html_())) else: html = '<strong><em>{}</em> solution</strong>'.format(self.status) return html diff --git a/cobra/flux_analysis/phenotype_phase_plane.py b/cobra/flux_analysis/phenotype_phase_plane.py index 72d36ac..5eacbd2 100644 --- a/cobra/flux_analysis/phenotype_phase_plane.py +++ b/cobra/flux_analysis/phenotype_phase_plane.py @@ -3,17 +3,15 @@ from __future__ import absolute_import, division from warnings import warn -from collections import defaultdict -from operator import itemgetter -from itertools import product from six import iteritems +from itertools import product from multiprocessing import Pool + import pandas as pd from optlang.interface import OPTIMAL - from numpy import ( nan, abs, arange, dtype, empty, int32, linspace, meshgrid, unravel_index, - zeros, array) + zeros, full) from cobra.solvers import get_solver_name, solver_dict import cobra.util.solver as sutil @@ -307,17 +305,18 @@ def calculate_phenotype_phase_plane( return data -def production_envelope(model, reactions, objective=None, c_source=None, - points=20, solver=None): +def production_envelope(model, reactions, objective=None, carbon_sources=None, + points=20, threshold=1e-7, solver=None): """Calculate the objective value conditioned on all combinations of fluxes for a set of chosen reactions - The production envelope can be used to analyze a models ability to + The production envelope can be used to analyze a model's ability to produce a given compound conditional on the fluxes for another set of - reaction, such as the uptake rates. The model is alternately optimize - with respect to minimizing and maximizing the objective and record the - obtained fluxes. Ranges to compute production is set to the effective - bounds, i.e. the minimum / maximum fluxes that can be obtained given + reactions, such as the uptake rates. The model is alternately optimized + with respect to minimizing and maximizing the objective and the + obtained fluxes are recorded. Ranges to compute production is set to the + effective + bounds, i.e., the minimum / maximum fluxes that can be obtained given current reaction bounds. Parameters @@ -325,21 +324,24 @@ def production_envelope(model, reactions, objective=None, c_source=None, model : cobra.Model The model to compute the production envelope for. reactions : list or string - A list of reactions, reaction identifiers or single reaction - objective : string, dict, model.solver.interface.Objective + A list of reactions, reaction identifiers or a single reaction. + objective : string, dict, model.solver.interface.Objective, optional The objective (reaction) to use for the production envelope. Use the - model's current objective is left missing. - c_source : cobra.Reaction or string - A reaction or reaction identifier that is the source of carbon for - computing carbon (mol carbon in output over mol carbon in input) and - mass yield (gram product over gram output). Only objectives with a - carbon containing input and output metabolite is supported. - points : int + model's current objective if left missing. + carbon_sources : list or string, optional + One or more reactions or reaction identifiers that are the source of + carbon for computing carbon (mol carbon in output over mol carbon in + input) and mass yield (gram product over gram output). Only objectives + with a carbon containing input and output metabolite is supported. + Will identify active carbon sources in the medium if none are specified. + points : int, optional The number of points to calculate production for. - solver : string + threshold : float, optional + A cut-off under which flux values will be considered to be zero. + solver : string, optional The solver to use - only here for consistency with older implementations (this argument will be removed in the future). The - solver should be set using `model.solver` directly. Only optlang + solver should be set using ``model.solver`` directly. Only optlang based solvers are supported. Returns @@ -349,18 +351,16 @@ def production_envelope(model, reactions, objective=None, c_source=None, - reaction id : one column per input reaction indicating the flux at each given point, - - flux: the objective flux + - carbon_source: identifiers of carbon exchange reactions + + A column for the maximum and minimum each for the following types: + - flux: the objective flux - carbon_yield: if carbon source is defined and the product is a single metabolite (mol carbon product per mol carbon feeding source) - - mass_yield: if carbon source is defined and the product is a single metabolite (gram product per 1 g of feeding source) - - direction: the direction of the optimization. - - Only points that give a valid solution are returned. - Examples -------- >>> import cobra.test @@ -374,168 +374,193 @@ def production_envelope(model, reactions, objective=None, c_source=None, 'based solver interfaces.') reactions = model.reactions.get_by_any(reactions) objective = model.solver.objective if objective is None else objective - result = None + data = dict() + if carbon_sources is None: + c_input = find_carbon_sources(model) + else: + c_input = model.reactions.get_by_any(carbon_sources) + if c_input is None: + data['carbon_source'] = None + elif hasattr(c_input, 'id'): + data['carbon_source'] = c_input.id + else: + data['carbon_source'] = ', '.join(rxn.id for rxn in c_input) + + size = points ** len(reactions) + for direction in ('minimum', 'maximum'): + data['flux_{}'.format(direction)] = full(size, nan, dtype=float) + data['carbon_yield_{}'.format(direction)] = full( + size, nan, dtype=float) + data['mass_yield_{}'.format(direction)] = full( + size, nan, dtype=float) + grid = pd.DataFrame(data) with model: model.objective = objective - if c_source is None: - c_input = get_c_input(model) - else: - c_input = model.reactions.get_by_any(c_source)[0] objective_reactions = list(sutil.linear_reaction_coefficients(model)) if len(objective_reactions) != 1: raise ValueError('cannot calculate yields for objectives with ' 'multiple reactions') - carbon_io = c_input, objective_reactions[0] + c_output = objective_reactions[0] min_max = fva(model, reactions, fraction_of_optimum=0) - grid = [linspace(min_max.minimum[rxn.id], min_max.maximum[rxn.id], - points, endpoint=True) for rxn in reactions] - grid_list = list(product(*grid)) - result = envelope_for_points(model, reactions, grid_list, carbon_io) - - return pd.DataFrame(result) - + min_max[min_max.abs() < threshold] = 0.0 + points = list(product(*[ + linspace(min_max.at[rxn.id, "minimum"], + min_max.at[rxn.id, "maximum"], + points, endpoint=True) for rxn in reactions])) + tmp = pd.DataFrame(points, columns=[rxn.id for rxn in reactions]) + grid = pd.concat([grid, tmp], axis=1, copy=False) + add_envelope(model, reactions, grid, c_input, c_output, threshold) + return grid + + +def add_envelope(model, reactions, grid, c_input, c_output, threshold): + if c_input is not None: + input_components = [reaction_elements(rxn) for rxn in c_input] + output_components = reaction_elements(c_output) + try: + input_weights = [reaction_weight(rxn) for rxn in c_input] + output_weight = reaction_weight(c_output) + except ValueError: + input_weights = [] + output_weight = [] + else: + input_components = [] + output_components = [] + input_weights = [] + output_weight = [] -def envelope_for_points(model, reactions, grid, carbon_io): - results = defaultdict(list) for direction in ('minimum', 'maximum'): - sense = "min" if direction == "minimum" else "max" - for point in grid: - with model: - model.solver.objective.direction = sense - for reaction, coordinate in zip(reactions, point): - reaction.bounds = coordinate, coordinate - model.slim_optimize() - if model.solver.status == OPTIMAL: - for reaction, coordinate in zip(reactions, point): - results[reaction.id].append(coordinate) - results['direction'].append(direction) - results['flux'].append(model.solver.objective.value) - if carbon_io[0] is not None: - results['carbon_yield'].append(carbon_yield(carbon_io)) - results['mass_yield'].append(mass_yield(carbon_io)) - for key, value in results.items(): - results[key] = array(value) - if carbon_io[0] is not None: - results['carbon_source'] = carbon_io[0].id - return results + with model: + model.objective_direction = direction + for i in range(len(grid)): + with model: + for rxn in reactions: + point = grid.at[i, rxn.id] + rxn.bounds = point, point + obj_val = model.slim_optimize() + if model.solver.status != OPTIMAL: + continue + grid.at[i, 'flux_{}'.format(direction)] = \ + 0.0 if abs(obj_val) < threshold else obj_val + if c_input is not None: + grid.at[i, 'carbon_yield_{}'.format(direction)] = \ + total_yield([rxn.flux for rxn in c_input], + input_components, + obj_val, + output_components) + grid.at[i, 'mass_yield_{}'.format(direction)] = \ + total_yield([rxn.flux for rxn in c_input], + input_weights, + obj_val, + output_weight) + + +def total_yield(input_fluxes, input_elements, output_flux, output_elements): + """ + Compute total output per input unit. + Units are typically mol carbon atoms or gram of source and product. -def carbon_yield(c_input_output): - """ mol product per mol carbon input + Parameters + ---------- + input_fluxes : list + A list of input reaction fluxes in the same order as the + ``input_components``. + input_elements : list + A list of reaction components which are in turn list of numbers. + output_flux : float + The output flux value. + output_elements : list + A list of stoichiometrically weighted output reaction components. Returns ------- float - the mol carbon atoms in the product (as defined by the model - objective) divided by the mol carbon in the input reactions (as - defined by the model medium) or zero in case of division by zero - arises + The ratio between output (mol carbon atoms or grams of product) and + input (mol carbon atoms or grams of source compounds). """ - c_input, c_output = c_input_output - if c_input is None: - return nan - carbon_input_flux = total_carbon_flux(c_input, consumption=True) - carbon_output_flux = total_carbon_flux(c_output, consumption=False) + carbon_input_flux = sum( + total_components_flux(flux, components, consumption=True) + for flux, components in zip(input_fluxes, input_elements)) + carbon_output_flux = total_components_flux( + output_flux, output_elements, consumption=False) try: return carbon_output_flux / carbon_input_flux except ZeroDivisionError: return nan -def mass_yield(c_input_output): - """Gram product divided by gram of carbon input source +def reaction_elements(reaction): + """ + Split metabolites into the atoms times their stoichiometric coefficients. Parameters ---------- - c_input_output : tuple - Two reactions, the one that feeds carbon to the system and the one - that produces carbon containing compound. + reaction : Reaction + The metabolic reaction whose components are desired. Returns ------- - float - gram product per 1 g of feeding source + list + Each of the reaction's metabolites' desired carbon elements (if any) + times that metabolite's stoichiometric coefficient. """ - c_input, c_output = c_input_output - if input is None: - return nan - try: - c_source, source_flux = single_flux(c_input, consumption=True) - c_product, product_flux = single_flux(c_output, consumption=False) - except ValueError: - return nan - mol_prod_mol_src = product_flux / source_flux - x = mol_prod_mol_src * c_product.formula_weight - return x / c_source.formula_weight + c_elements = [coeff * met.elements.get('C', 0) + for met, coeff in iteritems(reaction.metabolites)] + return [elem for elem in c_elements if elem != 0] -def total_carbon_flux(reaction, consumption=True): - """summed product carbon flux for a reaction +def reaction_weight(reaction): + """Return the metabolite weight times its stoichiometric coefficient.""" + if len(reaction.metabolites) != 1: + raise ValueError('Reaction weight is only defined for single ' + 'metabolite products or educts.') + met, coeff = next(iteritems(reaction.metabolites)) + return [coeff * met.formula_weight] + + +def total_components_flux(flux, components, consumption=True): + """ + Compute the total components consumption or production flux. Parameters ---------- - reaction : Reaction - the reaction to carbon return flux for - consumption : bool - flux for consumed metabolite, else produced + flux : float + The reaction flux for the components. + components : list + List of stoichiometrically weighted components. + consumption : bool, optional + Whether to sum up consumption or production fluxes. - Returns - ------- - float - reaction flux multiplied by number of carbon for the products of the - reaction """ direction = 1 if consumption else -1 - c_flux = [reaction.flux * coeff * met.elements.get('C', 0) * direction - for met, coeff in reaction.metabolites.items()] + c_flux = [elem * flux * direction for elem in components] return sum([flux for flux in c_flux if flux > 0]) -def single_flux(reaction, consumption=True): - """flux into single product for a reaction - - only defined for reactions with single products +def find_carbon_sources(model): + """ + Find all active carbon source reactions. Parameters ---------- - reaction : Reaction - the reaction to product flux for - consumption : bool - flux for consumed metabolite, else produced + model : Model + A genome-scale metabolic model. Returns ------- - tuple - metabolite, flux for the metabolite - """ - if len(list(reaction.metabolites)) != 1: - raise ValueError('product flux only defined for single metabolite ' - 'reactions') - met, coeff = next(iteritems(reaction.metabolites)) - direction = 1 if consumption else -1 - return met, reaction.flux * coeff * direction - - -def get_c_input(model): - """ carbon source reactions + list + The medium reactions with carbon input flux. - Returns - ------- - Reaction - The medium reaction with highest input carbon flux """ try: model.slim_optimize(error_value=None) except OptimizationError: - return None + return [] reactions = model.reactions.get_by_any(list(model.medium)) - reactions_fluxes = [(rxn, total_carbon_flux(rxn, consumption=True)) - for rxn in reactions] - source_reactions = [(rxn, c_flux) for rxn, c_flux - in reactions_fluxes if c_flux > 0] - try: - return max(source_reactions, key=itemgetter(1))[0] - except ValueError: - return None + reactions_fluxes = [ + (rxn, total_components_flux(rxn.flux, reaction_elements(rxn), + consumption=True)) for rxn in reactions] + return [rxn for rxn, c_flux in reactions_fluxes if c_flux > 0] diff --git a/cobra/flux_analysis/summary.py b/cobra/flux_analysis/summary.py index c75f4b6..a53b205 100644 --- a/cobra/flux_analysis/summary.py +++ b/cobra/flux_analysis/summary.py @@ -10,16 +10,10 @@ from tabulate import tabulate from cobra.flux_analysis.variability import flux_variability_analysis from cobra.util.solver import linear_reaction_coefficients +from cobra.util.util import format_long_string from cobra.core import get_solution -def format_long_string(string, max_length): - if len(string) > max_length: - string = string[:max_length - 3] - string += '...' - return string - - def metabolite_summary(met, solution=None, threshold=0.01, fva=False, floatfmt='.3g'): """Print a summary of the reactions which produce and consume this @@ -145,7 +139,7 @@ def model_summary(model, solution=None, threshold=1E-8, fva=None, """ objective_reactions = linear_reaction_coefficients(model) boundary_reactions = model.exchanges - summary_rxns = list(objective_reactions.keys()) + boundary_reactions + summary_rxns = set(objective_reactions.keys()).union(boundary_reactions) if solution is None: model.slim_optimize(error_value=None) diff --git a/cobra/manipulation/validate.py b/cobra/manipulation/validate.py index cf39969..e9f217f 100644 --- a/cobra/manipulation/validate.py +++ b/cobra/manipulation/validate.py @@ -50,10 +50,6 @@ def check_reaction_bounds(model): def check_metabolite_compartment_formula(model): errors = [] for met in model.metabolites: - if met.compartment is not None and \ - met.compartment not in model.compartments: - errors.append("Metabolite '%s' compartment '%s' not found" % - (met.id, met.compartment)) if met.formula is not None and len(met.formula) > 0: if not met.formula.isalnum(): errors.append("Metabolite '%s' formula '%s' not alphanumeric" % diff --git a/cobra/test/data/iJO1366.pickle b/cobra/test/data/iJO1366.pickle index 6346508..2e89067 100644 Binary files a/cobra/test/data/iJO1366.pickle and b/cobra/test/data/iJO1366.pickle differ diff --git a/cobra/test/data/mini.pickle b/cobra/test/data/mini.pickle index eb9d561..af9ef9c 100644 Binary files a/cobra/test/data/mini.pickle and b/cobra/test/data/mini.pickle differ diff --git a/cobra/test/data/raven.pickle b/cobra/test/data/raven.pickle index 9750952..adfaf9b 100644 Binary files a/cobra/test/data/raven.pickle and b/cobra/test/data/raven.pickle differ diff --git a/cobra/test/data/salmonella.pickle b/cobra/test/data/salmonella.pickle index 3848bdb..ba0568c 100644 Binary files a/cobra/test/data/salmonella.pickle and b/cobra/test/data/salmonella.pickle differ diff --git a/cobra/test/data/textbook_solution.pickle b/cobra/test/data/textbook_solution.pickle index f5e453d..a9633a3 100644 Binary files a/cobra/test/data/textbook_solution.pickle and b/cobra/test/data/textbook_solution.pickle differ diff --git a/cobra/test/test_flux_analysis.py b/cobra/test/test_flux_analysis.py index 773d10b..aef830e 100644 --- a/cobra/test/test_flux_analysis.py +++ b/cobra/test/test_flux_analysis.py @@ -620,6 +620,10 @@ class TestCobraFluxAnalysis: model.summary() self.check_in_line(out.getvalue(), expected_entries) + with model: + model.objective = model.exchanges[0] + model.summary() + @pytest.mark.parametrize("fraction", [0.95]) def test_model_summary_with_fva(self, model, opt_solver, fraction): if opt_solver == "optlang-gurobi": @@ -834,7 +838,7 @@ class TestProductionEnvelope: def test_envelope_one(self, model): df = production_envelope(model, ["EX_o2_e"]) - assert abs(sum(df.flux) - 9.342) < 0.001 + assert numpy.isclose(df["flux_maximum"].sum(), 9.342, atol=1e-3) def test_envelope_multi_reaction_objective(self, model): obj = {model.reactions.EX_ac_e: 1, @@ -842,12 +846,22 @@ class TestProductionEnvelope: with pytest.raises(ValueError): production_envelope(model, "EX_o2_e", obj) + @pytest.mark.parametrize("variables, num", [ + (["EX_glc__D_e"], 30), + (["EX_glc__D_e", "EX_o2_e"], 20), + (["EX_glc__D_e", "EX_o2_e", "EX_ac_e"], 10) + ]) + def test_multi_variable_envelope(self, model, variables, num): + df = production_envelope(model, variables, points=num) + assert len(df) == num ** len(variables) + def test_envelope_two(self, model): df = production_envelope(model, ["EX_glc__D_e", "EX_o2_e"], objective="EX_ac_e") - assert abs(numpy.sum(df.carbon_yield) - 83.579) < 0.001 - assert abs(numpy.sum(df.flux) - 1737.466) < 0.001 - assert abs(numpy.sum(df.mass_yield) - 82.176) < 0.001 + assert numpy.isclose(df["flux_maximum"].sum(), 1737.466, atol=1e-3) + assert numpy.isclose(df["carbon_yield_maximum"].sum(), 83.579, + atol=1e-3) + assert numpy.isclose(df["mass_yield_maximum"].sum(), 82.176, atol=1e-3) class TestReactionUtils: diff --git a/cobra/test/test_manipulation.py b/cobra/test/test_manipulation.py index 40b98d5..a7d37f3 100644 --- a/cobra/test/test_manipulation.py +++ b/cobra/test/test_manipulation.py @@ -225,10 +225,9 @@ class TestManipulation: assert rxns.EX_h_e.annotation["SBO"] == "SBO:0000628" def test_validate_formula_compartment(self, model): - model.metabolites[1].compartment = "fake" model.metabolites[1].formula = "(a*.bcde)" errors = check_metabolite_compartment_formula(model) - assert len(errors) == 2 + assert len(errors) == 1 def test_validate_mass_balance(self, model): assert len(check_mass_balance(model)) == 0 diff --git a/cobra/test/test_model.py b/cobra/test/test_model.py index 745d5d7..b272264 100644 --- a/cobra/test/test_model.py +++ b/cobra/test/test_model.py @@ -287,6 +287,11 @@ class TestReactions: assert gene is not model.genes.get_by_id(gene.id) assert gene.model is not model + def test_radd(self, model): + new = sum([model.reactions.PGI, model.reactions.EX_h2o_e]) + assert new._model is not model + assert len(new.metabolites) == 3 + def test_mul(self, model): new = model.reactions.PGI * 2 assert set(new.metabolites.values()) == {-2, 2} @@ -395,6 +400,15 @@ class TestCobraModel: def test_compartments(self, model): assert set(model.compartments) == {"c", "e"} + model = Model("test", "test") + met_c = Metabolite("a_c", compartment="c") + met_e = Metabolite("a_e", compartment="e") + rxn = Reaction("foo") + rxn.add_metabolites({met_e: -1, met_c: 1}) + model.add_reactions([rxn]) + assert model.compartments == {'c': '', 'e': ''} + model.compartments = {'c': 'cytosol'} + assert model.compartments == {'c': 'cytosol', 'e': ''} def test_add_reaction(self, model): old_reaction_count = len(model.reactions) @@ -766,6 +780,19 @@ class TestCobraModel: benchmark(benchmark_change_objective) + def test_get_objective_direction(self, model): + assert model.objective_direction == "max" + value = model.slim_optimize() + assert numpy.isclose(value, 0.874, 1e-3) + + def test_set_objective_direction(self, model): + with model: + model.objective_direction = "min" + assert model.objective_direction == "min" + value = model.slim_optimize() + assert value == 0.0 + assert model.objective_direction == "max" + def test_slim_optimize(self, model): with model: assert model.slim_optimize() > 0.872 diff --git a/cobra/util/util.py b/cobra/util/util.py index 3aeb377..6299a50 100644 --- a/cobra/util/util.py +++ b/cobra/util/util.py @@ -3,6 +3,13 @@ from __future__ import absolute_import +def format_long_string(string, max_length=50): + if len(string) > max_length: + string = string[:max_length - 3] + string += '...' + return string + + class AutoVivification(dict): """Implementation of perl's autovivification feature. Checkout http://stackoverflow.com/a/652284/280182 """ diff --git a/release-notes/0.8.2.md b/release-notes/0.8.2.md index d2ecece..a5251a2 100644 --- a/release-notes/0.8.2.md +++ b/release-notes/0.8.2.md @@ -2,8 +2,6 @@ ## Fixes -- the Solution class no longer contains links progenitor model's - reactions and metabolites - Guarantee that sampler._reproject always returns a feasible point and will not attempt to reproject already feasible points. [#564](https://github.com/opencobra/cobrapy/pull/564) @@ -15,3 +13,14 @@ model or `ValueError` is raised. - Fix use of underscores in key/value pairs in legacy sbml notes. [#547](https://github.com/opencobra/cobrapy/issues/547) + +## Backwards incompatible changes + +- the Solution class no longer contains links progenitor model's + reactions and metabolites. Removed since it those can change after + the solution has been computed making them erroneous. This of course + implies that `Solution` constructor changes to: + ``` + def __init__(self, objective_value, status, fluxes, reduced_costs=None, + shadow_prices=None, **kwargs): + ``` diff --git a/release-notes/0.9.0.md b/release-notes/0.9.0.md new file mode 100644 index 0000000..178a215 --- /dev/null +++ b/release-notes/0.9.0.md @@ -0,0 +1,30 @@ +# Release notes for cobrapy 0.9.0 + +## Fixes + +- `Model.compartment` is now a dynamic property fetching the + compartments from all metabolites therefore always + up-to-date. Assigning a dictionary to the same property updates the + internal dictionary of compartment descriptions. This change removes + the need for the check for missing compartments from + `validation.check_metabolite_compartment_formula`. +- Excessively long output of html representations in jupyter notebooks are now abbreviated [#577](https://github.com/opencobra/cobrapy/pull/577). +- Reaction forward and reverse variables are no longer cached with those object. No visible effect but simplifies the code. +- Fix bug in summary methods when used with exchange reaction sas objective. [#595](https://github.com/opencobra/cobrapy/pull/595). + +## New features + +- `Model.objective_direction` is a new revertible property to set maximization / minimization using the context manager. +- Change output of `production_envelope` to wide data frame format [#587](https://github.com/opencobra/cobrapy/pull/587). Also allow multiple carbon source reactions and better handling of zero-division exceptions. +- Enable summing lists of reactions, see [#596](https://github.com/opencobra/cobrapy/pull/596) + +## Deprecated features + +- `Model.get_metabolite_compartments` is deprecated (use + `Model.compartments` instead). +- `Reaction.get_compartments` is deprecated (use + `Reaction.compartments` instead). + +## Backwards incompatible changes + +- The format of the dataframe `production_envelope` changed listing max and min on different columns instead of the same. diff --git a/release-notes/next-release.md b/release-notes/next-release.md index c622d35..7d5acdd 100644 --- a/release-notes/next-release.md +++ b/release-notes/next-release.md @@ -5,3 +5,5 @@ ## New features ## Deprecated features + +## Backwards incompatible changes diff --git a/scripts/compare-benchmark.py b/scripts/compare-benchmark.py index c203dc4..686c659 100644 --- a/scripts/compare-benchmark.py +++ b/scripts/compare-benchmark.py @@ -4,6 +4,8 @@ import argparse import re from os.path import basename +pd.set_option('display.width', 200) + def benchmark_to_df(json_file): with open(json_file) as jf: @@ -18,9 +20,13 @@ def benchmark_to_df(json_file): if __name__ == "__main__": parser = argparse.ArgumentParser(description=""" - compare cobrapy benchmarks.""") - parser.add_argument('first', help='first branch to compare') - parser.add_argument('second', help='second branch to compare') + compare cobrapy benchmarks. + Run pytest with + pytest --benchmark-save=without-cache --benchmark-min-rounds=20 + then compare saved json files with this script. + """) + parser.add_argument('first', help='first json file') + parser.add_argument('second', help='second json file') args = parser.parse_args() first = benchmark_to_df(args.first) @@ -31,5 +37,5 @@ if __name__ == "__main__": both = pd.merge(first, second, how="inner", on="test", suffixes=(first_name, second_name)) both["fraction"] = both.iloc[:, 2] / both.iloc[:, 1] - both.sort_values(by="fraction") - print(both) + print(both.sort_values(by="fraction")) + diff --git a/setup.cfg b/setup.cfg index 6c7b437..73d72f1 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,5 +1,5 @@ [bumpversion] -current_version = 0.8.2 +current_version = 0.9.0 commit = True tag = True parse = (?P<major>\d+) diff --git a/setup.py b/setup.py index be73b98..bc1dcb6 100644 --- a/setup.py +++ b/setup.py @@ -144,7 +144,7 @@ except: setup( name="cobra", - version="0.8.2", + version="0.9.0", packages=find_packages(), setup_requires=setup_requirements, install_requires=["future", "swiglpk", "optlang>=1.2.1", -- Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-cobra.git _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
