Measured zero concentrations and isozyme coverage of the proteomics

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
 

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.                                                                                                                                                                                                                                                          |

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,

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",

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

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",
 ]

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.
 

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.