diff --git a/CHANGELOG.md b/CHANGELOG.md index 6e57433..065508a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,18 @@ ## Unreleased +* New: `level_complex_concs` and `level_complex_abundances` give every enzyme + complex one abundance, the median per copy over its measured subunits, so a + single badly quantified subunit no longer sets a complex's capacity. + `complex_families` groups isozyme variants into complexes and leaves shared + subunits and electron carriers as measured. See `docs/relaxation_methods.md`. +* New: `relax_enz_concs` relaxes enzyme caps (and, if need be, the protein + pool) by the least total fold increase that reaches the target growth, in one + LP, so the result does not depend on the order in which limiting enzymes are + found. +* `flexibilize_enz_concs` gained `per_complex=True`: complexes are ranked and + raised as one unit, keeping levelled subunits levelled. Off by default. + * Fix: `write_dlkcat_input` and the OpenKineticsPredictor input builder raised `TypeError: cannot use 'list' as a set element` on any model that had been saved and reloaded. cobra returns an annotation value as a list diff --git a/docs/relaxation_methods.md b/docs/relaxation_methods.md index 29eb2d8..3b823d5 100644 --- a/docs/relaxation_methods.md +++ b/docs/relaxation_methods.md @@ -1,4 +1,4 @@ -# `flexibilize_enz_concs` vs `relax_proteomics_greedy` +# Relaxing proteomics constraints: `flexibilize_enz_concs`, `relax_proteomics_greedy` and `relax_enz_concs` Both functions answer the same question โ€” *the measured proteomics constraints make the ecModel grow too slowly (or not at all); which enzyme bounds do we @@ -6,6 +6,8 @@ loosen, and by how much, to restore feasible growth?* โ€” but they are **not** duplicates. They use different selection signals, different relaxation step sizes, and produce results of very different stringency. This note records the differences and the recommendation (keep both; do not merge). +`relax_enz_concs`, a third method that solves the same question as one LP, and +the per-complex options for all of them are described at the end. ## Side by side @@ -50,3 +52,68 @@ follow-ups would make them more interchangeable: the over-relaxation point in the modeling-intent review, ยง3.3). 2. Align the failure convention โ€” currently one returns a `converged=False` result and the other raises; pick one so callers don't special-case each. + +## `relax_enz_concs`: the least total relaxation, in one LP + +Both methods above relax one enzyme at a time: rank, raise the top one, +re-solve, repeat. Which enzymes end up raised then depends on the ranking +signal and on the order in which enzymes become limiting. On a large model +with noisy proteomics, a small change to the ranking (for example, control +coefficient per mg against gain per fold) can change which enzymes are +raised and, downstream, predictions that depend on the remaining caps. + +`relax_enz_concs` asks the question once. Every measured cap becomes +`usage <= conc * (1 + x)`, the protein pool may grow by a fold too, and one LP +minimises the sum of the `x` subject to `growth >= exp_growth`: + +- `weight="relative"` (default) minimises summed fold increases, so an enzyme + measured at almost nothing is not relaxed freely because its absolute cost is + small; `weight="absolute"` minimises the mg added. +- The pool's fold costs `pool_weight` (default 1000) times an enzyme's, so the + pool grows only where raising enzymes cannot reach the target, the role of + the pool fallback in `flexibilize_enz_concs`. +- Folds count from the current caps, so a cap raised earlier is never lowered. +- If no relaxation reaches the target, it raises `ValueError`: growth is then + not limited by protein. + +It returns a `FlexEnzResult` like `flexibilize_enz_concs` (`frequence` is 1 for +every entry). Use it when the outcome should not depend on the order of +relaxation, and to compare with the iterative methods: if they disagree on +which enzymes need raising, conclusions that hinge on those enzymes are +fragile. + +## Enzyme complexes + +Absolute proteomics recovers the subunits of one complex unevenly: subunits +present 1:1 can come out several-fold apart while their replicates agree within +a few percent. Capped one by one, a complex can then carry no more flux than +its least abundant measurement allows, and relaxation methods that raise one +subunit at a time leave its subunits at unrelated levels, or raise a different +subunit in each condition. + +- `level_complex_concs` (after `fill_enz_concs`, before `constrain_enz_concs`) + gives every complex one abundance: the median per copy over its measured + subunits, scaled back to each subunit by its copies in `ec.rxn_enz_mat`; + unmeasured subunits are filled. `level_complex_abundances` does the same on + molar data before it is written into the model. +- `complex_families` defines the complexes. Two enzyme sets are one complex when + they share at least two subunits and half the larger set (isozyme variants, + the same complex on several reactions). A protein in three or more complexes, + in two with equal weight, or that also catalyses a reaction alone is *shared* + and keeps its own measurement: electron carriers (cytochrome c), thioredoxin, + or the E3 subunit of the 2-oxo-acid dehydrogenases, which would otherwise + chain unrelated complexes together. A subunit in two complexes belongs to the + one with more reactions (core subunits that also appear in an assembly step). +- `flexibilize_enz_concs(per_complex=True)` and `relax_enz_concs(per_complex=True)` + raise the measured subunits of a complex by one factor, so levelled + concentrations stay levelled. In `flexibilize_enz_concs` complexes are ranked + by the growth gained per fold of raising all their caps (the usage reactions' + reduced costs times the caps), which leaves out the protein pool's price and + does not depend on which of several simultaneously binding subunits the solver + puts the dual on. + +Levelling relies on the stoichiometry in `ec.rxn_enz_mat`; where +`apply_complex_data` found no exact ComplexPortal match, every subunit counts +as one copy. Per-complex relaxation without levelling raises every subunit by +the fold its worst-measured subunit needs, which inflates the well-measured +ones; use the two together. diff --git a/src/geckopy/__init__.py b/src/geckopy/__init__.py index 90d3e8a..ad0a949 100644 --- a/src/geckopy/__init__.py +++ b/src/geckopy/__init__.py @@ -111,16 +111,22 @@ # limit_proteins from geckopy.limit_proteins import ( + ComplexFamilies, FlexEnzResult, GreedyRelaxResult, + LevelComplexResult, RelaxationStep, apply_flux_data_constraints, calculate_f_factor, + complex_families, constrain_enz_concs, constrain_flux_data, # deprecated; alias of apply_flux_data_constraints fill_enz_concs, flexibilize_enz_concs, get_conc_control_coeffs, + level_complex_abundances, + level_complex_concs, + relax_enz_concs, relax_proteomics_greedy, ) @@ -206,15 +212,21 @@ "sensitivity_tuning", "truncate_values", # limit_proteins + "ComplexFamilies", "FlexEnzResult", "GreedyRelaxResult", + "LevelComplexResult", "RelaxationStep", "apply_flux_data_constraints", "calculate_f_factor", + "complex_families", "constrain_enz_concs", "fill_enz_concs", "flexibilize_enz_concs", "get_conc_control_coeffs", + "level_complex_abundances", + "level_complex_concs", + "relax_enz_concs", "relax_proteomics_greedy", # utilities "AddNewRxnsResult", diff --git a/src/geckopy/limit_proteins/__init__.py b/src/geckopy/limit_proteins/__init__.py index 79f5f1a..03b96dd 100644 --- a/src/geckopy/limit_proteins/__init__.py +++ b/src/geckopy/limit_proteins/__init__.py @@ -11,6 +11,14 @@ from .fill_enz_concs import fill_enz_concs from .flexibilize_enz_concs import FlexEnzResult, flexibilize_enz_concs from .get_conc_control_coeffs import get_conc_control_coeffs +from .level_complex_concs import ( + ComplexFamilies, + LevelComplexResult, + complex_families, + level_complex_abundances, + level_complex_concs, +) +from .relax_enz_concs import relax_enz_concs from .relax_proteomics_greedy import ( GreedyRelaxResult, RelaxationStep, @@ -18,15 +26,21 @@ ) __all__ = [ + "ComplexFamilies", "FlexEnzResult", + "LevelComplexResult", "GreedyRelaxResult", "RelaxationStep", "apply_flux_data_constraints", "calculate_f_factor", + "complex_families", "constrain_enz_concs", "constrain_flux_data", # deprecated; alias of apply_flux_data_constraints "fill_enz_concs", "flexibilize_enz_concs", "get_conc_control_coeffs", + "level_complex_abundances", + "level_complex_concs", + "relax_enz_concs", "relax_proteomics_greedy", ] diff --git a/src/geckopy/limit_proteins/flexibilize_enz_concs.py b/src/geckopy/limit_proteins/flexibilize_enz_concs.py index 04e5426..4905278 100644 --- a/src/geckopy/limit_proteins/flexibilize_enz_concs.py +++ b/src/geckopy/limit_proteins/flexibilize_enz_concs.py @@ -37,6 +37,7 @@ USAGE_PREFIX, ) from .get_conc_control_coeffs import get_conc_control_coeffs +from .level_complex_concs import ComplexFamilies, complex_families if TYPE_CHECKING: from ..ec_model.ec_model import EcModel @@ -103,6 +104,7 @@ def flexibilize_enz_concs( iter_per_enzyme: int = 5, bio_rxn: Optional[str] = None, verbose: bool = True, + per_complex: bool = False, ) -> FlexEnzResult: """Iteratively relax enzyme concentration constraints until the model reaches ``exp_growth``. @@ -153,6 +155,15 @@ def flexibilize_enz_concs( adapter; pass it (with ``exp_growth``) to run without an adapter. verbose Whether per-iteration progress is logged at INFO. + per_complex + Treat every complex (see ``complex_families``) as one enzyme: + rank enzymes and complexes by the growth gained per fold of + raising all their caps together (a levelled complex binds on all + subunits at once, so no single subunit shows a gain), raise every + measured subunit of the picked complex by the same factor, count + iterations per complex, and keep one factor per complex in the + refinement pass, so that concentrations levelled by + ``level_complex_concs`` stay levelled. Off by default. Returns ------- @@ -203,6 +214,23 @@ def flexibilize_enz_concs( proteins = [model.ec.enzymes[i] for i in measured_idx] n_proteins = len(proteins) frequence = np.zeros(n_proteins, dtype=int) + touched = np.zeros(n_proteins, dtype=bool) + position = {p: i for i, p in enumerate(proteins)} + families: Optional[ComplexFamilies] = complex_families(model) if per_complex else None + # Indices into ``proteins`` of each unit that is relaxed together: a + # complex's measured, non-shared subunits, or a single enzyme. + units: list[list[int]] = [] + unit_of = np.zeros(n_proteins, dtype=int) + assigned: set[int] = set() + for i, protein in enumerate(proteins): + if i in assigned: + continue + members = [position[m] for m in families.members(protein) if m in position] if families else [i] + for j in members: + unit_of[j] = len(units) + assigned.update(members) + units.append(members) + unit_frequence = np.zeros(len(units), dtype=int) initial_sol = model.optimize() if ( @@ -238,14 +266,17 @@ def flexibilize_enz_concs( ) flex_break = True break - _, control_coeffs = get_conc_control_coeffs( - model, - proteins=proteins, - fold_change=fold_change, - limit=_CONTROL_COEFF_LIMIT, - ) + if families is not None: + scores = _unit_scores(model, proteins, units, fold_change, _CONTROL_COEFF_LIMIT) + else: + _, scores = get_conc_control_coeffs( + model, + proteins=proteins, + fold_change=fold_change, + limit=_CONTROL_COEFF_LIMIT, + ) - if not (control_coeffs > 0).any(): + if not (scores > 0).any(): pool_old, pool_new, pool_pred_growth = _try_relax_pool( model, bio_rxn, exp_growth, pred_growth, verbose, ) @@ -253,23 +284,32 @@ def flexibilize_enz_concs( pred_growth = pool_pred_growth break - max_idx = int(np.argmax(control_coeffs)) + if families is not None: + unit = int(np.argmax(scores)) + else: + unit = int(unit_of[int(np.argmax(scores))]) + members = units[unit] + max_idx = members[0] frequence[max_idx] += 1 - if frequence[max_idx] > iter_per_enzyme: + unit_frequence[unit] += 1 + if unit_frequence[unit] > iter_per_enzyme: logger.warning( - "flexibilize_enz_concs: iter limit reached for protein %r; " + "flexibilize_enz_concs: iter limit reached for %s %r; " "consider revising its kcat.", - proteins[max_idx], + "complex" if len(members) > 1 else "protein", + [proteins[j] for j in members] if len(members) > 1 else proteins[max_idx], ) flex_break = True break - increase = fold_change * frequence[max_idx] - ec_idx = measured_idx[max_idx] - usage_rxn = model.reactions.get_by_id( - f"{USAGE_PREFIX}{proteins[max_idx]}" - ) - usage_rxn.upper_bound = float(model.ec.concs[ec_idx]) * (1.0 + increase) + increase = fold_change * unit_frequence[unit] + for i in members: + touched[i] = True + usage_rxn = model.reactions.get_by_id(f"{USAGE_PREFIX}{proteins[i]}") + usage_rxn.upper_bound = max( + usage_rxn.upper_bound, + float(model.ec.concs[measured_idx[i]]) * (1.0 + increase), + ) new_sol = model.optimize() pred_growth = float(new_sol.objective_value or 0.0) @@ -282,6 +322,7 @@ def flexibilize_enz_concs( return _build_result( model, proteins, measured_idx, frequence, flex_break, exp_growth, pool_old, pool_new, bio_rxn_id, + touched=touched, families=families, ) @@ -345,9 +386,13 @@ def _build_result( pool_old: Optional[float], pool_new: Optional[float], bio_rxn_id: str, + touched: Optional[np.ndarray] = None, + families: Optional[ComplexFamilies] = None, ) -> FlexEnzResult: """Build the FlexEnzResult and apply the post-loop UB refinement.""" flex_mask = frequence > 0 + if touched is not None: + flex_mask = flex_mask | touched flex_indices = np.where(flex_mask)[0] if len(flex_indices) == 0: return _maybe_add_pool(FlexEnzResult(), pool_old, pool_new) @@ -377,6 +422,8 @@ def _build_result( dtype=float, ) + if families is not None: + new_concs = _per_complex(flex_proteins, old_concs_arr, new_concs, families) keep_old_mask = new_concs < old_concs_arr for i, rxn_id in enumerate(usage_rxn_ids): ub = ( @@ -412,6 +459,90 @@ def _build_result( return _finalize_and_sort(_maybe_add_pool(result, pool_old, pool_new)) +def _unit_scores( + model: "EcModel", + proteins: list[str], + units: list[list[int]], + fold_change: float, + limit: float, +) -> np.ndarray: + """Growth gained per fold of raising each unit's caps together. + + With duals, ``sum(reduced cost of usage_prot_i * its cap)`` over the + unit's subunits that run at more than ``limit`` of their cap. Unlike + the ``prot_`` shadow price, the usage reaction's reduced cost + leaves out the protein pool's price, which would otherwise be counted + once per subunit; and the weighted sum is the same for every optimal + dual when a levelled complex binds on all subunits at once. When no + dual is positive (no duals exposed, or the protein pool binds and the + solver puts the whole price on it), one solve per candidate unit with + all its caps raised by ``fold_change`` and the pool lifted. + """ + from .get_conc_control_coeffs import _solution_is_optimal, _solver_supports_duals + + scores = np.zeros(len(units)) + sol = model.optimize() + if not _solution_is_optimal(sol): + return scores + candidates = [] + for k, members in enumerate(units): + reactions = [model.reactions.get_by_id(f"{USAGE_PREFIX}{proteins[j]}") for j in members] + active = [r for r in reactions if r.upper_bound > 0 and float(sol.fluxes[r.id]) / r.upper_bound > limit] + if active: + candidates.append((k, reactions, active)) + if _solver_supports_duals(model): + reduced = sol.reduced_costs + for k, _, active in candidates: + scores[k] = sum(abs(float(reduced[r.id])) * r.upper_bound for r in active) + scores[scores < 1e-12] = 0.0 + if (scores > 0).any(): + return scores + # No dual singles a unit out (none exposed, or the protein pool binds as + # well and carries the whole price): probe each candidate with the pool + # lifted, so the probe finds the caps that limit once the pool is relaxed. + pool_rxn = model.reactions.get_by_id(POOL_EXCHANGE_ID) + with model: + pool_rxn.upper_bound = max(pool_rxn.upper_bound, 1000.0) + base = model.optimize() + if not _solution_is_optimal(base): + return scores + baseline = float(base.objective_value) + for k, reactions, _ in candidates: + with model: + for r in reactions: + r.upper_bound = r.upper_bound * fold_change + probe = model.optimize() + if _solution_is_optimal(probe): + gain = float(probe.objective_value) - baseline + scores[k] = gain / (fold_change - 1.0) if gain > _GROWTH_DELTA_THRESHOLD else 0.0 + return scores + + +def _per_complex( + proteins: list[str], + old: np.ndarray, + new: np.ndarray, + families: ComplexFamilies, +) -> np.ndarray: + """One factor per complex: every flexibilized subunit takes the largest new/old of its complex.""" + out = new.copy() + index = {p: i for i, p in enumerate(proteins)} + seen: set[str] = set() + for protein in proteins: + if protein in seen: + continue + group = [m for m in families.members(protein) if m in index] + seen.update(group) + if len(group) < 2: + continue + ratios = [new[index[m]] / old[index[m]] for m in group if old[index[m]] > 0] + factor = max(ratios) if ratios else 1.0 + for m in group: + if old[index[m]] > 0: + out[index[m]] = old[index[m]] * factor + return out + + def _maybe_add_pool( result: FlexEnzResult, pool_old: Optional[float], diff --git a/src/geckopy/limit_proteins/level_complex_concs.py b/src/geckopy/limit_proteins/level_complex_concs.py new file mode 100644 index 0000000..d060f04 --- /dev/null +++ b/src/geckopy/limit_proteins/level_complex_concs.py @@ -0,0 +1,278 @@ +"""One abundance per enzyme complex. + +Absolute proteomics recovers the subunits of a complex unevenly, so +after ``fill_enz_concs`` the subunits of one complex can disagree several +fold about how much of it there is. Capped one by one, the complex can +then carry no more flux than its least abundant subunit allows. This +module groups the enzymes of every multi-subunit reaction into complexes +and gives each complex one abundance, scaled to every subunit by its +stoichiometry. + +A complex is a set of enzymes that catalyses a reaction together, taken +from ``ec.rxn_enz_mat``. Two enzyme sets are one complex when they share +at least two subunits and at least half of the larger set (isozyme +variants that swap a subunit, the same complex on different reactions). +A single shared protein never merges two sets, so electron carriers and +shared subunits (cytochrome c, thioredoxin, the E3 subunit of the +2-oxo-acid dehydrogenases) do not chain unrelated complexes together. + +A subunit in exactly two complexes belongs to the one whose enzyme sets +catalyse more reactions (the core subunits COX1-3 of complex IV also +appear in one assembly reaction). A subunit in three or more complexes, +in two with equal weight, or that also catalyses a reaction on its own +is *shared*: it is left as measured and does not take part in the +summary. +""" +from __future__ import annotations + +import logging +from dataclasses import dataclass, field +from typing import TYPE_CHECKING, Optional, Sequence + +import numpy as np + +if TYPE_CHECKING: + from ..ec_model.ec_model import EcModel + +logger = logging.getLogger(__name__) + +_STATISTICS = {"median": np.median, "mean": np.mean} + + +@dataclass +class ComplexFamilies: + """The complexes of an ecModel and their shared subunits. + + Attributes + ---------- + complexes + One dict per complex, ``{enzyme: copies}`` over every subunit + that appears in any of its enzyme sets. + owner + For every subunit that belongs to a complex, the index of that + complex in ``complexes``. + shared + Subunits of a complex that belong to none: carriers and subunits + shared between complexes. + """ + + complexes: list[dict[str, float]] = field(default_factory=list) + owner: dict[str, int] = field(default_factory=dict) + shared: set[str] = field(default_factory=set) + + def owned(self, index: int) -> dict[str, float]: + """``{enzyme: copies}`` of the subunits complex ``index`` owns.""" + return {e: c for e, c in self.complexes[index].items() if self.owner.get(e) == index} + + def members(self, enzyme: str) -> list[str]: + """The subunits of the complex ``enzyme`` belongs to, or ``[enzyme]``.""" + index = self.owner.get(enzyme) + return [enzyme] if index is None else list(self.owned(index)) + + +@dataclass +class LevelComplexResult: + """What ``level_complex_concs`` changed. + + Attributes + ---------- + adjusted + Measured subunits whose concentration was replaced by the complex's. + added + Unmeasured subunits given the complex's concentration. + shared + Subunits left as measured because they belong to more than one complex. + skipped + Complexes (as sorted subunit tuples) without any measured, non-shared subunit. + """ + + adjusted: list[str] = field(default_factory=list) + added: list[str] = field(default_factory=list) + shared: list[str] = field(default_factory=list) + skipped: list[tuple[str, ...]] = field(default_factory=list) + + +def complex_families( + model: "EcModel", min_overlap: float = 0.5, min_shared: int = 2, +) -> ComplexFamilies: + """Group the enzymes of multi-subunit reactions into complexes. + + Parameters + ---------- + model + An EcModel with ``ec.rxn_enz_mat`` (copies per subunit). + min_overlap + Two enzyme sets are one complex when they share at least this + fraction of the larger set's subunits... + min_shared + ...and at least this many subunits. + """ + matrix = model.ec.rxn_enz_mat.tocsr() + enzymes = list(model.ec.enzymes) + sets: dict[frozenset, dict[str, float]] = {} + reactions: dict[frozenset, int] = {} + alone: set[str] = set() + for row in range(matrix.shape[0]): + start, end = matrix.indptr[row], matrix.indptr[row + 1] + present = {enzymes[j]: float(v) for j, v in zip(matrix.indices[start:end], matrix.data[start:end]) if v} + if len(present) == 1: + alone.update(present) + elif present: + sets.setdefault(frozenset(present), present) + reactions[frozenset(present)] = reactions.get(frozenset(present), 0) + 1 + + keys = list(sets) + parent = list(range(len(keys))) + + def find(i: int) -> int: + while parent[i] != i: + parent[i] = parent[parent[i]] + i = parent[i] + return i + + by_enzyme: dict[str, list[int]] = {} + for i, key in enumerate(keys): + for enzyme in key: + by_enzyme.setdefault(enzyme, []).append(i) + for indices in by_enzyme.values(): + for x in range(len(indices)): + for y in range(x + 1, len(indices)): + i, j = indices[x], indices[y] + common = len(keys[i] & keys[j]) + if common >= min_shared and common >= min_overlap * max(len(keys[i]), len(keys[j])): + parent[find(i)] = find(j) + + grouped: dict[int, dict[str, float]] = {} + weight: dict[int, dict[str, int]] = {} + for i, key in enumerate(keys): + w = weight.setdefault(find(i), {}) + for enzyme in key: + w[enzyme] = w.get(enzyme, 0) + reactions[key] + complex_ = grouped.setdefault(find(i), {}) + for enzyme, copies in sets[key].items(): + if enzyme in complex_ and complex_[enzyme] != copies: + logger.warning( + "complex_families: %s has %g copies in one variant and %g in another; " + "the first is used", enzyme, complex_[enzyme], copies, + ) + complex_.setdefault(enzyme, copies) + roots = list(grouped) + complexes = [grouped[r] for r in roots] + + found: dict[str, list[tuple[int, int]]] = {} + for index, root in enumerate(roots): + for enzyme, w in weight[root].items(): + found.setdefault(enzyme, []).append((w, index)) + owner: dict[str, int] = {} + for enzyme, hits in found.items(): + if enzyme in alone or len(hits) > 2: + continue + hits.sort(reverse=True) + if len(hits) == 1 or hits[0][0] > hits[1][0]: + owner[enzyme] = hits[0][1] + shared = set(found) - set(owner) + return ComplexFamilies(complexes=complexes, owner=owner, shared=shared) + + +def level_complex_abundances( + model: "EcModel", + uniprot_ids: Sequence[str], + abundances: Sequence[float], + *, + statistic: str = "median", + fill_missing: bool = True, + families: Optional[ComplexFamilies] = None, +) -> tuple[list[str], np.ndarray, LevelComplexResult]: + """Level molar abundances to one abundance per complex. + + Parameters + ---------- + model + The EcModel whose complexes are used. + uniprot_ids, abundances + Measured molar abundances (any molar unit; per subunit, not per complex). + statistic + ``"median"`` or ``"mean"`` over the measured, non-shared subunits' + abundance divided by their copies. + fill_missing + Give unmeasured, non-shared subunits the complex's abundance. + families + Precomputed ``complex_families(model)``. + + Returns + ------- + ids, abundances, result + The levelled data (unmeasured subunits appended when filled) and + what changed. + """ + if statistic not in _STATISTICS: + raise ValueError(f"statistic must be one of {sorted(_STATISTICS)}, not {statistic!r}") + summarise = _STATISTICS[statistic] + families = families or complex_families(model) + values = {p: float(a) for p, a in zip(uniprot_ids, abundances)} + order = list(uniprot_ids) + result = LevelComplexResult(shared=sorted(families.shared & set(values))) + + for index, complex_ in enumerate(families.complexes): + members = families.owned(index) + measured = [values[e] / c for e, c in members.items() if np.isfinite(values.get(e, np.nan)) and values[e] > 0] + if not measured: + result.skipped.append(tuple(sorted(complex_))) + continue + per_copy = float(summarise(measured)) + for enzyme, copies in members.items(): + level = per_copy * copies + if enzyme in values and np.isfinite(values[enzyme]) and values[enzyme] > 0: + if not np.isclose(values[enzyme], level): + result.adjusted.append(enzyme) + values[enzyme] = level + elif fill_missing: + if enzyme not in values: + order.append(enzyme) + values[enzyme] = level + result.added.append(enzyme) + return order, np.array([values[p] for p in order], dtype=float), result + + +def level_complex_concs( + model: "EcModel", + *, + statistic: str = "median", + fill_missing: bool = True, +) -> LevelComplexResult: + """Give every complex in ``model.ec.concs`` one abundance. + + Run after ``fill_enz_concs`` and before ``constrain_enz_concs``. + Concentrations (mg/gDCW) are converted to molar amounts with + ``ec.mw`` for the levelling and back afterwards. Mutates + ``model.ec.concs``. + + Parameters + ---------- + model + EcModel with ``ec.concs`` filled. + statistic + ``"median"`` (default) or ``"mean"``. + fill_missing + Give unmeasured, non-shared subunits the complex's concentration. + """ + enzymes = list(model.ec.enzymes) + mw = np.asarray(model.ec.mw, dtype=float) + concs = np.asarray(model.ec.concs, dtype=float) + molar = concs / mw + measured = [i for i in range(len(enzymes)) if np.isfinite(molar[i])] + ids, levelled, result = level_complex_abundances( + model, [enzymes[i] for i in measured], molar[measured], + statistic=statistic, fill_missing=fill_missing, + ) + new = concs.copy() + index = {e: i for i, e in enumerate(enzymes)} + for enzyme, value in zip(ids, levelled): + i = index[enzyme] + new[i] = value * mw[i] + model.ec.concs = new + logger.info( + "level_complex_concs: %d subunits levelled, %d filled, %d shared left as measured", + len(result.adjusted), len(result.added), len(result.shared), + ) + return result diff --git a/src/geckopy/limit_proteins/relax_enz_concs.py b/src/geckopy/limit_proteins/relax_enz_concs.py new file mode 100644 index 0000000..be3b8ba --- /dev/null +++ b/src/geckopy/limit_proteins/relax_enz_concs.py @@ -0,0 +1,196 @@ +"""Loosen proteomics constraints by the least total amount, in one LP. + +``flexibilize_enz_concs`` and ``relax_proteomics_greedy`` relax enzyme +caps one at a time: find the most limiting enzyme, raise it, re-solve, +repeat. Which enzymes end up raised then depends on how "most limiting" +is ranked and on the order in which they are found; on a large model +with noisy proteomics, small changes to that ranking can change the +result. ``relax_enz_concs`` asks the question once instead: what is the +smallest set of fold increases of the measured caps (and, if necessary, +of the protein pool) that lets the model reach the target growth rate? + +Every measured enzyme's cap becomes ``usage <= conc * (1 + x)`` and the +sum of the ``x`` is minimised under ``growth >= exp_growth``. The protein +pool may grow by a fold too, at ``pool_weight`` times the price of an +enzyme's fold, so that it grows only where raising enzymes does not +suffice. With ``per_complex=True`` the measured subunits of a complex +(see ``complex_families``) share one ``x``, so concentrations levelled +by ``level_complex_concs`` stay levelled. +""" +from __future__ import annotations + +import logging +from typing import TYPE_CHECKING, Optional + +import numpy as np + +from ..ec_model.constants import POOL_EXCHANGE_ID, USAGE_PREFIX +from .flexibilize_enz_concs import ( + FlexEnzResult, + _finalize_and_sort, + _maybe_add_pool, +) +from .level_complex_concs import complex_families + +if TYPE_CHECKING: + from ..ec_model.ec_model import EcModel + +logger = logging.getLogger(__name__) + + +def relax_enz_concs( + model: "EcModel", + *, + exp_growth: Optional[float] = None, + bio_rxn: Optional[str] = None, + per_complex: bool = False, + weight: str = "relative", + pool_weight: float = 1000.0, + floor: float = 1e-9, +) -> FlexEnzResult: + """Relax enzyme caps by the least total amount that reaches ``exp_growth``. + + Parameters + ---------- + model + EcModel with ``model.ec.concs`` filled and the caps applied + (``fill_enz_concs`` and ``constrain_enz_concs``). The usage + upper bounds (and ``prot_pool_exchange``) are updated in place; + ``ec.concs`` itself is left unchanged, as in + ``flexibilize_enz_concs``. + exp_growth + Target growth rate. Defaults to ``params.gr_exp``. + bio_rxn + Biomass reaction id. Defaults to ``params.bio_rxn``. + per_complex + One fold per complex (its measured, owned subunits; see + ``complex_families``) instead of one per enzyme. + weight + ``"relative"`` minimises the summed fold increases, so an enzyme + measured at almost nothing is not relaxed freely because its + absolute cost is small; ``"absolute"`` minimises the mg/gDCW + added. + pool_weight + Price of one fold of the protein pool relative to one fold of an + enzyme (or, with ``weight="absolute"``, per mg). + floor + Concentration below which a measured value is treated as this + floor when scaling its fold, so a zero measurement can still be + relaxed (at a high relative price). + + Returns + ------- + FlexEnzResult + The relaxed enzymes, as from ``flexibilize_enz_concs`` + (``frequence`` is 1 for every entry), with a ``"prot_pool"`` + entry when the pool was raised. + + Raises + ------ + ValueError + If no enzyme is measured, the weighting is unknown, or no + relaxation reaches ``exp_growth`` (the target is then not limited + by protein). + """ + from optlang.symbolics import add + + from ..adapter import resolve_param + + if weight not in ("relative", "absolute"): + raise ValueError(f"weight must be 'relative' or 'absolute', not {weight!r}") + if exp_growth is None: + exp_growth = resolve_param(model, None, "gr_exp", purpose="relax_enz_concs needs `exp_growth`") + if exp_growth is None: + raise ValueError("exp_growth not provided and params.gr_exp is unset.") + exp_growth = float(exp_growth) + bio = model.reactions.get_by_id( + resolve_param(model, bio_rxn, "bio_rxn", purpose="relax_enz_concs needs the biomass reaction id") + ) + + concs = np.asarray(model.ec.concs, dtype=float) + measured = [i for i in range(len(concs)) if np.isfinite(concs[i])] + if not measured: + raise ValueError("No enzyme concentrations are measured. Run fill_enz_concs and constrain_enz_concs first.") + enzymes = list(model.ec.enzymes) + + old_caps = {i: model.reactions.get_by_id(f"{USAGE_PREFIX}{enzymes[i]}").upper_bound for i in measured} + # folds count from the current cap, so a cap raised before is never lowered + base = {i: min(old_caps[i], 1000.0) for i in measured} + + families = complex_families(model) if per_complex else None + position = {e: i for i, e in enumerate(enzymes)} + measured_set = set(measured) + units: list[list[int]] = [] + assigned: set[int] = set() + for i in measured: + if i in assigned: + continue + if families is None or not base[i] > 0: + unit = [i] + else: # a subunit measured at zero cannot share a fold with the others + unit = [position[m] for m in families.members(enzymes[i]) + if m in position and position[m] in measured_set + and position[m] not in assigned and base[position[m]] > 0] + assigned.update(unit) + units.append(unit) + + pool = model.reactions.get_by_id(POOL_EXCHANGE_ID) + pool_cap = pool.upper_bound + old_bio_bounds = bio.bounds + previous_objective = model.objective.expression + previous_direction = model.objective.direction + + folds, added = [], [] + for k, unit in enumerate(units): + x = model.problem.Variable(f"relax_unit_{k}", lb=0.0, ub=None) + added.append(x) + for i in unit: + rxn = model.reactions.get_by_id(f"{USAGE_PREFIX}{enzymes[i]}") + rxn.upper_bound = 1000.0 + added.append(model.problem.Constraint( + rxn.flux_expression - max(base[i], floor) * x, ub=base[i], name=f"relax_cap_{enzymes[i]}", + )) + folds.append((x, unit)) + pool_x = model.problem.Variable("relax_pool", lb=0.0, ub=None) + pool.upper_bound = 1000.0 + added += [pool_x, model.problem.Constraint(pool.flux_expression - pool_cap * pool_x, ub=pool_cap, name="relax_cap_pool")] + bio.bounds = (exp_growth, max(bio.upper_bound, exp_growth)) + model.add_cons_vars(added) + + try: + terms = [(1.0 if weight == "relative" else sum(max(base[i], floor) for i in unit)) * x for x, unit in folds] + terms.append(pool_weight * (1.0 if weight == "relative" else pool_cap) * pool_x) + model.objective = model.problem.Objective(add(terms), direction="min") + solution = model.optimize() + feasible = solution.status == "optimal" + values = model.solver.primal_values if feasible else {} + finally: + model.remove_cons_vars(added) + model.objective = model.problem.Objective(previous_objective, direction=previous_direction) + bio.bounds = (old_bio_bounds[0], max(old_bio_bounds[1], exp_growth)) + + if not feasible: + for i, cap in old_caps.items(): + model.reactions.get_by_id(f"{USAGE_PREFIX}{enzymes[i]}").upper_bound = cap + pool.upper_bound = pool_cap + raise ValueError( + f"relax_enz_concs: no relaxation of the measured enzymes or the protein pool reaches " + f"growth {exp_growth:g} ({solution.status}); the target is not limited by protein." + ) + + ids, old, new = [], [], [] + for x, unit in folds: + fold = float(values.get(x.name, 0.0)) + for i in unit: + cap = base[i] + max(base[i], floor) * fold + model.reactions.get_by_id(f"{USAGE_PREFIX}{enzymes[i]}").upper_bound = cap if fold > 1e-9 else old_caps[i] + if fold > 1e-9: + ids.append(enzymes[i]); old.append(base[i]); new.append(cap) + pool_fold = float(values.get(pool_x.name, 0.0)) + pool_new = pool_cap * (1 + pool_fold) if pool_fold > 1e-9 else None + pool.upper_bound = pool_new if pool_new is not None else pool_cap + logger.info("relax_enz_concs: %d enzymes relaxed%s", len(ids), + f", protein pool {pool_cap:g} -> {pool_new:g}" if pool_new is not None else "") + result = FlexEnzResult(uniprot_ids=ids, old_concs=np.array(old, dtype=float), + flex_concs=np.array(new, dtype=float), frequence=np.ones(len(ids), dtype=int)) + return _finalize_and_sort(_maybe_add_pool(result, pool_cap if pool_new is not None else None, pool_new)) diff --git a/tests/test_level_complex_concs.py b/tests/test_level_complex_concs.py new file mode 100644 index 0000000..1a05254 --- /dev/null +++ b/tests/test_level_complex_concs.py @@ -0,0 +1,262 @@ +"""Tests for complex_families, level_complex_concs and flexibilize_enz_concs(per_complex=True).""" +from pathlib import Path + +import cobra +import numpy as np +import pytest +from scipy import sparse + +from geckopy import EcModel, ModelAdapter +from geckopy.ec_model.ec_data import EcData +from geckopy.limit_proteins import ( + complex_families, + flexibilize_enz_concs, + level_complex_abundances, + level_complex_concs, +) + +ENZYMES = ["A", "B", "C", "D", "F", "S", "E", "G", "H"] +MW = np.array([100.0, 200.0, 50.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]) + + +def _ec_only(concs, rows): + """An EcModel carrying only ec data: ``rows`` is {rxn: {enzyme: copies}}.""" + model = EcModel("toy") + matrix = np.zeros((len(rows), len(ENZYMES))) + for r, enzymes in enumerate(rows.values()): + for enzyme, copies in enzymes.items(): + matrix[r, ENZYMES.index(enzyme)] = copies + n = len(rows) + model.ec = EcData( + rxns=list(rows), kcat=np.ones(n), source=[""] * n, notes=[""] * n, eccodes=[""] * n, + genes=[f"g_{e}" for e in ENZYMES], enzymes=list(ENZYMES), mw=MW.copy(), + sequence=[""] * len(ENZYMES), concs=np.asarray(concs, dtype=float), + rxn_enz_mat=sparse.csr_matrix(matrix), + ) + return model + + +ROWS = { + "R1_EXP_1": {"A": 1, "B": 2}, + "R1_EXP_2": {"A": 1, "B": 2, "C": 1}, # isozyme variant: same complex + "R2": {"D": 1, "F": 1, "S": 1}, + "R3": {"E": 1, "G": 1, "S": 1}, # S is shared with R2's complex + "R4": {"H": 1}, # a single enzyme: not a complex +} + + +def test_isozyme_variants_are_one_complex_and_a_subunit_in_two_is_shared(): + families = complex_families(_ec_only(np.full(len(ENZYMES), np.nan), ROWS)) + complexes = sorted(sorted(c) for c in families.complexes) + assert complexes == [["A", "B", "C"], ["D", "F", "S"], ["E", "G", "S"]] + assert families.shared == {"S"} + assert sorted(families.members("B")) == ["A", "B", "C"] + assert families.members("S") == ["S"] + assert families.members("H") == ["H"] + + +def test_an_enzyme_that_also_works_alone_is_shared(): + rows = dict(ROWS, R5={"A": 1}) + families = complex_families(_ec_only(np.full(len(ENZYMES), np.nan), rows)) + assert "A" in families.shared + + +def test_levelling_takes_the_median_per_copy_and_scales_by_stoichiometry(): + # molar per copy: A = 1, B = 4 / 2 = 2, C = 6; median 2. + model = _ec_only(np.full(len(ENZYMES), np.nan), ROWS) + ids, values, result = level_complex_abundances(model, ["A", "B", "C"], [1.0, 4.0, 6.0]) + levelled = dict(zip(ids, values)) + assert levelled == pytest.approx({"A": 2.0, "B": 4.0, "C": 2.0}) + assert sorted(result.adjusted) == ["A", "C"] + + +def test_mean_is_available(): + model = _ec_only(np.full(len(ENZYMES), np.nan), ROWS) + ids, values, _ = level_complex_abundances(model, ["A", "B", "C"], [1.0, 4.0, 6.0], statistic="mean") + assert dict(zip(ids, values))["A"] == pytest.approx(3.0) + + +def test_shared_subunits_stay_as_measured_and_missing_ones_are_filled(): + model = _ec_only(np.full(len(ENZYMES), np.nan), ROWS) + ids, values, result = level_complex_abundances(model, ["D", "S"], [1.0, 9.0]) + levelled = dict(zip(ids, values)) + assert levelled["S"] == 9.0 # shared: untouched and not in the summary + assert levelled["D"] == levelled["F"] == 1.0 + assert "F" in result.added + + +def test_filling_can_be_switched_off(): + model = _ec_only(np.full(len(ENZYMES), np.nan), ROWS) + ids, _, _ = level_complex_abundances(model, ["D"], [1.0], fill_missing=False) + assert "F" not in ids + + +def test_level_complex_concs_works_in_mass_units(): + concs = np.full(len(ENZYMES), np.nan) + concs[ENZYMES.index("A")] = 100.0 # 1 mmol per copy at MW 100 + concs[ENZYMES.index("B")] = 1600.0 # 8 mmol, 4 per copy at MW 200 + concs[ENZYMES.index("C")] = 100.0 # 2 mmol per copy at MW 50 + model = _ec_only(concs, ROWS) + level_complex_concs(model) + new = dict(zip(ENZYMES, model.ec.concs)) + # median per copy is 2 mmol: A = 200 mg, B = 2 copies * 2 * 200 = 800 mg, C = 2 * 50 = 100 mg. + assert new["A"] == pytest.approx(200.0) + assert new["B"] == pytest.approx(800.0) + assert new["C"] == pytest.approx(100.0) + assert np.isnan(new["H"]) + + +def test_an_unknown_statistic_is_rejected(): + model = _ec_only(np.full(len(ENZYMES), np.nan), ROWS) + with pytest.raises(ValueError, match="statistic"): + level_complex_abundances(model, ["A"], [1.0], statistic="max") + + +# --------------------------------------------------------------------------- # +# flexibilize_enz_concs(per_complex=True) +# --------------------------------------------------------------------------- # + +def _adapter(tmp_path: Path, gr_exp: float) -> ModelAdapter: + (tmp_path / "model_adapter.toml").write_text( + f'conv_gem = "dummy.xml"\norg_name = "test"\ngr_exp = {gr_exp}\nbio_rxn = "biomass"\n' + ) + return ModelAdapter.from_folder(tmp_path) + + +def _complex_toy(adapter): + """A -> B catalysed by a two-subunit complex X1 + X2; growth equals the flux.""" + model = EcModel("toy", adapter=adapter) + a_e, a_c, b_c = (cobra.Metabolite(m, compartment="c") for m in ("A_e", "A_c", "B_c")) + pool = cobra.Metabolite("prot_pool", compartment="c") + x1, x2 = cobra.Metabolite("prot_X1", compartment="c"), cobra.Metabolite("prot_X2", compartment="c") + model.add_metabolites([a_e, a_c, b_c, pool, x1, x2]) + reactions = { + "EX_A": ({a_e: -1}, (-1000, 0)), + "TR_A": ({a_e: -1, a_c: 1}, (0, 1000)), + "R_AB": ({a_c: -1, b_c: 1, x1: -1, x2: -1}, (0, 1000)), + "biomass": ({b_c: -1}, (0, 1000)), + "prot_pool_exchange": ({pool: 1}, (0, 1000)), + "usage_prot_X1": ({pool: -1, x1: 1}, (0, 1.0)), + "usage_prot_X2": ({pool: -1, x2: 1}, (0, 3.0)), + } + for rid, (mets, bounds) in reactions.items(): + r = cobra.Reaction(rid) + r.add_metabolites(mets) + r.bounds = bounds + model.add_reactions([r]) + model.objective = "biomass" + model.ec = EcData( + rxns=["R_AB"], kcat=np.array([1.0]), source=[""], notes=[""], eccodes=[""], + genes=["g_X1", "g_X2"], enzymes=["X1", "X2"], mw=np.array([100.0, 100.0]), + sequence=["", ""], concs=np.array([1.0, 3.0]), + rxn_enz_mat=sparse.csr_matrix([[1.0, 1.0]]), + ) + return model + + +def test_without_per_complex_only_the_limiting_subunit_is_raised(tmp_path): + model = _complex_toy(_adapter(tmp_path, gr_exp=2.0)) + result = flexibilize_enz_concs(model, verbose=False) + assert result.uniprot_ids == ["X1"] + assert model.reactions.usage_prot_X2.upper_bound == pytest.approx(3.0) + + +def test_per_complex_raises_every_subunit_by_the_same_factor(tmp_path): + model = _complex_toy(_adapter(tmp_path, gr_exp=2.0)) + result = flexibilize_enz_concs(model, verbose=False, per_complex=True) + assert sorted(result.uniprot_ids) == ["X1", "X2"] + factor_1 = model.reactions.usage_prot_X1.upper_bound / 1.0 + factor_2 = model.reactions.usage_prot_X2.upper_bound / 3.0 + assert factor_1 == pytest.approx(factor_2) + assert model.optimize().objective_value == pytest.approx(2.0, rel=5e-3) + + +def test_one_shared_protein_does_not_chain_complexes_together(): + rows = {"R1": {"A": 1, "B": 1}, "R2": {"A": 1, "C": 1}, "R3": {"A": 1, "D": 1, "E": 1, "F": 1}} + families = complex_families(_ec_only(np.full(len(ENZYMES), np.nan), rows)) + assert sorted(sorted(c) for c in families.complexes) == [["A", "B"], ["A", "C"], ["A", "D", "E", "F"]] + assert families.shared == {"A"} + + +def test_an_electron_carrier_in_several_complexes_is_shared_not_a_bridge(): + # B, C, D form a complex that takes carrier S or carrier H (two variants); + # S also pairs with E, like cytochrome c with a dehydrogenase. + rows = {"R1_EXP_1": {"B": 1, "C": 1, "D": 1, "S": 1}, "R1_EXP_2": {"B": 1, "C": 1, "D": 1, "H": 1}, + "R2": {"E": 1, "S": 1}} + families = complex_families(_ec_only(np.full(len(ENZYMES), np.nan), rows)) + assert sorted(sorted(c) for c in families.complexes) == [["B", "C", "D", "H", "S"], ["E", "S"]] + assert "S" in families.shared + assert "E" not in families.members("B") + + +def test_per_complex_relaxes_a_levelled_complex_whose_subunits_bind_together(tmp_path): + model = _complex_toy(_adapter(tmp_path, gr_exp=2.0)) + model.reactions.usage_prot_X2.upper_bound = 1.0 + model.ec.concs = np.array([1.0, 1.0]) + flexibilize_enz_concs(model, verbose=False, per_complex=True) + assert model.reactions.usage_prot_X1.upper_bound == pytest.approx(model.reactions.usage_prot_X2.upper_bound) + assert model.optimize().objective_value == pytest.approx(2.0, rel=5e-3) + + +def test_per_complex_never_lowers_a_cap_below_its_measurement(tmp_path): + model = _complex_toy(_adapter(tmp_path, gr_exp=2.0)) + flexibilize_enz_concs(model, verbose=False, per_complex=True) + assert model.reactions.usage_prot_X1.upper_bound >= 1.0 + assert model.reactions.usage_prot_X2.upper_bound >= 3.0 + + +def test_a_subunit_in_two_complexes_belongs_to_the_one_with_more_reactions(): + # A, B, C, D is a complex on three reactions; A, B with E, F, G, H on one + # (an assembly step): A and B belong to the first, not shared. + rows = {"R1": {"A": 1, "B": 1, "C": 1, "D": 1}, "R2": {"A": 1, "B": 1, "C": 1, "D": 1}, + "R3": {"A": 1, "B": 1, "C": 1, "D": 1}, "R4": {"A": 1, "B": 1, "E": 1, "F": 1, "G": 1, "H": 1}} + families = complex_families(_ec_only(np.full(len(ENZYMES), np.nan), rows)) + assert sorted(families.members("A")) == ["A", "B", "C", "D"] + assert sorted(families.members("E")) == ["E", "F", "G", "H"] + assert families.shared == set() + + +def test_a_subunit_in_two_equally_weighted_complexes_is_shared(): + rows = {"R1": {"A": 1, "B": 1, "C": 1, "D": 1}, "R2": {"A": 1, "B": 1, "E": 1, "F": 1, "G": 1, "H": 1}} + families = complex_families(_ec_only(np.full(len(ENZYMES), np.nan), rows)) + assert {"A", "B"} <= families.shared + + +def _series_toy(adapter): + """A -> B by complex (X1, X2, X3), B -> P by single enzyme E; the pool binds. + + At flux 1 every subunit uses 1 mg of a 1.25 mg cap (80%) and E uses + 1 mg of a 1 mg cap; the pool holds exactly 4 mg. Only E limits. + """ + model = EcModel("series", adapter=adapter) + mets = {n: cobra.Metabolite(n, compartment="c") for n in ("A", "B", "P", "pool", "prot_X1", "prot_X2", "prot_X3", "prot_E")} + model.add_metabolites(list(mets.values())) + spec = { + "EX_A": ({"A": 1}, (0, 1000)), + "R1": ({"A": -1, "B": 1, "prot_X1": -1, "prot_X2": -1, "prot_X3": -1}, (0, 1000)), + "R2": ({"B": -1, "P": 1, "prot_E": -1}, (0, 1000)), + "biomass": ({"P": -1}, (0, 1000)), + "prot_pool_exchange": ({"pool": 1}, (0, 4.0)), + **{f"usage_prot_{e}": ({"pool": -1, f"prot_{e}": 1}, (0, 1.25 if e != "E" else 1.0)) for e in ("X1", "X2", "X3", "E")}, + } + for rid, (st, bounds) in spec.items(): + r = cobra.Reaction(rid) + r.add_metabolites({mets[m]: c for m, c in st.items()}) + r.bounds = bounds + model.add_reactions([r]) + model.objective = "biomass" + model.ec = EcData( + rxns=["R1", "R2"], kcat=np.ones(2), source=[""] * 2, notes=[""] * 2, eccodes=[""] * 2, + genes=["g1", "g2", "g3", "gE"], enzymes=["X1", "X2", "X3", "E"], mw=np.full(4, 100.0), + sequence=[""] * 4, concs=np.array([1.25, 1.25, 1.25, 1.0]), + rxn_enz_mat=sparse.csr_matrix([[1.0, 1.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]), + ) + return model + + +def test_per_complex_does_not_pick_a_large_complex_for_the_pool_price(tmp_path): + model = _series_toy(_adapter(tmp_path, gr_exp=1.1)) + flexibilize_enz_concs(model, verbose=False, per_complex=True) + for e in ("X1", "X2", "X3"): + assert model.reactions.get_by_id(f"usage_prot_{e}").upper_bound == pytest.approx(1.25) + assert model.optimize().objective_value >= 1.1 * 0.9975 - 1e-9 diff --git a/tests/test_relax_enz_concs.py b/tests/test_relax_enz_concs.py new file mode 100644 index 0000000..2b09e8f --- /dev/null +++ b/tests/test_relax_enz_concs.py @@ -0,0 +1,109 @@ +"""Tests for relax_enz_concs: the least total relaxation, in one LP.""" +from pathlib import Path + +import cobra +import numpy as np +import pytest +from scipy import sparse + +from geckopy import EcModel, ModelAdapter +from geckopy.ec_model.ec_data import EcData +from geckopy.limit_proteins import relax_enz_concs + + +def _adapter(tmp_path: Path, gr_exp: float) -> ModelAdapter: + (tmp_path / "model_adapter.toml").write_text( + f'conv_gem = "dummy.xml"\norg_name = "test"\ngr_exp = {gr_exp}\nbio_rxn = "biomass"\n' + ) + return ModelAdapter.from_folder(tmp_path) + + +def _two_routes(adapter, caps, pool=1000.0): + """A -> B by enzyme E1 (2 mg per unit flux) or E2 (1 mg per unit); growth = flux to B. + + ``caps`` gives the measured concentration of each enzyme (mg). + """ + model = EcModel("routes", adapter=adapter) + mets = {n: cobra.Metabolite(n, compartment="c") for n in ("A", "B", "pool", "prot_E1", "prot_E2")} + model.add_metabolites(list(mets.values())) + spec = { + "EX_A": ({"A": 1}, (0, 10)), + "R1": ({"A": -1, "B": 1, "prot_E1": -2}, (0, 1000)), + "R2": ({"A": -1, "B": 1, "prot_E2": -1}, (0, 1000)), + "biomass": ({"B": -1}, (0, 1000)), + "prot_pool_exchange": ({"pool": 1}, (0, pool)), + "usage_prot_E1": ({"pool": -1, "prot_E1": 1}, (0, caps["E1"])), + "usage_prot_E2": ({"pool": -1, "prot_E2": 1}, (0, caps["E2"])), + } + for rid, (st, bounds) in spec.items(): + r = cobra.Reaction(rid) + r.add_metabolites({mets[m]: c for m, c in st.items()}) + r.bounds = bounds + model.add_reactions([r]) + model.objective = "biomass" + model.ec = EcData( + rxns=["R1", "R2"], kcat=np.ones(2), source=[""] * 2, notes=[""] * 2, eccodes=[""] * 2, + genes=["g1", "g2"], enzymes=["E1", "E2"], mw=np.full(2, 100.0), sequence=["", ""], + concs=np.array([caps["E1"], caps["E2"]], dtype=float), rxn_enz_mat=sparse.csr_matrix(np.eye(2)), + ) + return model + + +def test_the_least_summed_fold_increase_reaches_the_target(tmp_path): + # 8 units: x on E1 (2 mg each) and 8 - x on E2 (1 mg each), both capped at 2 mg. + # Folds (x - 1)+ and (3 - x/2)+ sum to 2 + x/2 on [1, 6], 3 - x/2 below: least is 2.5 at x = 1, + # E1 kept at 2 mg and E2 raised to 7 mg. + model = _two_routes(_adapter(tmp_path, gr_exp=8.0), {"E1": 2.0, "E2": 2.0}) + result = relax_enz_concs(model) + assert model.optimize().objective_value == pytest.approx(8.0) + folds = model.reactions.usage_prot_E1.upper_bound / 2.0 - 1 + model.reactions.usage_prot_E2.upper_bound / 2.0 - 1 + assert folds == pytest.approx(2.5) + assert model.reactions.usage_prot_E2.upper_bound == pytest.approx(7.0) + assert "prot_pool" not in result.uniprot_ids + + +def test_enough_capacity_relaxes_nothing(tmp_path): + model = _two_routes(_adapter(tmp_path, gr_exp=3.0), {"E1": 2.0, "E2": 2.0}) + result = relax_enz_concs(model) + assert result.uniprot_ids == [] + assert model.reactions.usage_prot_E1.upper_bound == 2.0 + + +def test_the_pool_grows_only_when_enzymes_cannot_help(tmp_path): + model = _two_routes(_adapter(tmp_path, gr_exp=8.0), {"E1": 100.0, "E2": 100.0}, pool=4.0) + result = relax_enz_concs(model) + assert result.uniprot_ids == ["prot_pool"] + assert model.reactions.prot_pool_exchange.upper_bound == pytest.approx(8.0) # all on E2 + + +def test_a_cap_raised_before_is_not_lowered(tmp_path): + model = _two_routes(_adapter(tmp_path, gr_exp=8.0), {"E1": 2.0, "E2": 2.0}) + model.reactions.usage_prot_E2.upper_bound = 6.0 # raised earlier + relax_enz_concs(model) + assert model.reactions.usage_prot_E2.upper_bound >= 6.0 + assert model.optimize().objective_value == pytest.approx(8.0) + + +def test_an_unreachable_target_is_an_error_and_leaves_the_caps(tmp_path): + model = _two_routes(_adapter(tmp_path, gr_exp=20.0), {"E1": 2.0, "E2": 2.0}) # supply caps growth at 10 + with pytest.raises(ValueError, match="not limited by protein"): + relax_enz_concs(model) + assert model.reactions.usage_prot_E1.upper_bound == 2.0 + assert model.reactions.prot_pool_exchange.upper_bound == 1000.0 + + +def test_per_complex_raises_the_subunits_of_a_complex_together(tmp_path): + from test_level_complex_concs import _adapter as adapter2, _complex_toy + + model = _complex_toy(adapter2(tmp_path, gr_exp=2.0)) + model.ec.concs = np.array([1.0, 1.0]) + model.reactions.usage_prot_X2.upper_bound = 1.0 + relax_enz_concs(model, per_complex=True) + assert model.reactions.usage_prot_X1.upper_bound == pytest.approx(2.0) + assert model.reactions.usage_prot_X2.upper_bound == pytest.approx(2.0) + + +def test_an_unknown_weighting_is_rejected(tmp_path): + model = _two_routes(_adapter(tmp_path, gr_exp=8.0), {"E1": 2.0, "E2": 2.0}) + with pytest.raises(ValueError, match="weight"): + relax_enz_concs(model, weight="max")