diff --git a/CHANGELOG.md b/CHANGELOG.md index fbcd3dc..f6f3242 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,11 @@ ## Unreleased * geckopy is licensed under the GNU General Public License v3 (GPL-3.0), as RAVEN and raven-toolbox. `LICENSE`, the `license` field and the classifier of `pyproject.toml`, and the README state it. Releases up to 4.0.0b2 were published under the MIT licence. +* `load_prot_data` takes `zeros_as_measured`: by default an abundance of 0 is read as not measured (NaN) before any filter, also with `filter_data=False`, so the mean of a protein with a zero replicate leaves the zeros out; with `True` a 0 counts as a measured replicate, is kept in the mean and is exempt from `min_val` and `cut_lowest`, as GECKO's `loadProtData` option `zerosAsMeasured`. +* `get_isozyme_coverage` reports the isozyme groups of `model.ec.rxns` that have both a measured and an unmeasured entry, and the enzymes that `zero_unmeasured_isozymes` can limit, as GECKO's `getIsozymeCoverage`. +* `zero_unmeasured_isozymes` lowers the usage upper bound of the unmeasured isozymes of a reaction that has a measured isozyme, to 0 or to a detection limit (`method="detection_limit"`), as GECKO's `zeroUnmeasuredIsozymes`; it lowers the bounds in place and returns a DataFrame, and a lower bound above the new upper bound is lowered with it. +* `level_complex_concs` and `level_complex_abundances` keep a subunit with a concentration of 0 at 0 instead of filling it with the level of its complex. +* `flexibilize_enz_concs` warns about enzymes with a concentration of 0, which it cannot raise, and does not report them as flexibilized; with `per_complex=True` the count of a complex is recorded on its first subunit with a concentration above 0. ## 4.0.0b2 — 2026-11-02 diff --git a/docs/migrating_from_gecko_matlab.md b/docs/migrating_from_gecko_matlab.md index ae8128d..10ae2f1 100644 --- a/docs/migrating_from_gecko_matlab.md +++ b/docs/migrating_from_gecko_matlab.md @@ -300,7 +300,9 @@ Each row maps a MATLAB function to its geckopy equivalent. The | `fillEnzConcs` | `fill_enz_concs` | — | | `flexibilizeEnzConcs` | `flexibilize_enz_concs` | **`[Py]`** Returns a `FlexEnzResult` dataclass. Raises `NotImplementedError` on gecko-light. | | `getConcControlCoeffs` | `get_conc_control_coeffs` | **`[3→4]`** Probe direction reads the forward `usage_prot_*` upper bound. **`[Py]`** Default path reads the reduced cost of each usage reaction's forward variable (one LP solve plus one check solve per nonzero value, against one solve per enzyme for the finite-difference probe), so an enzyme whose own bound does not bind gets 0 also while the protein pool binds; falls back to finite-difference when the solver exposes no variable duals or returns NaN. Gates on `sol.status == "optimal"`. | +| `getIsozymeCoverage` | `get_isozyme_coverage` | **`[Py]`** Returns the coverage `DataFrame` (columns `rxn`, `n_entries`, `n_measured`, `n_unmeasured`, `measured_enzymes`, `unmeasured_enzymes`, `shared_enzymes`) and the list of zeroable enzymes. Raises `NotImplementedError` on gecko-light. | | `updateProtPool` | *(not ported)* | Obsolete since GECKO 3.2.0. Use `set_prot_pool_size` instead. | +| `zeroUnmeasuredIsozymes` | `zero_unmeasured_isozymes` | **`[Py]`** Lowers the usage upper bounds in place and returns the changed enzymes as a `DataFrame` (`uniprot`, `new_upper_bound`). Method `'detectionLimit'` is `"detection_limit"`. | | — | `relax_proteomics_greedy` *(new)* | **`[Py]`** Alternative to `flexibilize_enz_concs`: relaxes a fixed top-*k* of measured concentrations in a greedy loop until a target growth is reached. Returns a `GreedyRelaxResult`. | ### 5.6 Utilities (`src/geckomat/utilities`) @@ -316,7 +318,7 @@ Each row maps a MATLAB function to its geckopy equivalent. The | `loadConventionalGEM` | `load_conventional_gem` | — | | `loadEcModel` | `load_ec_model` | **`[3→4]`** Auto-flips legacy reverse-direction `usage_prot_*` / `prot_pool_exchange` to the forward convention (warns once). Reads the new cobrapy-style YAML. **`[Py]`** Thin wrapper over `raven_toolbox.io.read_yaml_model` (which owns the YAML schema + typed `EcData`); SBML ecModel I/O dropped in `0.1.0a2`. | | `loadFluxData` | `load_flux_data` | **`[Py]`** Returns a `FluxData` dataclass. | -| `loadProtData` | `load_prot_data` | **`[Py]`** Returns a `ProtData` dataclass. | +| `loadProtData` | `load_prot_data` | **`[Py]`** Returns a `ProtData` dataclass. Option `zerosAsMeasured` is `zeros_as_measured`. | | `mapRxnsToConv` | `map_rxns_to_conv` | **`[Py]`** Returns a `MapRxnsResult` dataclass (was three parallel outputs). | | `plotEcFVA` | *(not ported)* | Use matplotlib / seaborn against the `ec_fva` DataFrame directly. | | `reportEnzymeUsage` | `report_enzyme_usage` | **`[3→4]`** Reads positive `usage_prot_*` flux. **`[Py]`** Returns an `EnzymeUsageReport` dataclass. | diff --git a/docs/relaxation_methods.md b/docs/relaxation_methods.md index b49781e..716e4e3 100644 --- a/docs/relaxation_methods.md +++ b/docs/relaxation_methods.md @@ -92,10 +92,7 @@ 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. + gives every complex one abundance: the median per copy over its subunits with an abundance above 0, scaled back to each subunit by its copies in `ec.rxn_enz_mat`; unmeasured subunits are filled, and a subunit measured at 0 stays 0. `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, diff --git a/src/geckopy/__init__.py b/src/geckopy/__init__.py index a23a0ff..1c7c3d4 100644 --- a/src/geckopy/__init__.py +++ b/src/geckopy/__init__.py @@ -138,10 +138,12 @@ fill_enz_concs, flexibilize_enz_concs, get_conc_control_coeffs, + get_isozyme_coverage, level_complex_abundances, level_complex_concs, relax_enz_concs, relax_proteomics_greedy, + zero_unmeasured_isozymes, ) # utilities @@ -258,10 +260,12 @@ "fill_enz_concs", "flexibilize_enz_concs", "get_conc_control_coeffs", + "get_isozyme_coverage", "level_complex_abundances", "level_complex_concs", "relax_enz_concs", "relax_proteomics_greedy", + "zero_unmeasured_isozymes", # utilities "AddNewRxnsResult", "EnzymeUsageReport", diff --git a/src/geckopy/databases/prot_data_loader.py b/src/geckopy/databases/prot_data_loader.py index e3a33fd..c89b79c 100644 --- a/src/geckopy/databases/prot_data_loader.py +++ b/src/geckopy/databases/prot_data_loader.py @@ -28,6 +28,7 @@ def load_prot_data( max_missing: Union[float, list[float]] = 2.0 / 3.0, cut_lowest: float = 5.0, add_stdevs: float = 1.0, + zeros_as_measured: bool = False, ) -> ProtData: """Load (and optionally filter) absolute proteomics data. @@ -38,6 +39,10 @@ def load_prot_data( condition 0's triplicate, next two are condition 1's duplicate). Levels are in mg/gDCW. + An abundance of 0 is read as not measured (NaN) before any filter, + unless ``zeros_as_measured`` is set. This applies to a file and to a + ``ProtData`` source, with ``filter_data`` True and False. + For each condition, the per-replicate measurements are reduced to a single value per protein: @@ -46,17 +51,23 @@ def load_prot_data( * **maxMissing:** rows with fewer than ``max_missing * n_replicates`` strictly-positive - measurements have all replicates set to NaN. + measurements (with ``zeros_as_measured``, measurements of 0 + or above) have all replicates set to NaN. * **RSD:** rows whose ``std/mean > max_rsd`` are set to NaN. * **Collapse:** per-row ``mean + add_stdevs * std``. - * **minVal:** collapsed values below ``min_val`` -> NaN. + * **minVal:** collapsed values below ``min_val`` -> NaN (with + ``zeros_as_measured``, a collapsed value of 0 is kept). * **Bottom cut:** among non-NaN values, the lowest - ``cut_lowest`` percent are set to NaN. + ``cut_lowest`` percent are set to NaN (with + ``zeros_as_measured``, the values of 0 are left out of the + count and of the cut). 2. ``filter_data=False`` just takes the per-row mean (NaN-safe). Rows that end up NaN in every condition are dropped, along with - their UniProt IDs. + their UniProt IDs. A protein that is not measured in a condition is + NaN in ``abundances`` and leaves the usage of its enzyme drawing from + the protein pool (see ``fill_enz_concs`` and ``constrain_enz_concs``). Ported from GECKO MATLAB: src/geckomat/utilities/loadProtData.m. @@ -71,13 +82,14 @@ def load_prot_data( equal the number of replicate columns in the file. filter_data Whether to apply the filtering pipeline. When False, every - other parameter is unused. + other parameter except ``zeros_as_measured`` is unused. min_val Minimum collapsed value to keep. max_rsd Maximum relative standard deviation per row. max_missing - Minimum fraction of strictly-positive replicates required. + Minimum fraction of strictly-positive replicates required + (replicates of 0 or above with ``zeros_as_measured``). May be a single float (applied to all conditions) or a list of length ``len(repl_per_cond)``. cut_lowest @@ -86,6 +98,20 @@ def load_prot_data( add_stdevs Number of standard deviations to add to the per-row mean when collapsing. + zeros_as_measured + Whether an abundance of 0 is a measurement. If False, a 0 is read + as not measured (NaN) before any filter, also in a ``ProtData`` + source: the mean of a protein over its replicates leaves out the + zeros, and a protein with only zeros is removed, so it draws from + the protein pool. If True, a 0 counts as a measured replicate and + is kept in the mean; a protein whose replicates are all 0 is + returned as 0, is not removed by ``min_val`` or ``cut_lowest``, and + limits the usage of its enzyme to 0 in ``constrain_enz_concs``. + ``max_rsd`` and ``max_missing`` still apply to a protein with some + zero replicates. Proteomics reports a protein below the detection + limit as 0 or as missing. With True, a 0 is read as a protein that + is absent, for example in a deletion strain. Applies with + ``filter_data`` True and False. Returns ------- @@ -105,6 +131,8 @@ def load_prot_data( total_replicates = sum(repl_per_cond) uniprot_ids, abundances = _load_raw(source, total_replicates) + if not zeros_as_measured: + abundances[abundances == 0] = np.nan # Drop empty UniProt IDs. keep = [bool(uid) for uid in uniprot_ids] @@ -128,6 +156,7 @@ def load_prot_data( max_missing=_max_missing_for(max_missing, i), cut_lowest=cut_lowest, add_stdevs=add_stdevs, + zeros_as_measured=zeros_as_measured, ) else: collapsed = _nanmean_axis1(cond_block) @@ -252,6 +281,7 @@ def _filter_condition( max_missing: float, cut_lowest: float, add_stdevs: float, + zeros_as_measured: bool, ) -> np.ndarray: """Run the MATLAB filtering pipeline on one condition's (n_proteins, n_replicates) block. Returns a 1-D collapsed @@ -261,7 +291,8 @@ def _filter_condition( # 1) maxMissing -- only when there's more than one replicate. if n_repl > 1: - positive_count = np.sum(cond > 0, axis=1) + present = cond >= 0 if zeros_as_measured else cond > 0 + positive_count = np.sum(present, axis=1) threshold = max_missing * n_repl too_few = positive_count < threshold cond[too_few, :] = np.nan @@ -283,16 +314,20 @@ def _filter_condition( # 4) minVal filter. too_small = collapsed < min_val + if zeros_as_measured: + too_small &= collapsed != 0 collapsed[too_small] = np.nan # 5) Bottom-cut_lowest% filter. valid = ~np.isnan(collapsed) + if zeros_as_measured: + valid &= collapsed != 0 n_valid = int(valid.sum()) if n_valid > 0 and cut_lowest > 0: cut_n = int(math.floor(n_valid * cut_lowest / 100.0)) if cut_n > 0: valid_indices = np.where(valid)[0] - order = valid_indices[np.argsort(collapsed[valid_indices])] + order = valid_indices[np.argsort(collapsed[valid_indices], kind="stable")] collapsed[order[:cut_n]] = np.nan return collapsed diff --git a/src/geckopy/limit_proteins/__init__.py b/src/geckopy/limit_proteins/__init__.py index 03b96dd..032bcd7 100644 --- a/src/geckopy/limit_proteins/__init__.py +++ b/src/geckopy/limit_proteins/__init__.py @@ -11,6 +11,7 @@ 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 .get_isozyme_coverage import get_isozyme_coverage from .level_complex_concs import ( ComplexFamilies, LevelComplexResult, @@ -24,6 +25,7 @@ RelaxationStep, relax_proteomics_greedy, ) +from .zero_unmeasured_isozymes import zero_unmeasured_isozymes __all__ = [ "ComplexFamilies", @@ -39,8 +41,10 @@ "fill_enz_concs", "flexibilize_enz_concs", "get_conc_control_coeffs", + "get_isozyme_coverage", "level_complex_abundances", "level_complex_concs", "relax_enz_concs", "relax_proteomics_greedy", + "zero_unmeasured_isozymes", ] diff --git a/src/geckopy/limit_proteins/constrain_enz_concs.py b/src/geckopy/limit_proteins/constrain_enz_concs.py index a1e1442..1b543cc 100644 --- a/src/geckopy/limit_proteins/constrain_enz_concs.py +++ b/src/geckopy/limit_proteins/constrain_enz_concs.py @@ -40,6 +40,13 @@ def constrain_enz_concs( ``prot_pool`` coefficient of -1, so a prior lower bound above the new upper bound does not conflict with it. + A NaN entry means that the enzyme was not measured, and its usage + draws from the protein pool. An entry of 0 means that the enzyme was + measured as absent: its usage upper bound is 0, which removes every + flux through the enzyme. ``flexibilize_enz_concs`` multiplies the + concentration and cannot raise a bound of 0; + ``relax_proteomics_greedy`` sets the bound back to its default. + Ported from GECKO MATLAB: src/geckomat/limit_proteins/constrainEnzConcs.m. diff --git a/src/geckopy/limit_proteins/fill_enz_concs.py b/src/geckopy/limit_proteins/fill_enz_concs.py index ad86e17..10bbce8 100644 --- a/src/geckopy/limit_proteins/fill_enz_concs.py +++ b/src/geckopy/limit_proteins/fill_enz_concs.py @@ -40,6 +40,13 @@ def fill_enz_concs( ``model.ec.enzymes`` are silently ignored. Enzymes absent from ``prot_data`` keep their NaN. + A NaN entry in ``model.ec.concs`` means that the protein was not + measured, and the usage of the enzyme draws from the protein pool. An + entry of 0 means that the protein was measured as absent, and + ``constrain_enz_concs`` limits the usage of the enzyme to 0. + ``load_prot_data`` returns a measured 0 with the option + ``zeros_as_measured``, and reads a 0 as NaN otherwise. + For 2-D ``prot_data.abundances`` (multiple conditions / samples), ``data_col`` selects which column to use. diff --git a/src/geckopy/limit_proteins/flexibilize_enz_concs.py b/src/geckopy/limit_proteins/flexibilize_enz_concs.py index 79fb0b2..73b68d2 100644 --- a/src/geckopy/limit_proteins/flexibilize_enz_concs.py +++ b/src/geckopy/limit_proteins/flexibilize_enz_concs.py @@ -130,6 +130,13 @@ def flexibilize_enz_concs( over-relaxing enzymes that the prior iterative pass increased more than necessary. + The usage upper bound of an enzyme is raised as a factor of its + concentration, so an enzyme with a concentration of 0 in + ``model.ec.concs`` stays at 0 and is never flexibilized. If such + enzymes exist, a warning lists them: a growth rate that is not + reached can then be caused by these enzymes, and + ``flexibilize_enz_concs`` reports no limiting enzyme for them. + Ported from GECKO MATLAB: src/geckomat/limit_proteins/flexibilizeEnzConcs.m. @@ -164,7 +171,9 @@ def flexibilize_enz_concs( 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. + ``level_complex_concs`` stay levelled. The ``frequence`` of a + complex is recorded on its first subunit with a concentration above + 0. Off by default. Returns ------- @@ -214,6 +223,19 @@ def flexibilize_enz_concs( "No enzyme concentrations are measured (model.ec.concs is " "all-NaN). Run fill_enz_concs and constrain_enz_concs first." ) + measured_zero = [ + model.ec.enzymes[i] for i in measured_idx if model.ec.concs[i] == 0 + ] + if measured_zero: + listed = ", ".join(measured_zero[:10]) + if len(measured_zero) > 10: + listed += f" and {len(measured_zero) - 10} more" + logger.warning( + "flexibilize_enz_concs: the enzymes %s have a concentration of 0 " + "in model.ec.concs. Their usage upper bound stays 0, because it " + "is raised as a factor of the concentration.", + listed, + ) proteins = [model.ec.enzymes[i] for i in measured_idx] n_proteins = len(proteins) frequence = np.zeros(n_proteins, dtype=int) @@ -296,7 +318,12 @@ def flexibilize_enz_concs( else: unit = int(unit_of[int(np.argmax(scores))]) members = units[unit] - max_idx = members[0] + # The count of a unit is recorded on its first member with a + # concentration above 0 + max_idx = next( + (j for j in members if model.ec.concs[measured_idx[j]] > 0), + members[0], + ) frequence[max_idx] += 1 unit_frequence[unit] += 1 if unit_frequence[unit] > iter_per_enzyme: @@ -406,6 +433,10 @@ def _build_result( flex_mask = frequence > 0 if touched is not None: flex_mask = flex_mask | touched + # An enzyme with a concentration of 0 keeps its bound and is not reported + flex_mask = flex_mask & np.array( + [model.ec.concs[j] != 0 for j in measured_idx], dtype=bool, + ) flex_indices = np.where(flex_mask)[0] if len(flex_indices) == 0: return _maybe_add_pool(FlexEnzResult(), pool_old, pool_new) diff --git a/src/geckopy/limit_proteins/get_isozyme_coverage.py b/src/geckopy/limit_proteins/get_isozyme_coverage.py new file mode 100644 index 0000000..44ef8bf --- /dev/null +++ b/src/geckopy/limit_proteins/get_isozyme_coverage.py @@ -0,0 +1,139 @@ +"""Report isozyme groups with partial proteomics coverage. + +Ported from GECKO MATLAB: +src/geckomat/limit_proteins/getIsozymeCoverage.m. +""" +from __future__ import annotations + +import re +from typing import TYPE_CHECKING + +import numpy as np +import pandas as pd + +if TYPE_CHECKING: + from ..ec_model.ec_model import EcModel + +_COLUMNS = [ + "rxn", "n_entries", "n_measured", "n_unmeasured", + "measured_enzymes", "unmeasured_enzymes", "shared_enzymes", +] + + +def get_isozyme_coverage(model: "EcModel") -> tuple[pd.DataFrame, list[str]]: + """Report isozyme groups with partial proteomics coverage. + + The reactions in ``model.ec.rxns`` that share an identifier once the + suffix ``_EXP_`` is removed are the isozymes of one reaction of the + conventional GEM. An entry is measured when at least one of its enzymes + has a concentration above 0 in ``model.ec.concs``, and unmeasured when + none of its enzymes has a concentration (all NaN). A group has partial + coverage when it has both a measured and an unmeasured entry. An + unmeasured entry next to a measured one can be a protein below the + detection limit as well as a protein that is absent, so the report makes + no statement on which. The kcat of an entry is not considered. A missing + or short ``model.ec.concs`` counts as NaN for the enzymes it does not + cover. + + Ported from GECKO MATLAB: + src/geckomat/limit_proteins/getIsozymeCoverage.m. + + Parameters + ---------- + model + Full (not gecko-light) EcModel with concentrations in + ``model.ec.concs`` (see ``fill_enz_concs``). Not modified. + + Returns + ------- + coverage : pandas.DataFrame + One row per group with partial coverage, sorted by reaction, with + the columns: + + * ``rxn``: the reaction identifier without the ``_EXP_`` suffix. + * ``n_entries``: number of entries of the group in ``model.ec.rxns``. + * ``n_measured``: number of measured entries. + * ``n_unmeasured``: number of unmeasured entries. + * ``measured_enzymes``: list of the enzymes of the measured entries + that have a concentration above 0. + * ``unmeasured_enzymes``: list of the enzymes of the unmeasured + entries. + * ``shared_enzymes``: list of the unmeasured enzymes that also + catalyse an entry outside the unmeasured entries of groups with + partial coverage. ``zero_unmeasured_isozymes`` leaves them + unchanged. + zeroable : list of str + The enzymes that ``zero_unmeasured_isozymes`` acts on, in model + order: the enzymes in ``unmeasured_enzymes`` that are not in + ``shared_enzymes`` of any group. + + Raises + ------ + NotImplementedError + If the model is a gecko-light model. + ValueError + If ``model.ec.rxn_enz_mat`` is missing or empty. + """ + if model.ec.gecko_light: + raise NotImplementedError( + "get_isozyme_coverage does not work on gecko-light models." + ) + matrix = getattr(model.ec, "rxn_enz_mat", None) + if matrix is None or 0 in matrix.shape: + raise ValueError("The model has no ec.rxn_enz_mat. Run make_ec_model first.") + + enzymes = list(model.ec.enzymes) + n_enz = len(enzymes) + concs = np.full(n_enz, np.nan) + given_concs = getattr(model.ec, "concs", None) + given = np.asarray( + given_concs if given_concs is not None else [], dtype=float, + ).ravel() + n_given = min(len(given), n_enz) + concs[:n_given] = given[:n_given] + matrix = (matrix != 0).astype(float).tocsr() + + positive = concs > 0 + missing = np.isnan(concs) + entry_measured = matrix @ positive.astype(float) > 0 + entry_unmeasured = matrix @ (~missing).astype(float) == 0 + base_rxn = [re.sub(r"_EXP_\d+", "", rxn) for rxn in model.ec.rxns] + groups, group_idx = np.unique(np.array(base_rxn, dtype=object), return_inverse=True) + group_idx = np.asarray(group_idx).ravel() + + partial = np.zeros(len(groups), dtype=bool) + for g in range(len(groups)): + in_group = group_idx == g + partial[g] = entry_measured[in_group].any() and entry_unmeasured[in_group].any() + + # An unmeasured entry of a group with partial coverage + in_partial = entry_unmeasured & partial[group_idx] + # An enzyme is zeroable when it is unmeasured and every entry it + # catalyses is such an entry + catalyses = matrix.T.tocsr() + n_catalysed = np.asarray(catalyses.sum(axis=1)).ravel() + n_in_partial = np.asarray(catalyses @ in_partial.astype(float)).ravel() + zeroable_mask = missing & (n_catalysed > 0) & (n_in_partial == n_catalysed) + zeroable = [enzymes[i] for i in np.where(zeroable_mask)[0]] + + rows = np.where(partial)[0] + data: dict[str, list] = {c: [] for c in _COLUMNS} + for g in rows: + in_group = group_idx == g + measured_rows = in_group & entry_measured + unmeasured_rows = in_group & entry_unmeasured + measured_mask = (np.asarray(matrix[measured_rows].sum(axis=0)).ravel() > 0) & positive + unmeasured_mask = (np.asarray(matrix[unmeasured_rows].sum(axis=0)).ravel() > 0) & missing + data["rxn"].append(str(groups[g])) + data["n_entries"].append(int(in_group.sum())) + data["n_measured"].append(int(measured_rows.sum())) + data["n_unmeasured"].append(int(unmeasured_rows.sum())) + data["measured_enzymes"].append([enzymes[i] for i in np.where(measured_mask)[0]]) + data["unmeasured_enzymes"].append([enzymes[i] for i in np.where(unmeasured_mask)[0]]) + data["shared_enzymes"].append( + [enzymes[i] for i in np.where(unmeasured_mask & ~zeroable_mask)[0]] + ) + coverage = pd.DataFrame(data, columns=_COLUMNS) + for column in ("n_entries", "n_measured", "n_unmeasured"): + coverage[column] = coverage[column].astype(int) + return coverage, zeroable diff --git a/src/geckopy/limit_proteins/level_complex_concs.py b/src/geckopy/limit_proteins/level_complex_concs.py index d060f04..9274c35 100644 --- a/src/geckopy/limit_proteins/level_complex_concs.py +++ b/src/geckopy/limit_proteins/level_complex_concs.py @@ -83,7 +83,8 @@ class LevelComplexResult: 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. + Complexes (as sorted subunit tuples) without a non-shared subunit with a + concentration above 0. """ adjusted: list[str] = field(default_factory=list) @@ -192,10 +193,11 @@ def level_complex_abundances( 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. + ``"median"`` or ``"mean"`` over the abundance divided by the + copies of the non-shared subunits with an abundance above 0. fill_missing - Give unmeasured, non-shared subunits the complex's abundance. + Give unmeasured, non-shared subunits the complex's abundance. A + subunit with an abundance of 0 is measured as absent and stays 0. families Precomputed ``complex_families(model)``. @@ -226,7 +228,7 @@ def level_complex_abundances( if not np.isclose(values[enzyme], level): result.adjusted.append(enzyme) values[enzyme] = level - elif fill_missing: + elif fill_missing and values.get(enzyme) != 0: if enzyme not in values: order.append(enzyme) values[enzyme] = level @@ -254,7 +256,14 @@ def level_complex_concs( statistic ``"median"`` (default) or ``"mean"``. fill_missing - Give unmeasured, non-shared subunits the complex's concentration. + Give unmeasured, non-shared subunits the complex's concentration. A + subunit with a concentration of 0 is measured as absent and stays 0. + + Returns + ------- + LevelComplexResult + What changed. ``ec.concs`` has one abundance per complex, except for + subunits at 0, shared subunits and skipped complexes. """ enzymes = list(model.ec.enzymes) mw = np.asarray(model.ec.mw, dtype=float) diff --git a/src/geckopy/limit_proteins/relax_proteomics_greedy.py b/src/geckopy/limit_proteins/relax_proteomics_greedy.py index 424cf08..1c30f80 100644 --- a/src/geckopy/limit_proteins/relax_proteomics_greedy.py +++ b/src/geckopy/limit_proteins/relax_proteomics_greedy.py @@ -63,6 +63,9 @@ def relax_proteomics_greedy( ) -> GreedyRelaxResult: """Greedy relaxation of proteomics constraints until growth >= minimal_growth. + An enzyme with a concentration of 0 in ``model.ec.concs`` is a candidate + like any other and is relaxed in the same way. + Parameters ---------- model diff --git a/src/geckopy/limit_proteins/zero_unmeasured_isozymes.py b/src/geckopy/limit_proteins/zero_unmeasured_isozymes.py new file mode 100644 index 0000000..603a955 --- /dev/null +++ b/src/geckopy/limit_proteins/zero_unmeasured_isozymes.py @@ -0,0 +1,125 @@ +"""Limit unmeasured isozymes that have a measured sibling. + +Ported from GECKO MATLAB: +src/geckomat/limit_proteins/zeroUnmeasuredIsozymes.m. +""" +from __future__ import annotations + +import math +import numbers +from typing import TYPE_CHECKING, Optional + +import numpy as np +import pandas as pd + +from ..ec_model.constants import USAGE_PREFIX +from .get_isozyme_coverage import get_isozyme_coverage + +if TYPE_CHECKING: + from ..ec_model.ec_model import EcModel + +_METHODS = ("zero", "detection_limit") + + +def zero_unmeasured_isozymes( + model: "EcModel", + *, + method: str = "zero", + limit: Optional[float] = None, +) -> pd.DataFrame: + """Limit unmeasured isozymes that have a measured sibling. + + For a reaction with a measured and an unmeasured isozyme (see + ``get_isozyme_coverage``), the usage upper bound of the enzymes of the + unmeasured isozyme is lowered, instead of leaving them to draw from the + protein pool without a limit. The function lowers the usage upper bound + in the model and does not change ``model.ec.concs``. An enzyme that also + catalyses a reaction without a measured sibling is left unchanged, + because limiting it would limit that reaction as well. + + Not detecting a protein is not the same as the protein being absent: a + protein below the detection limit of the measurement is missing from the + data as well. Method ``"zero"`` treats the unmeasured isozyme as absent + and sets its usage upper bound to 0, which removes every flux through + the enzyme and can lower the growth rate or make the model infeasible. + Method ``"detection_limit"`` treats it as present at or below the + detection limit and sets the upper bound to the detection limit. An + upper bound that is already below the value to set is not raised. A + lower bound above the upper bound to set is lowered to it. + + An enzyme lowered by the function has a NaN concentration and is a + candidate in ``relax_proteomics_greedy`` like any other enzyme with a + bound below the default. + + Run after ``constrain_enz_concs``, which resets the usage upper bounds of + all enzymes. The function is not called by other functions of geckopy. + + Ported from GECKO MATLAB: + src/geckomat/limit_proteins/zeroUnmeasuredIsozymes.m (method + ``'detectionLimit'`` there is ``"detection_limit"`` here). + + Parameters + ---------- + model + Full (not gecko-light) EcModel with concentrations in + ``model.ec.concs``, constrained by ``constrain_enz_concs``. The usage + upper bounds are lowered in place. + method + ``"zero"`` or ``"detection_limit"``. + limit + The detection limit in mg/gDCW for ``method="detection_limit"``. + Defaults to the smallest concentration above 0 in + ``model.ec.concs``. + + Returns + ------- + pandas.DataFrame + One row per enzyme that the function acted on, with the columns + ``uniprot`` and ``new_upper_bound``. + + Raises + ------ + ValueError + If ``method`` is unknown, if the detection limit is negative, NaN or + missing (no concentration above 0 and no ``limit``), or if a usage + reaction of an enzyme to act on is missing. + NotImplementedError + If the model is a gecko-light model. + """ + if method not in _METHODS: + raise ValueError(f"method must be 'zero' or 'detection_limit', got {method!r}.") + + _, zeroable = get_isozyme_coverage(model) + + if method == "zero": + new_ub = [0.0] * len(zeroable) + else: + if limit is None: + concs = getattr(model.ec, "concs", None) + concs = np.asarray([] if concs is None else concs, dtype=float) + above = concs[concs > 0] + limit = float(above.min()) if above.size else math.nan + if ( + isinstance(limit, bool) + or not isinstance(limit, numbers.Real) + or math.isnan(limit) + or limit < 0 + ): + raise ValueError("The detection limit must be a nonnegative number.") + new_ub = [float(limit)] * len(zeroable) + + reaction_ids = {r.id for r in model.reactions} + absent = [e for e in zeroable if f"{USAGE_PREFIX}{e}" not in reaction_ids] + if absent: + raise ValueError(f"The model has no usage reaction for {', '.join(absent)}.") + + applied = [] + for enzyme, value in zip(zeroable, new_ub): + rxn = model.reactions.get_by_id(f"{USAGE_PREFIX}{enzyme}") + value = min(value, rxn.upper_bound) + rxn.bounds = (min(rxn.lower_bound, value), value) + applied.append(value) + return pd.DataFrame({ + "uniprot": list(zeroable), + "new_upper_bound": np.array(applied, dtype=float), + }) diff --git a/tests/test_isozyme_coverage.py b/tests/test_isozyme_coverage.py new file mode 100644 index 0000000..6703351 --- /dev/null +++ b/tests/test_isozyme_coverage.py @@ -0,0 +1,311 @@ +"""Tests for get_isozyme_coverage, zero_unmeasured_isozymes and the measured-zero handling of flexibilize_enz_concs.""" +import logging + +import cobra +import numpy as np +import pytest +from scipy import sparse + +from geckopy import EcModel +from geckopy.ec_model.ec_data import EcData +from geckopy.limit_proteins import ( + constrain_enz_concs, + flexibilize_enz_concs, + get_isozyme_coverage, + zero_unmeasured_isozymes, +) + + +# --------------------------------------------------------------------------- # +# Model builders +# --------------------------------------------------------------------------- # + +def _isozyme_model(concs) -> EcModel: + """Five enzymes with usage reactions (upper bound 1000) and the entries + X, X_EXP_1: P1, P2; Y, Y_EXP_1: P3, P4; Z: P4; W, W_EXP_1: P5, P1.""" + enzymes = ["P1", "P2", "P3", "P4", "P5"] + model = EcModel("isozymes") + reactions = [] + for enzyme in enzymes: + rxn = cobra.Reaction(f"usage_prot_{enzyme}") + rxn.bounds = (0.0, 1000.0) + reactions.append(rxn) + model.add_reactions(reactions) + model.ec = EcData( + rxns=["X", "X_EXP_1", "Y", "Y_EXP_1", "Z", "W", "W_EXP_1"], + enzymes=enzymes, + concs=np.array(concs, dtype=float), + rxn_enz_mat=sparse.csr_matrix(np.array([ + [1, 0, 0, 0, 0], + [0, 1, 0, 0, 0], + [0, 0, 1, 0, 0], + [0, 0, 0, 1, 0], + [0, 0, 0, 1, 0], + [0, 0, 0, 0, 1], + [1, 0, 0, 0, 0], + ], dtype=float)), + ) + return model + + +def _usage_ub(model, enzymes): + return [model.reactions.get_by_id(f"usage_prot_{e}").upper_bound for e in enzymes] + + +def _zero_toy(concs) -> EcModel: + """A -> B by E1 or by E2 (R_AB1, R_AB2); growth (R4) equals the flux.""" + model = EcModel("toy") + mets = {n: cobra.Metabolite(n, compartment="c") for n in + ["A_e", "A_c", "B_c", "prot_pool", "prot_E1", "prot_E2"]} + model.add_metabolites(list(mets.values())) + spec = [ + ("EX_A", {"A_e": -1}, -1000, 0), + ("TR_A", {"A_e": -1, "A_c": 1}, 0, 1000), + ("R_AB1", {"A_c": -1, "B_c": 1, "prot_E1": -1}, 0, 1000), + ("R_AB2", {"A_c": -1, "B_c": 1, "prot_E2": -1}, 0, 1000), + ("R4", {"B_c": -1}, 0, 1000), + ("prot_pool_exchange", {"prot_pool": 1}, 0, 1000), + ("usage_prot_E1", {"prot_pool": -1, "prot_E1": 1}, 0, 1000), + ("usage_prot_E2", {"prot_pool": -1, "prot_E2": 1}, 0, 1000), + ] + for rxn_id, stoich, lb, ub in spec: + rxn = cobra.Reaction(rxn_id) + rxn.add_metabolites({mets[k]: float(v) for k, v in stoich.items()}) + rxn.bounds = (float(lb), float(ub)) + model.add_reactions([rxn]) + model.objective = "R4" + model.ec = EcData( + rxns=["R_AB1", "R_AB2"], + kcat=np.array([1.0, 1.0]), + source=["", ""], + notes=["", ""], + eccodes=["", ""], + genes=["g_E1", "g_E2"], + enzymes=["E1", "E2"], + mw=np.array([100.0, 100.0]), + sequence=["", ""], + concs=np.array(concs, dtype=float), + rxn_enz_mat=sparse.csr_matrix(np.eye(2)), + ) + constrain_enz_concs(model) + return model + + +# --------------------------------------------------------------------------- # +# get_isozyme_coverage +# --------------------------------------------------------------------------- # + +def test_coverage_reports_partial_groups(): + """W, X and Y have a measured and an unmeasured entry; Z has one entry + only. P4 also catalyses Z, so zeroing it would limit Z as well.""" + model = _isozyme_model([3, np.nan, 2, np.nan, np.nan]) + coverage, zeroable = get_isozyme_coverage(model) + assert list(coverage.columns) == [ + "rxn", "n_entries", "n_measured", "n_unmeasured", + "measured_enzymes", "unmeasured_enzymes", "shared_enzymes", + ] + assert list(coverage["rxn"]) == ["W", "X", "Y"] + assert list(coverage["n_entries"]) == [2, 2, 2] + assert list(coverage["n_measured"]) == [1, 1, 1] + assert list(coverage["n_unmeasured"]) == [1, 1, 1] + assert list(coverage["measured_enzymes"]) == [["P1"], ["P1"], ["P3"]] + assert list(coverage["unmeasured_enzymes"]) == [["P5"], ["P2"], ["P4"]] + assert list(coverage["shared_enzymes"]) == [[], [], ["P4"]] + assert zeroable == ["P2", "P5"] + + +def test_coverage_counts_a_measured_subunit_as_measured_entry(): + """The entry V has the subunits P1 (measured) and P2 (not measured): it is + a measured entry, and P2 is not an unmeasured isozyme.""" + model = _isozyme_model([3, np.nan, np.nan, np.nan, np.nan]) + model.ec.rxns = ["V", "V_EXP_1"] + model.ec.rxn_enz_mat = sparse.csr_matrix(np.array([ + [1, 1, 0, 0, 0], + [0, 0, 0, 1, 1], + ], dtype=float)) + coverage, zeroable = get_isozyme_coverage(model) + assert list(coverage["rxn"]) == ["V"] + assert list(coverage["measured_enzymes"]) == [["P1"]] + assert list(coverage["unmeasured_enzymes"]) == [["P4", "P5"]] + assert zeroable == ["P4", "P5"] + + +def test_coverage_ignores_measured_zero(): + """A concentration of 0 is neither a measured sibling nor an unmeasured entry.""" + model = _isozyme_model([0, np.nan, 2, np.nan, np.nan]) + coverage, zeroable = get_isozyme_coverage(model) + assert list(coverage["rxn"]) == ["Y"] + assert zeroable == [] + + +def test_coverage_empty_without_measurements(): + model = _isozyme_model([np.nan] * 5) + coverage, zeroable = get_isozyme_coverage(model) + assert len(coverage) == 0 + assert zeroable == [] + model.ec.concs = np.empty(0) + assert len(get_isozyme_coverage(model)[0]) == 0 + + +def test_coverage_short_concs_count_as_nan(): + model = _isozyme_model([3]) + coverage, zeroable = get_isozyme_coverage(model) + full, full_zeroable = get_isozyme_coverage(_isozyme_model([3] + [np.nan] * 4)) + assert coverage.equals(full) + assert zeroable == full_zeroable + + +def test_coverage_keeps_enzyme_of_measured_entry_shared(): + """P2 is not measured, and it is in the measured entry X (with P1) as well + as in the unmeasured entry X_EXP_1, so it is shared and not zeroable.""" + model = _isozyme_model([3, np.nan, 2, np.nan, np.nan]) + dense = model.ec.rxn_enz_mat.toarray() + dense[0, :] = [1, 1, 0, 0, 0] + model.ec.rxn_enz_mat = sparse.csr_matrix(dense) + coverage, zeroable = get_isozyme_coverage(model) + assert coverage.loc[coverage["rxn"] == "X", "shared_enzymes"].iloc[0] == ["P2"] + assert zeroable == ["P5"] + + +def test_coverage_without_concs_is_empty(): + model = _isozyme_model([3, np.nan, 2, np.nan, np.nan]) + del model.ec.concs + coverage, zeroable = get_isozyme_coverage(model) + assert len(coverage) == 0 + assert zeroable == [] + + +def test_coverage_does_not_modify_model(): + model = _isozyme_model([3, np.nan, 2, np.nan, np.nan]) + concs = model.ec.concs.copy() + get_isozyme_coverage(model) + assert np.array_equal(model.ec.concs, concs, equal_nan=True) + assert _usage_ub(model, ["P1", "P2", "P3", "P4", "P5"]) == [1000.0] * 5 + + +def test_coverage_rejects_light_model(): + model = _isozyme_model([3, np.nan, 2, np.nan, np.nan]) + model.ec.gecko_light = True + with pytest.raises(NotImplementedError, match="gecko-light"): + get_isozyme_coverage(model) + + +def test_coverage_rejects_missing_rxn_enz_mat(): + model = _isozyme_model([3, np.nan, 2, np.nan, np.nan]) + model.ec.rxn_enz_mat = sparse.csr_matrix((0, 0)) + with pytest.raises(ValueError, match="rxn_enz_mat"): + get_isozyme_coverage(model) + + +# --------------------------------------------------------------------------- # +# zero_unmeasured_isozymes +# --------------------------------------------------------------------------- # + +def test_zero_sets_usage_bound_to_zero(): + model = _isozyme_model([3, np.nan, 2, np.nan, np.nan]) + concs = model.ec.concs.copy() + changed = zero_unmeasured_isozymes(model) + assert list(changed.columns) == ["uniprot", "new_upper_bound"] + assert list(changed["uniprot"]) == ["P2", "P5"] + assert list(changed["new_upper_bound"]) == [0.0, 0.0] + assert _usage_ub(model, ["P1", "P2", "P3", "P4", "P5"]) == [1000, 0, 1000, 1000, 0] + assert np.array_equal(model.ec.concs, concs, equal_nan=True) + + +def test_detection_limit(): + model = _isozyme_model([3, np.nan, 2, np.nan, np.nan]) + # The default limit is the smallest concentration above 0 + default = _isozyme_model([3, np.nan, 2, np.nan, np.nan]) + zero_unmeasured_isozymes(default, method="detection_limit") + assert _usage_ub(default, ["P2", "P5"]) == [2, 2] + limited = _isozyme_model([3, np.nan, 2, np.nan, np.nan]) + zero_unmeasured_isozymes(limited, method="detection_limit", limit=0.5) + assert _usage_ub(limited, ["P2", "P5", "P4"]) == [0.5, 0.5, 1000] + # An upper bound below the limit is not raised + model.reactions.get_by_id("usage_prot_P2").upper_bound = 0.1 + zero_unmeasured_isozymes(model, method="detection_limit", limit=0.5) + assert _usage_ub(model, ["P2", "P5"]) == [0.1, 0.5] + + +def test_zero_leaves_shared_enzyme_unchanged(): + model = _isozyme_model([0, np.nan, 2, np.nan, np.nan]) + changed = zero_unmeasured_isozymes(model) + # Y is partial but P4 is shared with Z, so nothing is zeroable + assert len(changed) == 0 + assert _usage_ub(model, ["P1", "P2", "P3", "P4", "P5"]) == [1000.0] * 5 + + +def test_zero_rejects_unknown_method_and_bad_limit(): + model = _isozyme_model([3, np.nan, 2, np.nan, np.nan]) + with pytest.raises(ValueError, match="method"): + zero_unmeasured_isozymes(model, method="half") + with pytest.raises(ValueError, match="detection limit"): + zero_unmeasured_isozymes(model, method="detection_limit", limit=-1) + with pytest.raises(ValueError, match="detection limit"): + zero_unmeasured_isozymes(model, method="detection_limit", limit=float("nan")) + no_concs = _isozyme_model([np.nan] * 5) + with pytest.raises(ValueError, match="detection limit"): + zero_unmeasured_isozymes(no_concs, method="detection_limit") + + +def test_zero_rejects_light_model(): + model = _isozyme_model([3, np.nan, 2, np.nan, np.nan]) + model.ec.gecko_light = True + with pytest.raises(NotImplementedError, match="gecko-light"): + zero_unmeasured_isozymes(model) + + +def test_zero_requires_usage_reactions(): + model = _isozyme_model([3, np.nan, 2, np.nan, np.nan]) + # P2 and P5 are zeroable; the reaction of the last one is missing + model.remove_reactions(["usage_prot_P5"]) + with pytest.raises(ValueError, match="usage reaction"): + zero_unmeasured_isozymes(model) + assert _usage_ub(model, ["P1", "P2", "P3", "P4"]) == [1000.0] * 4 + + +# --------------------------------------------------------------------------- # +# flexibilize_enz_concs with a measured zero +# --------------------------------------------------------------------------- # + +def test_flexibilize_leaves_measured_zero_at_zero(caplog): + """E1 is measured at 1 and limits growth; E2 is measured at 0 on a + parallel route. E1 is raised, E2 stays at 0 and the function warns.""" + # An enzyme at 0 is excluded from the result by the mask in `_build_result` (see test_level_complex_concs.py for the per-complex case). + model = _zero_toy([1, 0]) + with caplog.at_level(logging.WARNING): + result = flexibilize_enz_concs(model, exp_growth=2.0, bio_rxn="R4", verbose=False) + assert "enzymes E2 have a concentration of 0" in caplog.text + assert result.uniprot_ids == ["E1"] + assert _usage_ub(model, ["E2"]) == [0.0] + assert model.optimize().objective_value == pytest.approx(2.0, rel=5e-3) + + +def test_flexibilize_cannot_raise_measured_zero(): + """The only route is measured at 0: no enzyme is limiting, the bound stays 0.""" + # An enzyme at 0 is excluded from the result by the mask in `_build_result` (see test_level_complex_concs.py for the per-complex case). + model = _zero_toy([np.nan, 0]) + model.reactions.get_by_id("R_AB1").upper_bound = 0.0 + result = flexibilize_enz_concs(model, exp_growth=2.0, bio_rxn="R4", verbose=False) + assert result.uniprot_ids == [] + assert _usage_ub(model, ["E2"]) == [0.0] + assert model.optimize().objective_value == pytest.approx(0.0, abs=1e-9) + + +def test_flexibilize_warning_lists_at_most_ten_enzymes(caplog): + model = _zero_toy([1, 0]) + names = [f"Q{i}" for i in range(12)] + model.ec.enzymes = ["E1"] + names + model.ec.concs = np.array([1.0] + [0.0] * 12) + model.ec.mw = np.ones(13) + with caplog.at_level(logging.WARNING): + flexibilize_enz_concs(model, exp_growth=0.5, bio_rxn="R4", verbose=False) + assert "Q9" in caplog.text and "Q10" not in caplog.text + assert "and 2 more" in caplog.text + + +def test_flexibilize_does_not_warn_without_measured_zero(caplog): + model = _zero_toy([1, 1]) + with caplog.at_level(logging.WARNING): + flexibilize_enz_concs(model, exp_growth=2.0, bio_rxn="R4", verbose=False) + assert "concentration of 0" not in caplog.text diff --git a/tests/test_level_complex_concs.py b/tests/test_level_complex_concs.py index 1a05254..0b1a06e 100644 --- a/tests/test_level_complex_concs.py +++ b/tests/test_level_complex_concs.py @@ -85,6 +85,20 @@ def test_shared_subunits_stay_as_measured_and_missing_ones_are_filled(): assert "F" in result.added +def test_a_measured_zero_is_not_filled_with_the_complex_level(): + """B is measured at 0 in a complex with A and C: it stays 0.""" + model = EcModel("toy") + model.ec = EcData( + rxns=["R1"], kcat=np.ones(1), source=[""], notes=[""], eccodes=[""], + genes=["g_A", "g_B", "g_C"], enzymes=["A", "B", "C"], mw=np.full(3, 100.0), + sequence=[""] * 3, concs=np.array([2.0, 0.0, 4.0]), + rxn_enz_mat=sparse.csr_matrix(np.ones((1, 3))), + ) + result = level_complex_concs(model) + np.testing.assert_allclose(model.ec.concs, [3.0, 0.0, 3.0], atol=1e-12) + assert 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) @@ -260,3 +274,53 @@ def test_per_complex_does_not_pick_a_large_complex_for_the_pool_price(tmp_path): 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 + + +def test_per_complex_does_not_report_a_measured_zero_subunit(tmp_path): + """X0 is measured at 0 in the complex of X1 and X2: it keeps its bound and + is not reported, and no ratio is NaN.""" + model = _complex_toy(_adapter(tmp_path, gr_exp=2.0)) + pool = model.metabolites.get_by_id("prot_pool") + x0 = cobra.Metabolite("prot_X0", compartment="c") + usage = cobra.Reaction("usage_prot_X0") + usage.add_metabolites({pool: -1.0, x0: 1.0}) + usage.bounds = (0.0, 0.0) + model.add_reactions([usage]) + model.ec.enzymes = ["X1", "X2", "X0"] + model.ec.genes = ["g_X1", "g_X2", "g_X0"] + model.ec.mw = np.array([100.0, 100.0, 100.0]) + model.ec.sequence = ["", "", ""] + model.ec.concs = np.array([1.0, 3.0, 0.0]) + model.ec.rxn_enz_mat = sparse.csr_matrix([[1.0, 1.0, 1.0]]) + result = flexibilize_enz_concs(model, per_complex=True, verbose=False) + assert "X0" not in result.uniprot_ids + assert not np.isnan(result.ratio_incr).any() + assert model.reactions.usage_prot_X0.upper_bound == 0.0 + + +def test_a_complex_with_only_zero_subunits_is_skipped(): + model = EcModel("toy") + model.ec = EcData( + rxns=["R1"], kcat=np.ones(1), source=[""], notes=[""], eccodes=[""], + genes=["g_A", "g_B", "g_C"], enzymes=["A", "B", "C"], mw=np.full(3, 100.0), + sequence=[""] * 3, concs=np.array([0.0, np.nan, np.nan]), + rxn_enz_mat=sparse.csr_matrix(np.ones((1, 3))), + ) + result = level_complex_concs(model) + np.testing.assert_array_equal(model.ec.concs, [0.0, np.nan, np.nan]) + assert len(result.skipped) == 1 + assert result.added == [] + + +def test_per_complex_records_the_count_on_the_first_subunit_above_zero(tmp_path): + """X1 is measured at 0 with an open usage bound, X2 limits the complex. X2 + is raised and carries the count; X1 is not reported.""" + model = _complex_toy(_adapter(tmp_path, gr_exp=2.0)) + model.ec.concs = np.array([0.0, 1.0]) + model.reactions.usage_prot_X1.upper_bound = 1000.0 + model.reactions.usage_prot_X2.upper_bound = 1.0 + result = flexibilize_enz_concs(model, per_complex=True, verbose=False) + assert result.uniprot_ids == ["X2"] + assert list(result.frequence) == [1] + assert not np.isnan(result.ratio_incr).any() + assert model.reactions.usage_prot_X1.upper_bound == 1000.0 diff --git a/tests/test_load_prot_data.py b/tests/test_load_prot_data.py index 54bbb9d..e55d7b0 100644 --- a/tests/test_load_prot_data.py +++ b/tests/test_load_prot_data.py @@ -4,7 +4,7 @@ import numpy as np import pytest -from geckopy.databases import load_prot_data +from geckopy.databases import ProtData, load_prot_data # --------------------------------------------------------------------------- # @@ -274,3 +274,125 @@ def test_loaded_data_can_feed_calculate_f_factor(tmp_path): f = calculate_f_factor(model, pd) # P1 in model contributes 10; total = 15. f = 10/15 ~= 0.667. assert f == pytest.approx(10.0 / 15.0) + + +# --------------------------------------------------------------------------- # +# zeros_as_measured +# --------------------------------------------------------------------------- # + +_ZERO_ROWS = [ + ["uniprot", "r1", "r2", "r3"], + ["Z1", "0", "0", "0"], + ["Z2", "0", "4", "6"], + ["Z3", "5", "5", "5"], +] + + +@pytest.mark.parametrize("filter_data", [True, False]) +def test_zero_is_not_measured_by_default(tmp_path, filter_data): + """Z1 is 0 in every replicate; Z2 is 0 in one replicate. Both are read as + not measured: Z1 is dropped, and the mean of Z2 is over the two values.""" + p = tmp_path / "z.tsv" + _write_tsv(p, _ZERO_ROWS) + pd = load_prot_data( + p, [3], filter_data=filter_data, add_stdevs=0.0, cut_lowest=0.0, + ) + assert pd.uniprot_ids == ["Z2", "Z3"] + np.testing.assert_allclose(pd.abundances, [[5.0], [5.0]], atol=1e-12) + + +@pytest.mark.parametrize("filter_data", [True, False]) +def test_zeros_as_measured_keeps_zero(tmp_path, filter_data): + """With zeros_as_measured Z1 stays 0 and the 0 of Z2 enters its mean.""" + p = tmp_path / "z.tsv" + _write_tsv(p, _ZERO_ROWS) + pd = load_prot_data( + p, [3], filter_data=filter_data, add_stdevs=0.0, cut_lowest=0.0, + zeros_as_measured=True, + ) + assert pd.uniprot_ids == ["Z1", "Z2", "Z3"] + np.testing.assert_allclose(pd.abundances, [[0.0], [10.0 / 3.0], [5.0]], atol=1e-12) + + +def test_measured_zero_is_exempt_from_min_val_and_cut_lowest(tmp_path): + """P1 is 0; P2 to P10 are 2 to 10. min_val 2.5 removes P2, then + cut_lowest 20 percent removes one of the eight positive values (P3) and + not P1.""" + p = tmp_path / "z.tsv" + rows = [["uniprot", "r1", "r2", "r3"], ["P1", "0", "0", "0"]] + rows += [[f"P{k}", str(float(k)), str(float(k)), str(float(k))] for k in range(2, 11)] + _write_tsv(p, rows) + pd = load_prot_data( + p, [3], add_stdevs=0.0, min_val=2.5, cut_lowest=20.0, zeros_as_measured=True, + ) + assert pd.uniprot_ids == ["P1"] + [f"P{k}" for k in range(4, 11)] + np.testing.assert_allclose( + pd.abundances[:, 0], [0.0] + [float(k) for k in range(4, 11)], atol=1e-12, + ) + default = load_prot_data( + p, [3], add_stdevs=0.0, min_val=2.5, cut_lowest=20.0, + ) + assert default.uniprot_ids == [f"P{k}" for k in range(4, 11)] + + +def test_prot_data_source_reads_zero_as_not_measured_by_default(): + source = ProtData( + uniprot_ids=["Z1", "Z2", "Z3"], + abundances=np.array([[0.0, 0.0, 0.0], [0.0, 4.0, 6.0], [5.0, 5.0, 5.0]]), + ) + original = source.abundances.copy() + pd = load_prot_data(source, [3], add_stdevs=0.0, cut_lowest=0.0) + assert pd.uniprot_ids == ["Z2", "Z3"] + np.testing.assert_allclose(pd.abundances, [[5.0], [5.0]], atol=1e-12) + pd = load_prot_data(source, [3], add_stdevs=0.0, cut_lowest=0.0, zeros_as_measured=True) + assert pd.uniprot_ids == ["Z1", "Z2", "Z3"] + np.testing.assert_allclose(pd.abundances, [[0.0], [10.0 / 3.0], [5.0]], atol=1e-12) + np.testing.assert_array_equal(source.abundances, original) + + +def test_max_missing_counts_a_zero_replicate_only_with_zeros_as_measured(tmp_path): + """A = [0, 5, 5] has two positive replicates of three, so max_missing 1.0 + removes it by default; with zeros_as_measured all three are present, the + RSD is 2.887 / 3.333 = 0.87 and the mean is 10/3.""" + p = tmp_path / "m.tsv" + _write_tsv(p, [ + ["uniprot", "r1", "r2", "r3"], + ["A", "0", "5", "5"], + ["B", "5", "5", "5"], + ]) + kwargs = dict(max_missing=1.0, add_stdevs=0.0, cut_lowest=0.0) + pd = load_prot_data(p, [3], **kwargs) + assert pd.uniprot_ids == ["B"] + pd = load_prot_data(p, [3], zeros_as_measured=True, **kwargs) + assert pd.uniprot_ids == ["A", "B"] + np.testing.assert_allclose(pd.abundances, [[10.0 / 3.0], [5.0]], atol=1e-12) + + +def test_zero_gives_nan_in_one_condition_and_a_value_in_the_other(tmp_path): + p = tmp_path / "c.tsv" + _write_tsv(p, [ + ["uniprot", "r1", "r2", "r3", "r4"], + ["P", "0", "0", "3", "5"], + ["Q", "2", "2", "2", "2"], + ]) + kwargs = dict(add_stdevs=0.0, cut_lowest=0.0) + pd = load_prot_data(p, [2, 2], **kwargs) + assert pd.uniprot_ids == ["P", "Q"] + np.testing.assert_allclose(pd.abundances, [[np.nan, 4.0], [2.0, 2.0]], atol=1e-12) + pd = load_prot_data(p, [2, 2], zeros_as_measured=True, **kwargs) + np.testing.assert_allclose(pd.abundances, [[0.0, 4.0], [2.0, 2.0]], atol=1e-12) + + +@pytest.mark.parametrize("zeros_as_measured", [False, True]) +def test_cut_lowest_removes_the_first_of_tied_values(tmp_path, zeros_as_measured): + """Ten values 1,1,2,2,...,5,5; cut_lowest 10 percent removes one value, + the first of the two lowest (P1).""" + p = tmp_path / "t.tsv" + _write_tsv(p, [ + ["uniprot", "r1"], + *[[f"P{i}", str(float(v))] for i, v in enumerate([1, 1, 2, 2, 3, 3, 4, 4, 5, 5], 1)], + ]) + pd = load_prot_data( + p, [1], add_stdevs=0.0, cut_lowest=10.0, zeros_as_measured=zeros_as_measured, + ) + assert pd.uniprot_ids == [f"P{i}" for i in range(2, 11)]