Let enzymes with an unusable measurement be treated as unmeasured

README.md

@@ -84,6 +84,8 @@ Options change what the models are asked to account for:
   pool (`protein`, the default), or the smallest total flux at the dilution rate
   over every reaction (`flux`) or over the metabolic reactions only
   (`metabolic-flux`). The flux objectives need `--fit-rates`.
+- `--unmeasured UNIPROT ...` treats those enzymes as unmeasured: they are not capped at
+  their measured abundance and draw on the shared protein pool. Default: none.
 - `--scale-protein` rescales biomass to the measured protein content. Off by
   default: the model has no carbon to spare at the measured glucose uptake, and
   CN4 then falls short of its dilution rate.
@@ -147,8 +149,16 @@ the enzyme-constrained solution disposed of surplus carbon through unmeasured
 exits instead of respiring it, so the flux distribution did not describe the
 measured physiology. `--uptake-flex 1.08` gives CN4 the glucose it needs.
 
+`--unmeasured Q12189 P38858` releases RKI1 and SOL3, the enzymes that cap the
+pentose phosphate pathway. RKI1 is measured at about 0.01 mg/gDW at CN38, and the
+model already needs all of it for nucleotide precursors, so the pathway carries no
+flux however it is asked to be used; SOL3 is the next cap. Released, they draw on the
+shared pool, the protein pool does not grow, and the pathway carries flux at every
+condition. Drop the option to build strictly from the measured abundances.
+
 ```bash
-python -m overflow.build --fit-rates --rate-tolerance 0.08 --uptake-flex 1.08
+python -m overflow.build --fit-rates --rate-tolerance 0.08 --uptake-flex 1.08 \
+    --unmeasured Q12189 P38858
 python -m overflow.build_ribosome
 python -m overflow.analyze_usage
 python -m overflow.run_sampling --procs 12

src/overflow/build.py

@@ -13,7 +13,7 @@
 import time
 from dataclasses import dataclass, field
 from pathlib import Path
-from typing import Optional
+from typing import Optional, Sequence
 
 import numpy as np
 import pandas as pd
@@ -101,6 +101,7 @@ def build_condition(
     rate_tolerance: float = 0.05,
     uptake_flex: float = 1.05,
     objective: str = "protein",
+    unmeasured: Sequence[str] = (),
     verbose: bool = True,
 ) -> BuildResult:
     """Build one proteome-constrained ecModel.
@@ -122,6 +123,9 @@ def build_condition(
     protein pool (``"protein"``), or the smallest total flux at the
     dilution rate over every reaction (``"flux"``) or over the metabolic
     reactions only (``"metabolic-flux"``); both need ``fit_rates``.
+
+    ``unmeasured`` lists UniProt IDs whose measurements are not used as
+    caps, so those enzymes draw on the shared protein pool.
     """
     if objective not in OBJECTIVES:
         raise ValueError(f"objective must be one of {OBJECTIVES}, not {objective!r}")
@@ -150,7 +154,7 @@ def build_condition(
 
     proteomics = condition_prot_data(
         condition, model, table=table, masses=masses,
-        fix_complex_subunits=fix_complex_subunits,
+        fix_complex_subunits=fix_complex_subunits, unmeasured=unmeasured,
     )
     result.f_factor = proteomics.f_factor
     result.n_kept = proteomics.filtered.n_kept
@@ -357,6 +361,14 @@ def main(argv: Optional[list[str]] = None) -> int:
              "raise abundances by the least that makes that feasible",
     )
     parser.add_argument("--rate-tolerance", type=float, default=0.05)
+    parser.add_argument(
+        "--unmeasured",
+        nargs="+",
+        default=[],
+        metavar="UNIPROT",
+        help="treat these enzymes as unmeasured: they are not capped at their "
+             "measured abundance and draw on the shared protein pool",
+    )
     parser.add_argument(
         "--objective",
         choices=OBJECTIVES,
@@ -405,7 +417,7 @@ def main(argv: Optional[list[str]] = None) -> int:
             fit_rates=args.fit_rates,
             rate_tolerance=args.rate_tolerance,
             uptake_flex=args.uptake_flex,
-            objective=args.objective,
+            objective=args.objective, unmeasured=args.unmeasured,
         )
         write_outputs(result, args.models_dir, args.results_dir)
         results.append(result)

src/overflow/proteomics.py

@@ -8,7 +8,7 @@
 
 from dataclasses import dataclass
 from pathlib import Path
-from typing import TYPE_CHECKING, Optional
+from typing import TYPE_CHECKING, Optional, Sequence
 
 import numpy as np
 import pandas as pd
@@ -229,6 +229,13 @@ def rescaled_p_tot(
     return condition.p_tot * (kept / total)
 
 
+def drop_proteins(ids, values, drop: Sequence[str]):
+    """``ids`` and ``values`` without the proteins in ``drop``."""
+    drop = set(drop)
+    keep = [i for i, protein in enumerate(ids) if protein not in drop]
+    return [ids[i] for i in keep], np.asarray(values)[keep]
+
+
 @dataclass(frozen=True)
 class ConditionProteomics:
     """Everything the model build needs from one condition's proteomics."""
@@ -247,13 +254,18 @@ def condition_prot_data(
     masses: Optional[dict[str, float]] = None,
     flex_factor: float = 1.96,
     fix_complex_subunits: bool = True,
+    unmeasured: Sequence[str] = (),
 ) -> ConditionProteomics:
     """Assemble the proteomics for one condition.
 
     The f-factor and the protein-content rescaling come from the
     unfiltered replicate means, as they describe the proteome as a
     whole. The abundances that constrain enzymes come from the filtered
     set, optionally after levelling the respiratory complexes.
+
+    Proteins in ``unmeasured`` are left out of the abundances, so the model
+    treats them as unmeasured: they draw on the shared protein pool instead
+    of being capped at a measurement that is too low to carry their flux.
     """
     from geckopy.databases import ProtData
 
@@ -277,6 +289,7 @@ def condition_prot_data(
     if fix_complex_subunits:
         complex_fix = fix_complexes(model, OXPHOS_RXNS, kept_ids, kept_values)
         kept_ids, kept_values = complex_fix.uniprot_ids, complex_fix.abundances
+    kept_ids, kept_values = drop_proteins(kept_ids, kept_values, unmeasured)
 
     return ConditionProteomics(
         prot_data=ProtData(

tests/test_proteomics.py

@@ -11,6 +11,7 @@
 from overflow.config import CONDITION_ORDER, ROOT, load_conditions
 from overflow.proteomics import (
     condition_prot_data,
+    drop_proteins,
     f_factor,
     filter_prot_data,
     mean_abundances,
@@ -210,3 +211,26 @@ def test_f_factor_agrees_with_geckopy(ec_model, table, masses):
     ours = f_factor(ec_model.ec.enzymes, ids, mass)
     theirs = calculate_f_factor(ec_model, ProtData(uniprot_ids=ids, abundances=mass))
     assert ours == pytest.approx(theirs, rel=1e-12)
+
+
+def test_dropping_proteins_keeps_the_others_and_their_values():
+    ids, values = drop_proteins(["A", "B", "C"], np.array([1.0, 2.0, 3.0]), ["B", "Z"])
+    assert ids == ["A", "C"]
+    assert list(values) == [1.0, 3.0]
+
+
+def test_nothing_is_dropped_by_default():
+    ids, values = drop_proteins(["A", "B"], np.array([1.0, 2.0]), ())
+    assert ids == ["A", "B"] and list(values) == [1.0, 2.0]
+
+
+@pytest.mark.slow
+def test_an_unmeasured_enzyme_is_left_out_of_the_abundances(ec_model, table, masses):
+    condition = load_conditions()["CN38"]
+    kwargs = dict(table=table, masses=masses, fix_complex_subunits=False)
+    with_rki1 = condition_prot_data(condition, ec_model, **kwargs)
+    without = condition_prot_data(condition, ec_model, unmeasured=["Q12189"], **kwargs)
+    assert "Q12189" in with_rki1.prot_data.uniprot_ids
+    assert "Q12189" not in without.prot_data.uniprot_ids
+    assert len(without.prot_data.uniprot_ids) == len(with_rki1.prot_data.uniprot_ids) - 1
+    assert without.f_factor == with_rki1.f_factor and without.p_tot == with_rki1.p_tot