diff --git a/.github/workflows/parity-nightly.yml b/.github/workflows/parity-nightly.yml index e15af307..49255d99 100644 --- a/.github/workflows/parity-nightly.yml +++ b/.github/workflows/parity-nightly.yml @@ -57,7 +57,7 @@ jobs: cache: pip cache-dependency-path: pyproject.toml - - run: pip install --upgrade pip + - run: python -m pip install --upgrade pip - run: pip install -e ".[dev]" gurobipy # Fail here, loudly, rather than let the tests skip and the job go green diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index 78a25284..8cbfbf77 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -31,7 +31,7 @@ jobs: python-version: "3.12" cache: pip cache-dependency-path: pyproject.toml - - run: pip install --upgrade pip + - run: python -m pip install --upgrade pip - run: pip install build - run: python -m build diff --git a/CHANGELOG.md b/CHANGELOG.md index 220d2dba..b629dd22 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,18 @@ Milestones in the raven-toolbox port. For function-level status see [docs/raven_migration.md](https://github.com/CompMetEng/raven-toolbox/blob/develop/docs/reference/migration.md); for open work see [docs/todo.md](https://github.com/CompMetEng/raven-toolbox/blob/develop/docs/reference/todo.md). +## Unreleased + +* `init.ftinit` and `init.run_ftinit` give each scored reaction its own big-M, its flux + capacity `max(min(ub, big_m), min(-lb, big_m), force_on)` with an infinite bound counted as + `big_m`, as RAVEN's `ftINITInternalAlg` does. `big_m` is the cap on that capacity; its + default is 1000 (was 100). +* `tasks.find_task_essential_reactions` (and so `init.prep_init_model`) re-solves a knock-out + from a cold start when the warm-started solve does not report a feasible status, and + counts the reaction as essential only if that solve is not feasible either, as RAVEN's + `solveLP` does. Gurobi could report a feasible knock-out as infeasible when started from + the previous knock-out's basis, which added reactions to the essential set. + ## 3.0.0b1 — 2026-09-07 First beta toward 3.0, matching the upcoming MATLAB RAVEN 3.0 beta. diff --git a/scripts/analyze_init_params.py b/scripts/analyze_init_params.py index 0dfba826..0f4436c0 100644 --- a/scripts/analyze_init_params.py +++ b/scripts/analyze_init_params.py @@ -20,7 +20,7 @@ * ``ftinit_milp`` — single staged-MILP step (step 0 of series '1+1') on the merged model. Cheap (~30-200 s each); the core sweep for ``mip_gap``/``big_m``/``force_on``. * ``prep_scale`` — rescaleModelForINIT on/off and its ``max_stoich_diff``, fed into the - same step-0 MILP. Shows why scaling is needed for a fixed big-M. + same step-0 MILP. Shows why scaling is needed when the big-M is capped. * ``ftinit_full`` — the whole ``ftinit`` pipeline (both steps + gap-fill). Sweeps ``mip_gap``/``big_m``. Expensive (~200 s+/config). @@ -185,7 +185,7 @@ def cfg(label, **kw): save() return store[key] - # mip_gap sweep (big_m=100, force_on=0.1) + # mip_gap sweep (force_on=0.1; big_m=100 caps the per-reaction big-M) res = [cfg(f"gap={g}", mip_gap=g, big_m=100.0, force_on=0.1, time_limit=900) for g in MIP_GAPS] _finalize(res) doc += _print_table("ftINIT step-0: mip_gap (big_m=100, force_on=0.1)", res) @@ -194,9 +194,10 @@ def cfg(label, **kw): res = [cfg(f"big_m={int(b)}", mip_gap=0.001, big_m=b, force_on=0.1, time_limit=900) for b in BIG_MS] _finalize(res) doc += _print_table("ftINIT step-0: big_m (gap=0.001, force_on=0.1)", res, - "big_m caps a scored reaction's flux; large values weaken the LP relaxation.") + "big_m caps the per-reaction big-M (a scored reaction's flux capacity); " + "large values weaken the LP relaxation.") - # force_on sweep (gap=0.001, big_m=100) — changes the model (connectivity threshold) + # force_on sweep (gap=0.001, big_m=100 as the big-M cap); changes the model (connectivity threshold) res = [cfg(f"force_on={fo}", mip_gap=0.001, big_m=100.0, force_on=fo, time_limit=900) for fo in FORCE_ONS] _finalize(res) doc += _print_table("ftINIT step-0: force_on (gap=0.001, big_m=100)", res, @@ -236,7 +237,7 @@ def cfg(label, msd): res.append(cfg("scale=off", None)) _finalize(res) doc += _print_table("prep scaling: rescaleModelForINIT max_stoich_diff (+scale off), big_m=100", res, - "With big_m=100 fixed, scale=off / poor conditioning is expected to be " + "With the big-M capped at big_m=100, scale=off / poor conditioning is expected to be " "infeasible or far slower — that is the reason scaling is on by default.") return doc @@ -250,7 +251,7 @@ def cfg(label, **kw): print(f"[ftinit_full] {label} ...", flush=True) t = time.time() try: - out = ftinit(prep, rxn_scores, gene_scores=gene_scores, series="1+1", **kw) + out = ftinit(prep, rxn_scores=rxn_scores, gene_scores=gene_scores, series="1+1", **kw) r = Result(label=label, seconds=time.time() - t, status="ok", objective=0.0, n_kept=len(out.reactions), reactions=sorted(x.id for x in out.reactions)) diff --git a/scripts/analyze_init_robustness.py b/scripts/analyze_init_robustness.py index 6b982bbd..9f9ef8d5 100644 --- a/scripts/analyze_init_robustness.py +++ b/scripts/analyze_init_robustness.py @@ -188,7 +188,7 @@ def cached(key, fn): def model_for(e, **kw): g = gene_scores_from_expression(e, 1.0) r = score_reactions_from_genes(ref, g, no_gene_score=kw.get("no_gene_score", -2.0)) - return ftinit(prep, r, gene_scores=g, series="1+1", + return ftinit(prep, rxn_scores=r, gene_scores=g, series="1+1", force_on=kw.get("force_on", 0.1), mip_gap=MIP_GAP, time_limit=TIME_LIMIT) phases = set(args.phase.split(",")) diff --git a/scripts/analyze_init_solvers.py b/scripts/analyze_init_solvers.py index 84f27fb3..268437bd 100644 --- a/scripts/analyze_init_solvers.py +++ b/scripts/analyze_init_solvers.py @@ -101,7 +101,7 @@ def run(solver: str) -> Result: local_prep.min_model.solver = solver g = gene_scores_from_expression(expr, 1.0) r = score_reactions_from_genes(ref, g) - model = ftinit(local_prep, r, gene_scores=g, series="1+1", + model = ftinit(local_prep, rxn_scores=r, gene_scores=g, series="1+1", mip_gap=args.mip_gap, time_limit=args.time_limit) rset = sorted(x.id for x in model.reactions) res = Result(solver, time.time() - t, "ok", len(rset), rset) diff --git a/src/raven_toolbox/analysis/fseof.py b/src/raven_toolbox/analysis/fseof.py index 23c2d86b..b3e8f605 100644 --- a/src/raven_toolbox/analysis/fseof.py +++ b/src/raven_toolbox/analysis/fseof.py @@ -1,21 +1,19 @@ -"""Flux Scanning based on Enforced Objective Flux — FSEOF (port + redesign). +"""Flux Scanning based on Enforced Objective Flux (FSEOF). FSEOF (Choi et al., Appl Environ Microbiol 2010) finds metabolic-engineering targets for over-producing a metabolite: enforce an increasing flux toward the target product -while optimising growth, and watch how each reaction's flux responds. This is a port -of RAVEN's ``FSEOF`` with a substantially richer, more robust output (RAVEN's -weaknesses are noted in IMPROVEMENTS, FS1–FS4): - -* **Robust trend, not strict monotonicity.** Each reaction's flux is regressed against - the enforced product flux across the scan; the **slope** is the response and the - **correlation** (|r|) is a quality score. A reaction is a target if it tracks the - product cleanly (|r| ≥ ``correlation_threshold``), so one noisy step from LP - alternative optima does not discard it (and pFBA per step keeps the scan stable). -* **Direction classification RAVEN lacks.** Targets are labelled ``amplify`` (|flux| - rises with the product → over-express), ``knockdown`` (|flux| falls), or ``knockout`` - (|flux| → ~0 → delete). RAVEN only ever reports the amplification targets. -* **Gene-level view** via :attr:`FSEOFResult.gene_targets`, and the full flux scan is - retained in :attr:`FSEOFResult.scan` — all as DataFrames, not a printed TSV. +while optimising growth, and watch how each reaction's flux responds. This follows +RAVEN's ``FSEOF`` and returns tables instead of a printed list: + +* **Trend of the absolute flux.** Each reaction's absolute flux is regressed against the + enforced product flux across the scan; the **slope** is the response and the + **correlation** (|r|) is a quality score. A reaction is a target when its track is not + constant and |r| >= ``correlation_threshold``. Every step is solved with pFBA. +* **Direction classification.** Targets are labelled ``amplify`` (the absolute flux + rises with the product), ``knockdown`` (it falls) or ``knockout`` (it falls to at most + 1e-8), like the Amplify, Knockdown and Knockout labels of RAVEN's ``FSEOF``. +* **Gene-level view** via :attr:`FSEOFResult.gene_targets`, and the full flux scan in + :attr:`FSEOFResult.scan`, both as DataFrames. """ from __future__ import annotations @@ -97,7 +95,17 @@ def fseof( product flux for it to be reported as a target; raise it to keep only the most cleanly-trending reactions. flux_eps: - Tolerance below which a flux or slope is treated as zero. + Tolerance below which the theoretical maximum of the target flux is treated as zero + (the scan then raises). Reaction fluxes are classified as in RAVEN's ``FSEOF``: a + reaction is a target when some absolute flux exceeds 1e-10, the standard deviation of + its absolute flux is above 0 and the correlation reaches ``correlation_threshold``; + it is ``knockout`` when the slope is negative and the last absolute flux is at most + 1e-8, ``knockdown`` for a negative slope otherwise. + + Differences from ``FSEOF.m``: the scan stops at the first infeasible enforced level, + so the last absolute flux used for the knockout label is that of the last feasible + level; ``targets`` is sorted by ``|slope| * |r|`` and not grouped by label; and + ``flux_eps`` only gates the check of the theoretical maximum. min_target: Floor of the enforced-flux range, evenly spaced up to the scan's ceiling (``max_fraction`` of the theoretical maximum) in ``n_steps`` steps. @@ -127,40 +135,39 @@ def fseof( try: columns[level] = pfba(model).fluxes except OptimizationError: - break # enforced flux became infeasible — stop scanning + break # enforced flux became infeasible: stop scanning enforced.append(level) if len(enforced) < 2: raise RuntimeError("FSEOF needs at least two feasible enforced-flux levels.") scan = pd.DataFrame(columns) - targets = _classify(model, scan, np.asarray(enforced), correlation_threshold, flux_eps) + targets = _classify(model, scan, np.asarray(enforced), correlation_threshold) return FSEOFResult(scan=scan, enforced=enforced, targets=targets) -def _classify(model, scan, enforced, corr_threshold, flux_eps) -> pd.DataFrame: +_TRACK_EPS = 1e-10 # a track must exceed this somewhere to be fitted, as in FSEOF.m +_KNOCKOUT_EPS = 1e-8 # last absolute flux at or below this is a knockout, as in FSEOF.m + + +def _classify(model, scan, enforced, corr_threshold) -> pd.DataFrame: rows = [] for rxn in model.reactions: flux = scan.loc[rxn.id, enforced.tolist() if hasattr(enforced, "tolist") else enforced] flux = flux.to_numpy(dtype=float) initial, final = flux[0], flux[-1] - if flux.std() < flux_eps: # flat -> no response + # Slope and correlation are taken on absolute flux, as in RAVEN's FSEOF, so a + # reversible reaction running backwards counts by its magnitude. + magnitude = np.abs(flux) + if not (magnitude > _TRACK_EPS).any() or magnitude.std() <= 0: continue - fit = linregress(enforced, flux) + fit = linregress(enforced, magnitude) slope, corr = float(fit.slope), float(fit.rvalue) - if abs(corr) < corr_threshold or abs(slope) < flux_eps: + if abs(corr) < corr_threshold or slope == 0: continue - # Classify on the slope of |flux| vs the enforced product flux — the - # criterion the docstring states (|flux| rises = amplify, etc.), not an - # endpoint-only ``abs(final) vs abs(initial)`` comparison, which can - # mislabel a track whose first/last values straddle a peak/trough even - # though its overall trend runs the other way. Keep ``knockout`` for - # tracks the regression drives essentially to zero. - abs_fit = linregress(enforced, np.abs(flux)) - abs_slope = float(abs_fit.slope) - if abs(final) < flux_eps and abs_slope < 0: - ttype = "knockout" - elif abs_slope > 0: + if slope > 0: ttype = "amplify" + elif abs(final) <= _KNOCKOUT_EPS: + ttype = "knockout" else: ttype = "knockdown" rows.append({ diff --git a/src/raven_toolbox/analysis/reporter.py b/src/raven_toolbox/analysis/reporter.py index 4d9e9113..ea15b6df 100644 --- a/src/raven_toolbox/analysis/reporter.py +++ b/src/raven_toolbox/analysis/reporter.py @@ -47,8 +47,17 @@ def _gene_z(pvalues: dict[str, float]) -> dict[str, float]: def _reporter_one(model: cobra.Model, gene_z: dict[str, float], test: str) -> ReporterResult: z_values = np.fromiter(gene_z.values(), dtype=float) - mu = float(z_values.mean()) if z_values.size else 0.0 - sigma = float(z_values.std(ddof=0)) if z_values.size else 0.0 + # Moments are taken on z-scores centred on c, the first z-score plus the mean of the + # z-scores relative to it (two passes). A spread of a few ulps between otherwise identical + # z-scores is then not lost in the rounding of the mean, also when the first z-score is an outlier. + z_first = float(z_values[0]) if z_values.size else 0.0 + centre = z_first + float((z_values - z_first).mean()) if z_values.size else 0.0 + centred = z_values - centre + mu_c = float(centred.mean()) if z_values.size else 0.0 + sigma = float(centred.std(ddof=0)) if z_values.size else 0.0 + # The spread is tested for being 0 up to 100 * eps(max|z|), as in RAVEN. + spread = float(z_values.max() - z_values.min()) if z_values.size else 0.0 + has_spread = spread > 100 * float(np.spacing(np.abs(z_values).max())) if z_values.size else False rows = [] for met in model.metabolites: @@ -57,9 +66,9 @@ def _reporter_one(model: cobra.Model, gene_z: dict[str, float], test: str) -> Re continue zs = np.array([gene_z[g] for g in neighbours]) n = zs.size - raw = zs.sum() / math.sqrt(n) + raw_c = (zs - centre).sum() / math.sqrt(n) # Exact background correction for sampling-with-replacement (see module doc). - corrected = (raw - math.sqrt(n) * mu) / sigma if sigma > 0 else 0.0 + corrected = (raw_c - math.sqrt(n) * mu_c) / sigma if has_spread else 0.0 if not math.isfinite(corrected): # never on clamped z-scores; guards the p-value corrected = 0.0 rows.append( @@ -70,7 +79,7 @@ def _reporter_one(model: cobra.Model, gene_z: dict[str, float], test: str) -> Re "p_value": float(1.0 - norm.cdf(corrected)), "n_genes": n, "mean_z": float(zs.mean()), - "std_z": float(zs.std(ddof=1)) if n > 1 else float("nan"), + "std_z": float(zs.std(ddof=1)) if n > 1 else 0.0, } ) table = pd.DataFrame(rows, columns=["metabolite", "name", "z_score", "p_value", "n_genes", "mean_z", "std_z"]) @@ -86,11 +95,14 @@ def reporter_metabolites( ) -> list[ReporterResult]: """Compute Reporter Metabolites from per-gene differential-expression p-values. - ``gene_pvalues`` maps gene id → p-value (genes not in the model, or with a NaN or - out-of-``[0, 1]`` p-value, are dropped — a stray invalid p-value would otherwise - turn the whole result NaN). If ``gene_fold_changes`` (gene id → log fold change) - is given, two extra results are returned for the up- (fc ≥ 0) and down- (fc < 0) - regulated gene subsets, in addition to ``"all"``. + ``gene_pvalues`` maps gene id → p-value. Genes not in the model and genes with a + NaN (or ``None``) p-value are dropped. Any other p-value outside ``[0, 1]``, + including ±inf and for genes that are not in the model, raises ``ValueError`` + (RAVEN raises ``RAVEN:badInput``). If ``gene_fold_changes`` (gene id → log fold change) + has at least one non-zero value, two extra results are returned for the up- (fc ≥ 0) + and down- (fc < 0) regulated gene subsets, in addition to ``"all"``; with all fold + changes zero (or none given) only ``"all"`` is returned, as in RAVEN. A metabolite + with a single scored neighbouring gene has ``std_z`` 0. The background correction uses the population standard deviation (``ddof=0``) of the scored Z-scores, which the closed form requires. The @@ -104,16 +116,25 @@ def reporter_metabolites( z-scores smaller by a factor ``√(N/(N−1))`` for N scored genes; the ranking is the same. """ + given = {g: float(p) for g, p in gene_pvalues.items() if p is not None and not math.isnan(p)} + invalid = [(g, p) for g, p in given.items() if not 0.0 <= p <= 1.0] + if invalid: + gene, value = invalid[0] + msg = ( + f"The P-values must lie between 0 and 1, but {len(invalid)} of them do not. " + f"The first is {value:g} (gene {gene})" + ) + raise ValueError(msg) model_genes = {g.id for g in model.genes} - scored = { - g: float(p) - for g, p in gene_pvalues.items() - if g in model_genes and p is not None and not math.isnan(p) and 0.0 <= p <= 1.0 - } + scored = {g: p for g, p in given.items() if g in model_genes} gene_z = _gene_z(scored) results = [_reporter_one(model, gene_z, "all")] - if gene_fold_changes is not None: + # RAVEN's any(geneFoldChanges) runs after genes outside the model or without a + # p-value are dropped, and ignores NaN. + if gene_fold_changes is not None and any( + v != 0 and not math.isnan(v) for g, v in gene_fold_changes.items() if g in scored + ): up = {g: z for g, z in gene_z.items() if gene_fold_changes.get(g, 0.0) >= 0} down = {g: z for g, z in gene_z.items() if gene_fold_changes.get(g, 0.0) < 0} results.append(_reporter_one(model, up, "up")) diff --git a/src/raven_toolbox/biomass/gam.py b/src/raven_toolbox/biomass/gam.py index 092a16f5..7b6fc2ab 100644 --- a/src/raven_toolbox/biomass/gam.py +++ b/src/raven_toolbox/biomass/gam.py @@ -22,11 +22,13 @@ def set_gam( Scales every metabolite participating in the biomass pseudoreaction whose ``name`` is in ``cofactor_met_names`` (e.g. ATP, ADP, H2O, - H+, phosphate) to ``±value``, preserving the sign of its current + H+, phosphate) to exactly ``±value``, preserving the sign of its current coefficient. If both ``ngam_rxn`` and ``ngam_value`` are given, the NGAM - reaction's bounds are fixed at ``(ngam_value, ngam_value)``. + reaction's bounds are fixed at ``(ngam_value, ngam_value)``. An empty + ``ngam_rxn`` counts as not given, and ``ngam_value`` without an ``ngam_rxn`` is an + error. A given ``ngam_rxn`` must be in the model, also when no ``ngam_value`` is given. Parameters ---------- @@ -50,21 +52,30 @@ def set_gam( ------- cobra.Model The (mutated) model, for chaining. + + Raises + ------ + KeyError + If ``biomass_rxn``, or a given ``ngam_rxn``, is not in the model. Both reactions + are looked up before anything is changed. + ValueError + If ``ngam_value`` is given without an ``ngam_rxn``. """ + if ngam_value is not None and not ngam_rxn: + raise ValueError("ngam_value needs an ngam_rxn: the NGAM reaction to set is not given") rxn = model.reactions.get_by_id(biomass_rxn) + ngam = model.reactions.get_by_id(ngam_rxn) if ngam_rxn else None cofactor_set = set(cofactor_met_names) - deltas: dict[cobra.Metabolite, float] = {} + targets: dict[cobra.Metabolite, float] = {} for met, coef in rxn.metabolites.items(): if met.name in cofactor_set: target = (1 if coef > 0 else -1) * float(value) - # combine=True adds delta → new total = target - deltas[met] = target - coef - if deltas: - rxn.add_metabolites(deltas, combine=True) + targets[met] = target + if targets: + rxn.add_metabolites(targets, combine=False) # replaces the coefficient, so it is exactly target - if ngam_rxn is not None and ngam_value is not None: - ngam = model.reactions.get_by_id(ngam_rxn) + if ngam is not None and ngam_value is not None: ngam.bounds = (float(ngam_value), float(ngam_value)) return model diff --git a/src/raven_toolbox/biomass/scale.py b/src/raven_toolbox/biomass/scale.py index 603d86f1..98a313ae 100644 --- a/src/raven_toolbox/biomass/scale.py +++ b/src/raven_toolbox/biomass/scale.py @@ -59,10 +59,25 @@ def rescale_pseudoreaction( pseudoreaction) is left untouched. After the rescaling the coefficient of ``config.proton_met`` is recomputed so the reaction's total ionic charge sums to zero. + + Raises + ------ + ValueError + If a metabolite of the reaction other than the proton has no charge + (``charge is None``), because the charge balance is then unknown rather than + zero. Nothing is changed in that case. + KeyError + If the pseudoreaction or the proton metabolite is not in the model. """ comp = config.get(component_name) rxn = _find_pseudoreaction(model, comp) proton_met = model.metabolites.get_by_id(config.proton_met) + uncharged = [m.id for m in rxn.metabolites if m is not proton_met and m.charge is None] + if uncharged: + raise ValueError( + f"Cannot rebalance the charge of {rxn.id!r}: the following metabolites have no " + f"charge, so the charge balance is unknown rather than zero: {', '.join(uncharged)}." + ) # Step 1: scale substrate coefficients (everything but the component # product). Use a fresh dict to avoid mutating during iteration. @@ -79,7 +94,7 @@ def rescale_pseudoreaction( # express the same with two `_set_coefficient` calls. _set_coefficient(rxn, proton_met, 0.0) total_charge = sum( - (m.charge or 0) * coef for m, coef in rxn.metabolites.items() + m.charge * coef for m, coef in rxn.metabolites.items() if m is not proton_met ) _set_coefficient(rxn, proton_met, -total_charge) @@ -102,11 +117,17 @@ def scale_biomass( an in-place ``model`` mutation, not on a stale snapshot. With ``balance_out``, that component is rescaled so the components sum - to exactly ``target_total`` g/gDW (default 1), whatever the total was - before. It must be large enough to absorb the difference: a + to ``target_total`` g/gDW (default 1), whatever the total was before. + The sum equals ``target_total`` when the rescaled pseudoreaction needs + no change of its H+ coefficient; otherwise the H+ rebalance adds to its + fraction and the total exceeds ``target_total`` by that amount. The + balancing component must be large enough to absorb the difference: a ``ValueError`` is raised, before the model is changed, when its fraction is 0 or would become 0 or negative. ``target_total`` has no - effect without ``balance_out``. + effect without ``balance_out``. Both rescalings are tried on a scratch copy + of the pseudoreactions first, so an error leaves the model unchanged. The + scratch copy is needed because the H+ rebalance adds to the fractions, which + are therefore not linear in the scaling factor. """ fractions = sum_biomass(model, config) current = fractions[component_name] @@ -122,15 +143,18 @@ def scale_biomass( f"Cannot balance with {balance_out!r}: current fraction is 0 " "(pseudoreaction missing or empty)." ) - # Fractions are linear in the substrate coefficients, so the total - # after the first rescaling is known before the model is changed. - balance_new = balance_current + target_total - (fractions["total"] - current + new_value) + scratch = _scratch_copy(model, config) + rescale_pseudoreaction(scratch, config, component_name, new_value / current) + after = sum_biomass(scratch, config) + balance_current = after[balance_out] + balance_new = balance_current + target_total - after["total"] if balance_new <= 0: raise ValueError( f"Cannot balance {component_name!r} = {new_value:g} with " f"{balance_out!r}: it would have to go from {balance_current:.4g} " f"to {balance_new:.4g} g/gDW to bring the total to {target_total:g}." ) + rescale_pseudoreaction(scratch, config, balance_out, balance_new / balance_current) factor = new_value / current rescale_pseudoreaction(model, config, component_name, factor) @@ -145,6 +169,34 @@ def scale_biomass( # --- internals --------------------------------------------------------- +def _scratch_copy(model: cobra.Model, config: BiomassConfig) -> cobra.Model: + """A model holding copies of the biomass pseudoreactions and the proton metabolite.""" + scratch = cobra.Model("scratch") + mets: dict[str, cobra.Metabolite] = {} + + def copy_met(m: cobra.Metabolite) -> cobra.Metabolite: + if m.id not in mets: + mets[m.id] = cobra.Metabolite( + m.id, formula=m.formula, name=m.name, charge=m.charge, compartment=m.compartment + ) + return mets[m.id] + + for comp in config.components: + try: + rxn = _find_pseudoreaction(model, comp) + except KeyError: + continue + if rxn.id in scratch.reactions: + continue + new = cobra.Reaction(rxn.id, name=rxn.name) + scratch.add_reactions([new]) + new.add_metabolites({copy_met(m): c for m, c in rxn.metabolites.items()}) + proton = copy_met(model.metabolites.get_by_id(config.proton_met)) + if proton not in scratch.metabolites: + scratch.add_metabolites([proton]) + return scratch + + def _component_fraction(model: cobra.Model, comp: BiomassComponent) -> float: """Compute g/gDW of a single component from its pseudoreaction.""" try: diff --git a/src/raven_toolbox/comparison/compare.py b/src/raven_toolbox/comparison/compare.py index f2844d32..5638f57d 100644 --- a/src/raven_toolbox/comparison/compare.py +++ b/src/raven_toolbox/comparison/compare.py @@ -228,10 +228,10 @@ def compare_models( else list(tasks)) task_ids = [t.id for t in task_list] task_df = pd.DataFrame(False, index=task_ids, columns=model_ids, dtype=bool) - for mid, m in zip(model_ids, models_list, strict=True): + for col, (mid, m) in enumerate(zip(model_ids, models_list, strict=True)): results = check_tasks(m, task_list) - for r in results: - task_df.at[r.id, mid] = bool(r.passed) + for row, r in enumerate(results): # by position: task ids may repeat + task_df.iat[row, col] = bool(r.passed) if not r.passed and r.error: failed.setdefault(mid, []).append(f"{r.id}: {r.error}") diff --git a/src/raven_toolbox/curation/__init__.py b/src/raven_toolbox/curation/__init__.py index 3ed6d35f..4543d33e 100644 --- a/src/raven_toolbox/curation/__init__.py +++ b/src/raven_toolbox/curation/__init__.py @@ -25,14 +25,20 @@ - **genes_df**: ``genes, geneShortNames`` + MIRIAM columns. Match key is ``genes``. - **rxns_df**: ``rxnNames, grRules, lb, ub, rev, subSystems, eccodes, - rxnNotes, rxnReferences, rxnConfidenceScores`` + MIRIAM columns. - Match key is the reaction's *stoichiometry* — same metabolites and - coefficients ⇒ same reaction. + rxnNotes, rxnReferences, rxnConfidenceScores`` + MIRIAM columns, and an + optional ``rxnIdx``. Match key is the reaction's *stoichiometry*; the same + metabolites and coefficients ⇒ same reaction. A new reaction without ``lb`` + or ``ub`` is unbounded in that direction, as in RAVEN. - **rxns_coeffs_df**: ``rxnNames, metNames, comps, coefficient``. One - row per ``(reaction, metabolite)`` pair. The ``rxnNames`` column - links each coefficient back to a row in ``rxns_df``. An optional - ``index`` first column from the legacy yeast-GEM schema is silently - ignored. + row per ``(reaction, metabolite)`` pair. When both reaction tables have a + ``rxnIdx`` column, a coefficient row belongs to the reaction with the same + index and name, so reactions that share a name stay apart; otherwise the + ``rxnNames`` column links each coefficient to a row in ``rxns_df``. + +Headers are matched case-insensitively and read under the spellings RAVEN's +``curateModelFromTables`` accepts, for example ``metFormulas`` for ``formula``, +``metCharges`` for ``charge``, ``geneShortName`` for ``geneShortNames`` and +``index`` for ``rxnIdx``. Everything after the core columns in any of the four tables is interpreted as a MIRIAM annotation: the column header becomes the diff --git a/src/raven_toolbox/curation/batch.py b/src/raven_toolbox/curation/batch.py index 70439319..d12d8bf0 100644 --- a/src/raven_toolbox/curation/batch.py +++ b/src/raven_toolbox/curation/batch.py @@ -1,13 +1,19 @@ """Generic batch-curation engine driven by tabular data.""" from __future__ import annotations +import csv +import math +import re import warnings +from collections.abc import Mapping from dataclasses import dataclass, field from pathlib import Path import cobra +import numpy as np import pandas as pd +from raven_toolbox.utils.ids import next_ids from raven_toolbox.utils.parse import subsystem_to_str #: Core columns recognised in ``mets_df``. Anything else is treated as a @@ -23,12 +29,39 @@ "subSystems", "eccodes", "rxnNotes", "rxnReferences", "rxnConfidenceScores", ) #: Required columns in ``rxns_coeffs_df``. (Linkage column ``rxnNames`` -#: + one row per ``(reaction, metabolite)`` pair.) A leading ``index`` -#: column from the legacy yeast-GEM schema is silently ignored. +#: + one row per ``(reaction, metabolite)`` pair.) An optional ``rxnIdx`` column +#: (also read as ``index``) numbers the reactions; when both tables have it, a +#: coefficient row belongs to the reaction with the same index and name. DEFAULT_CORE_RXN_COEFFS_COLUMNS: tuple[str, ...] = ( "rxnNames", "metNames", "comps", "coefficient", ) +#: Header spellings read for each column, matched case-insensitively and mapped +#: to the canonical name (the same spellings ``curateModelFromTables`` accepts). +_MET_ALIASES = { + "metNames": ("metNames", "metName"), + "comps": ("comps", "comp"), + "formula": ("formula", "formulas", "metFormula", "metFormulas"), + "charge": ("charge", "charges", "metCharge", "metCharges"), + "inchi": ("inchi", "inchis"), + "metNotes": ("metNotes", "metNote"), +} +_GENE_ALIASES = { + "genes": ("genes", "gene"), + "geneShortNames": ("geneShortNames", "geneShortName"), +} +_RXN_ALIASES = { + "rxnIdx": ("rxnIdx", "index"), + "rxnNames": ("rxnNames", "rxnName"), +} +_COEFF_ALIASES = { + "rxnIdx": ("rxnIdx", "index"), + "rxnNames": ("rxnNames", "rxnName"), + "metNames": ("metNames", "metName"), + "comps": ("comps", "comp"), + "coefficient": ("coefficient", "coefficients", "coeffs", "stoichCoeffs"), +} + @dataclass class CurationResult: @@ -64,6 +97,7 @@ def batch_curate( rxns_coeffs_df: pd.DataFrame | None = None, met_id_prefix: str = "M_", rxn_id_prefix: str = "R_", + overwrite_with_empty: bool = True, ) -> CurationResult: """Add or update metabolites / reactions / genes from data tables. @@ -84,18 +118,84 @@ def batch_curate( the cobrapy / BiGG default). New entity ids are formed by finding the largest existing zero-padded suffix matching the prefix and incrementing from there. + overwrite_with_empty + How an empty cell is applied to an entity that already exists in the model, as + ``curateModelFromTables``'s ``overwriteWithEmpty``. With True (default) a column + present in the table overwrites the stored value, so an empty cell clears it: formula, + gene short name, reaction name, subsystem, EC code, notes and InChI become empty. For + charge and confidence score an empty cell becomes 0 when at least one cell of the + column has a value, and a column without any value is not written. For a MIRIAM + column an empty cell removes the namespace from the existing entity when at least + one existing entity of the table has a value in that column. With False an empty cell + leaves the stored value as it is. Independent of this option, an empty cell never + changes a gene rule or a lower or upper bound (no bound given). A new entity takes + the cells that have a value; an empty charge or confidence score is 0 when the column + has a value in some row of the table. + Must be a ``bool`` or ``numpy.bool_``, otherwise :class:`TypeError` is raised. + + Numeric cells (charge, confidence score, lower and upper bound, reversibility, ``rxnIdx`` + and coefficient) are read as RAVEN reads them: an empty cell or ``NaN`` (any case, with an + optional sign) is not given, and other text must be a plain decimal number or ``Inf``. A + comma, a hexadecimal number or any expression raises :class:`ValueError`, and nothing is + evaluated. ``true`` and ``false`` in any case are 1 and 0 in the lower bound, upper bound and + ``rev`` cells and are not numbers in the other columns. A cell is trimmed as MATLAB's + ``strtrim`` trims it, after dropping the backspaces at its start as ``textscan`` does; the + non-breaking space, NUL and the next line character stay in the text. + + A charge must be a finite whole number (``2``, ``-3``, ``2.0`` and ``2e0`` are accepted, + ``0.5``, ``2.5`` and ``Inf`` raise :class:`ValueError`; nothing is rounded). A ``rev`` cell + must be 0, 1, ``true`` or ``false`` (any other number raises, although the value is not + otherwise used). A ``rxnIdx`` cell must be a finite whole number of at least 1; each index + may be used by one row of the reaction table (the coefficient table repeats it, one row per + metabolite), otherwise :class:`ValueError` names the index and the data rows. The two + reaction tables are matched on the index and the name, and the index is compared by value + (``1``, ``1.0`` and ``1e0`` are the same). A coefficient must be a finite number other than + zero; two rows of one reaction for the same metabolite add their coefficients. A row of a + reaction table without any content is ignored (in the coefficient table the content is the + index, reaction name, metabolite name and compartment); a row with content needs a + ``rxnIdx`` when the table has that column, and a coefficient row needs a coefficient. A + metabolite row needs a compartment. A reaction of one table without a row in the other + raises :class:`ValueError`. An empty MIRIAM column removes no annotation. + + A lower bound above the upper bound (also against the stored bound of an existing + reaction, since both bounds are set together) raises :class:`ValueError`. RAVEN's + ``curateModelFromTables`` writes such bounds as given and ``solveLP`` rejects the model when + it is run ("Invalid set of bounds for reaction(s): ..."); a cobra reaction cannot hold them. + + Every check above, the required columns, the metabolites of the coefficient rows, the bounds + (several rows that update one reaction are checked one after the other), a gene listed twice + in ``genes_df`` and the gene rules run before the model is touched, so a table that fails + leaves the model unchanged. What is not checked beforehand is an error that cobra itself + raises while the tables are applied. + + A ``grRules`` cell must be a well-formed gene rule, in the grammar of RAVEN's + ``parseGrRule``: genes joined by ``and``/``or`` in any case or by ``&``, ``&&``, ``|`` and + ``||``, with balanced parentheses. A rule with an operator that lacks an operand, empty + parentheses or two genes without an operator between them raises :class:`ValueError`, as it + raises ``RAVEN:badInput`` there. The rule of a new reaction that names a gene which neither + the model nor ``genes_df`` has also raises :class:`ValueError`; RAVEN's ``addRxns`` is + called without permission to add genes and stops with an error of its own. The rule of an + existing reaction may name a gene the model does not have. A gene id that ``genes_df`` + lists more than once raises :class:`ValueError`, whether the model has the gene or not. + + A reaction table with a header and no row is empty: when both reaction tables are empty + nothing happens, and a coefficient row without a reaction row (or the reverse) raises + :class:`ValueError`. RAVEN matches the stoichiometry by ``rxnIdx`` alone, Python by the + index and the name: for the rows ``1 A`` and ``1 B`` with coefficient rows of ``A`` and ``B`` + RAVEN adds two reactions that both take all coefficients of index 1. Returns ------- A :class:`CurationResult` summarising the changes. """ + if not isinstance(overwrite_with_empty, (bool, np.bool_)): + raise TypeError( + "overwrite_with_empty must be a bool or numpy.bool_, " + f"got {type(overwrite_with_empty).__name__}" + ) + overwrite_with_empty = bool(overwrite_with_empty) result = CurationResult() - if mets_df is not None: - _apply_mets(model, mets_df, met_id_prefix, result) - if genes_df is not None: - _apply_genes(model, genes_df, result) - if (rxns_df is None) != (rxns_coeffs_df is None): raise ValueError( "rxns_df and rxns_coeffs_df must be provided together; got " @@ -103,8 +203,26 @@ def batch_curate( f"rxns_coeffs_df=" f"{'set' if rxns_coeffs_df is not None else 'None'}." ) + + # Validation phase: reads the tables and the model and changes nothing. + if mets_df is not None: + mets_df = _prepare_mets(mets_df) + if genes_df is not None: + genes_df = _prepare_genes(genes_df) + _check_genes(model, genes_df) + plan = None if rxns_df is not None: - _apply_reactions(model, rxns_df, rxns_coeffs_df, rxn_id_prefix, result) + rxns_df, rxns_coeffs_df = _prepare_reaction_tables(rxns_df, rxns_coeffs_df) + _check_reaction_tables(model, mets_df, rxns_df, rxns_coeffs_df) + plan = _plan_reactions(model, rxns_df, rxns_coeffs_df, genes_df) + + # Apply phase. + if mets_df is not None: + _apply_mets(model, mets_df, met_id_prefix, result, overwrite_with_empty) + if genes_df is not None: + _apply_genes(model, genes_df, result, overwrite_with_empty) + if plan is not None: + _apply_reactions(model, plan, rxns_df, rxn_id_prefix, result, overwrite_with_empty) return result @@ -117,15 +235,24 @@ def batch_curate_from_tsv( rxns_coeffs_tsv: str | Path | None = None, met_id_prefix: str = "M_", rxn_id_prefix: str = "R_", + overwrite_with_empty: bool = True, ) -> CurationResult: """File-path convenience wrapper for :func:`batch_curate`. - Each path is optional. TSVs are read with pandas' default - ``read_csv(sep='\\t')``; empty cells become ``NaN`` (not the - empty string), which the engine treats as "skip this field". + Each path is optional. A table file is read as RAVEN's ``curateModelFromTables`` reads it: + the first line is the header (also when it is blank, which leaves no usable column) and + its names are trimmed and lose one pair of surrounding double quotes; a column with an empty + header is ignored; when two columns match the same core name the first one is used and the + others are ignored, and a repeated name of another column gets ``.1``, ``.2`` ... appended + (RAVEN uses both columns). The rest is tab-delimited text with every cell kept as text + (``NA``, ``null`` or ``None`` stay text and ``22`` stays ``22``), a double quote around a + cell, spaces before a cell dropped, a trailing tab on a row ignored, and a row with more + non-empty cells than the header has columns refused with :class:`ValueError` (RAVEN wraps + the extra cell into a new row). A line without any text is skipped; RAVEN makes an empty row + from a line that holds only ``""`` or from a row that ends in two or more tabs. """ def _read(path): - return pd.read_csv(path, sep="\t") if path is not None else None + return None if path is None else _read_table(path) return batch_curate( model, @@ -135,9 +262,48 @@ def _read(path): rxns_coeffs_df=_read(rxns_coeffs_tsv), met_id_prefix=met_id_prefix, rxn_id_prefix=rxn_id_prefix, + overwrite_with_empty=overwrite_with_empty, ) +# --- reading table files ------------------------------------------------ + +def _read_table(path) -> pd.DataFrame: + """A tab-delimited table file as a data frame of text, read like ``readCurationTable``. + + The first line of the file is the header, also when it is blank: it is split at the tabs, + each name is trimmed and one pair of double quotes around it is removed, and a column with + an empty name is ignored. MATLAB's ``textscan`` reads the rest: it drops the spaces before + a cell (also before its quote) and reads a quoted cell up to the closing quote; the ``csv`` + module does the same, accepts a cell longer than 131072 characters (the limit is lifted + while the file is read) and keeps a NUL character. Repeated header names stay repeated in + the data frame. + """ + path = Path(path) + limit = csv.field_size_limit(2**31 - 1) + try: + with open(path, newline="", encoding="utf-8-sig") as handle: + header_line = handle.readline() + if header_line == "": + raise ValueError(f"The file {str(path)!r} is empty") + reader = csv.reader(handle, delimiter="\t", quotechar='"', skipinitialspace=True) + lines = [(reader.line_num + 1, row) for row in reader] + finally: + csv.field_size_limit(limit) + names = [re.sub(r'^"(.*)"$', r"\1", cell.strip(_BLANKS)) + for cell in header_line.replace("\ufeff", "").rstrip("\n").split("\t")] + width = len(names) + data = [] + for n, row in lines: + if len(row) == 0 or (len(row) == 1 and not _trim(row[0])): + continue + if len(row) > width and any(_trim(cell) for cell in row[width:]): + raise ValueError( + f"Line {n} of {path.name!r} has {len(row)} cells but the header has {width} columns") + data.append(list(row[:width]) + [""] * (width - len(row))) + return pd.DataFrame(data, columns=names, dtype=object) + + # --- metabolites ------------------------------------------------------ def _apply_mets( @@ -145,23 +311,28 @@ def _apply_mets( df: pd.DataFrame, id_prefix: str, result: CurationResult, + overwrite_with_empty: bool, ) -> None: miriam_cols = [c for c in df.columns if c not in DEFAULT_CORE_MET_COLUMNS] + any_charge = "charge" in df.columns and bool(df["charge"].map(_has_value).any()) name_index = _name_compartment_index(model) new_rows: list[pd.Series] = [] + updates: list[tuple[cobra.Metabolite, pd.Series]] = [] for _, row in df.iterrows(): - name = row["metNames"] - comp = row["comps"] - existing = name_index.get((name, comp)) + existing = name_index.get((_text(row["metNames"]), _text(row["comps"]))) if existing is not None: - _update_metabolite(existing, row, miriam_cols) - result.updated_metabolites.append(existing.id) + updates.append((existing, row)) else: new_rows.append(row) + active = _miriam_with_values(miriam_cols, [r for _, r in updates]) + for existing, row in updates: + _update_metabolite(existing, row, miriam_cols, active, overwrite_with_empty, any_charge) + result.updated_metabolites.append(existing.id) + if new_rows: - new_mets = _build_new_metabolites(model, new_rows, miriam_cols, id_prefix) + new_mets = _build_new_metabolites(model, new_rows, miriam_cols, id_prefix, any_charge) model.add_metabolites(new_mets) result.added_metabolites.extend(m.id for m in new_mets) @@ -175,17 +346,286 @@ def _apply_mets( ) +def _canonical_columns(df: pd.DataFrame, aliases: Mapping[str, tuple[str, ...]], lower_rest: bool = False) -> pd.DataFrame: + """Rename the columns of ``df`` that match an alias (case-insensitive) to the canonical name. + + As in RAVEN's ``matchColumns``, the first column that matches a canonical name is the core + column and every later column that matches it is ignored; a column with an empty header is + ignored. With ``lower_rest`` a column that matches a core reaction column in another case + is renamed to that spelling. The remaining columns keep their names, a repeated name gets + ``.1``, ``.2`` ... appended. + """ + lookup = {a.lower(): canonical for canonical, names in aliases.items() for a in names} + if lower_rest: + lookup.update({c.lower(): c for c in DEFAULT_CORE_RXN_COLUMNS}) + keep: list[int] = [] + names: list = [] + taken: set[str] = set() + for position, col in enumerate(df.columns): + if isinstance(col, str) and _trim(col) == "": + continue + canonical = lookup.get(str(col).strip().lower()) + if canonical is not None: + if canonical in taken: + continue + taken.add(canonical) + names.append(canonical) + else: + names.append(col) + keep.append(position) + out = df.iloc[:, keep].copy() + used = set(names) + seen: set = set() + for i, name in enumerate(names): + if name in seen: + k = 1 + while f"{name}.{k}" in used: + k += 1 + names[i] = f"{name}.{k}" + used.add(names[i]) + seen.add(name) + out.columns = names + return out + + +def _require_columns(df: pd.DataFrame, columns, label: str) -> None: + missing = [c for c in columns if c not in df.columns] + if missing: + raise ValueError(f"{label} is missing required columns: {missing}") + + +def _prepare_mets(df: pd.DataFrame) -> pd.DataFrame: + """Metabolite table with canonical column names and a numeric charge column; every cell is + checked before anything is applied.""" + df = _canonical_columns(df, _MET_ALIASES) + _require_columns(df, ("metNames", "comps"), "mets_df") + for position, comp in enumerate(df["comps"], start=1): + if not _has_value(comp): + raise ValueError(f"Row {position} of the metsInfo table has no compartment") + return _parse_numeric_columns( + df, ("charge",), "metsInfo", checks={"charge": _check_whole_charge}).astype(object) + + +def _prepare_genes(df: pd.DataFrame) -> pd.DataFrame: + df = _canonical_columns(df, _GENE_ALIASES) + _require_columns(df, ("genes",), "genes_df") + return df.astype(object) + + +def _check_genes(model: cobra.Model, genes_df: pd.DataFrame) -> None: + """Raise :class:`ValueError` for a gene id that the table lists more than once. + + This holds for a gene the model has as well as for a new one: two rows for one gene would + say two things about it (RAVEN would add a second gene with the same id). + """ + seen: set[str] = set() + repeated: list[str] = [] + for value in genes_df["genes"]: + gid = _text(value) + if gid in seen and gid not in repeated: + repeated.append(gid) + seen.add(gid) + if repeated: + raise ValueError(f"These genes are listed more than once in genes_df: {repeated}") + + +def _row_has_content(row, columns=None) -> bool: + values = row if columns is None else [row[c] for c in columns] + return any(_has_value(v) for v in values) + + +_RULE_TOKEN = re.compile(r"\(|\)|[^\s()]+") +_AND_WORDS = frozenset({"and", "&", "&&"}) +_OR_WORDS = frozenset({"or", "|", "||"}) +#: Spellings of the operators that cobra reads as they are. +_COBRA_OPERATORS = frozenset({"and", "AND", "&", "or", "OR", "|"}) + + +def _rule_tokens(rule: str) -> list[tuple[str, str]]: + """Split a gene rule into ``(kind, text)`` pairs; kind is ``(``, ``)``, ``and``, ``or`` or ``gene``. + + The operators are ``and``/``or`` in any case and the symbols ``&``, ``&&``, ``|`` and ``||``, + as RAVEN's ``parseGrRule`` reads them; every other word is a gene id. + """ + tokens = [] + for match in _RULE_TOKEN.finditer(rule): + word = match.group() + low = word.lower() + if word in ("(", ")"): + tokens.append((word, word)) + elif low in _AND_WORDS: + tokens.append(("and", word)) + elif low in _OR_WORDS: + tokens.append(("or", word)) + else: + tokens.append(("gene", word)) + return tokens + + +def _rule_genes(rule: str) -> list[str]: + """Gene ids of a gene rule; :class:`ValueError` when the rule is not well formed. + + The grammar is that of RAVEN's ``parseGrRule``: ``expr := term (or term)*``, + ``term := factor (and factor)*``, ``factor := gene | ( expr )``. + """ + tokens = _rule_tokens(rule) + genes: list[str] = [] + + def factor(pos: int) -> int: + if pos >= len(tokens): + raise ValueError("the rule ends after an operator") + kind, text = tokens[pos] + if kind == "(": + pos = expr(pos + 1) + if pos >= len(tokens) or tokens[pos][0] != ")": + raise ValueError("unbalanced parentheses") + return pos + 1 + if kind == "gene": + genes.append(text) + return pos + 1 + raise ValueError(f'expected a gene but found "{text}"') + + def term(pos: int) -> int: + pos = factor(pos) + while pos < len(tokens) and tokens[pos][0] == "and": + pos = factor(pos + 1) + return pos + + def expr(pos: int) -> int: + pos = term(pos) + while pos < len(tokens) and tokens[pos][0] == "or": + pos = term(pos + 1) + return pos + + if not tokens: + return genes + end = expr(0) + if end < len(tokens): + raise ValueError(f'unexpected "{tokens[end][1]}"') + return genes + + +def _row_rule_genes(row: pd.Series) -> list[str]: + """Gene ids of the ``grRules`` cell of a reaction row; :class:`ValueError` when it is malformed.""" + rule = row.get("grRules") + if not _has_value(rule): + return [] + try: + return _rule_genes(str(rule)) + except ValueError as exc: + raise ValueError( + f'The grRules value "{rule}" of reaction {_text(row["rxnNames"])!r} is not a valid gene rule: {exc}' + ) from None + + +def _stored_gene_rule(rule) -> str: + """The rule as cobra reads it: operator spellings cobra does not know (``And``, ``&&``, + ``||``) are written as ``and``/``or``; everything else is kept as it is.""" + text = str(rule) + pieces: list[str] = [] + last = 0 + for match in _RULE_TOKEN.finditer(text): + word = match.group() + low = word.lower() + if word in _COBRA_OPERATORS or word in ("(", ")"): + continue + if low in _AND_WORDS: + replacement = "and" + elif low in _OR_WORDS: + replacement = "or" + else: + continue + pieces.append(text[last:match.start()]) + pieces.append(replacement) + last = match.end() + pieces.append(text[last:]) + return "".join(pieces) + + +def _check_gene_rules(rxns_df: pd.DataFrame) -> None: + """Raise :class:`ValueError` for a ``grRules`` cell that is not a well-formed gene rule. + + The grammar is RAVEN's (``parseGrRule``), so both tools refuse the same rules: an operator + without two operands, unbalanced or empty parentheses, and two genes without an operator + between them. cobra would log a malformed rule and store an empty one. + """ + if "grRules" not in rxns_df.columns: + return + for _, row in rxns_df.iterrows(): + _row_rule_genes(row) + + +def _prepare_reaction_tables(rxns_df: pd.DataFrame, coeffs_df: pd.DataFrame): + """Reaction and coefficient tables with canonical column names and numeric columns. + + A row without any content is ignored (for the coefficient table only the reaction index, + reaction name, metabolite name and compartment count as content). A row with content needs + a reaction index when the table has that column, and a coefficient row needs a coefficient. + """ + rxns_df = _canonical_columns(rxns_df, _RXN_ALIASES, lower_rest=True) + _require_columns(rxns_df, ("rxnNames",), "rxns_df") + coeffs_df = _canonical_columns(coeffs_df, _COEFF_ALIASES) + _validate_rxns_coeffs(coeffs_df) + + rxns_df = rxns_df.reset_index(drop=True).astype(object) + rxns_df = rxns_df.loc[pd.Series( + [_row_has_content(r) for r in rxns_df.itertuples(index=False)], index=rxns_df.index, + dtype=bool)] + coeff_content = [c for c in ("rxnIdx", "rxnNames", "metNames", "comps") if c in coeffs_df.columns] + coeffs_df = coeffs_df.reset_index(drop=True).astype(object) + coeffs_df = coeffs_df.loc[pd.Series( + [_row_has_content(r, coeff_content) for _, r in coeffs_df.iterrows()], + index=coeffs_df.index, dtype=bool)] + + rxns_df = _parse_numeric_columns( + rxns_df, ("rxnIdx", "lb", "ub", "rev", "rxnConfidenceScores"), "rxnsInfo", + checks={"rxnIdx": _check_rxn_index, "rev": _check_boolean_rev}, + boolean_columns=("lb", "ub", "rev")) + coeffs_df = _parse_numeric_columns( + coeffs_df, ("rxnIdx", "coefficient"), "rxnsCoeffs", + checks={"rxnIdx": _check_rxn_index, "coefficient": _check_coefficient}) + for df, label in ((rxns_df, "rxnsInfo"), (coeffs_df, "rxnsCoeffs")): + if "rxnIdx" in df.columns and df["rxnIdx"].isna().any(): + row = int(df.index[df["rxnIdx"].isna()][0]) + 1 + raise ValueError( + f"Data row {row} of the {label} table has values but its rxnIdx is empty or NaN") + if "rxnIdx" in rxns_df.columns: + repeated = rxns_df["rxnIdx"].duplicated(keep=False) + if repeated.any(): + index = rxns_df["rxnIdx"][repeated].iloc[0] + rows = ", ".join(str(int(i) + 1) for i in rxns_df.index[rxns_df["rxnIdx"] == index]) + raise ValueError( + f"The rxnIdx {index:g} is used by data rows {rows} of the rxnsInfo table; " + "each reaction needs its own index") + rxns_df = rxns_df.astype(object) + coeffs_df = coeffs_df.astype(object) + _check_gene_rules(rxns_df) + if coeffs_df["coefficient"].isna().any(): + row = int(coeffs_df.index[coeffs_df["coefficient"].isna()][0]) + 1 + raise ValueError( + f"The coefficient in row {row} of the rxnsCoeffs table is empty or NaN; " + "every row with a reaction, metabolite or compartment needs a coefficient") + return rxns_df, coeffs_df + + def _name_compartment_index(model: cobra.Model) -> dict[tuple[str, str], cobra.Metabolite]: return {(m.name, m.compartment): m for m in model.metabolites} -def _update_metabolite(met: cobra.Metabolite, row: pd.Series, miriam_cols: list[str]) -> None: - if _has_value(row.get("formula")): - met.formula = str(row["formula"]) - if _has_value(row.get("charge")): - met.charge = _coerce_int(row["charge"]) - _apply_met_notes(met, row) - _apply_miriam(met, row, miriam_cols) +def _update_metabolite( + met: cobra.Metabolite, + row: pd.Series, + miriam_cols: list[str], + active_miriam: set[str], + overwrite_with_empty: bool, + any_charge: bool, +) -> None: + if _is_written(row, "formula", overwrite_with_empty): + met.formula = _text(row["formula"]) + if _is_written(row, "charge", overwrite_with_empty) and any_charge: + met.charge = _coerce_int(row["charge"]) if _has_value(row["charge"]) else 0 + _apply_met_notes(met, row, overwrite_with_empty) + _apply_miriam(met, row, miriam_cols, active_miriam, overwrite_with_empty) def _build_new_metabolites( @@ -193,19 +633,22 @@ def _build_new_metabolites( rows: list[pd.Series], miriam_cols: list[str], id_prefix: str, + any_charge: bool, ) -> list[cobra.Metabolite]: ids = _generate_ids(model.metabolites, id_prefix, len(rows)) new_mets: list[cobra.Metabolite] = [] for new_id, row in zip(ids, rows, strict=True): met = cobra.Metabolite( id=new_id, - name=str(row["metNames"]), - compartment=str(row["comps"]), + name=_text(row["metNames"]), + compartment=_text(row["comps"]), ) if _has_value(row.get("formula")): met.formula = str(row["formula"]) if _has_value(row.get("charge")): met.charge = _coerce_int(row["charge"]) + elif any_charge: + met.charge = 0 _apply_met_notes(met, row) _apply_miriam(met, row, miriam_cols) new_mets.append(met) @@ -214,20 +657,27 @@ def _build_new_metabolites( # --- genes ------------------------------------------------------------ -def _apply_genes(model: cobra.Model, df: pd.DataFrame, result: CurationResult) -> None: +def _apply_genes( + model: cobra.Model, df: pd.DataFrame, result: CurationResult, overwrite_with_empty: bool +) -> None: miriam_cols = [c for c in df.columns if c not in DEFAULT_CORE_GENE_COLUMNS] existing_genes = {g.id: g for g in model.genes} new_rows: list[pd.Series] = [] + updates: list[tuple[cobra.Gene, pd.Series]] = [] for _, row in df.iterrows(): - gid = str(row["genes"]) + gid = _text(row["genes"]) existing = existing_genes.get(gid) if existing is not None: - _update_gene(existing, row, miriam_cols) - result.updated_genes.append(gid) + updates.append((existing, row)) else: new_rows.append(row) + active = _miriam_with_values(miriam_cols, [r for _, r in updates]) + for existing, row in updates: + _update_gene(existing, row, miriam_cols, active, overwrite_with_empty) + result.updated_genes.append(existing.id) + if new_rows: # cobrapy doesn't have a direct "add a free-standing gene" API; # genes are typically added via reactions. Use cobra.Gene + an @@ -235,7 +685,7 @@ def _apply_genes(model: cobra.Model, df: pd.DataFrame, result: CurationResult) - from cobra.core.gene import Gene for row in new_rows: - g = Gene(str(row["genes"])) + g = Gene(_text(row["genes"])) if _has_value(row.get("geneShortNames")): g.name = str(row["geneShortNames"]) _apply_miriam(g, row, miriam_cols) @@ -251,50 +701,110 @@ def _apply_genes(model: cobra.Model, df: pd.DataFrame, result: CurationResult) - ) -def _update_gene(gene, row: pd.Series, miriam_cols: list[str]) -> None: - if _has_value(row.get("geneShortNames")): - gene.name = str(row["geneShortNames"]) - _apply_miriam(gene, row, miriam_cols) +def _update_gene( + gene, row: pd.Series, miriam_cols: list[str], active_miriam: set[str], overwrite_with_empty: bool +) -> None: + if _is_written(row, "geneShortNames", overwrite_with_empty): + gene.name = _text(row["geneShortNames"]) + _apply_miriam(gene, row, miriam_cols, active_miriam, overwrite_with_empty) # --- reactions -------------------------------------------------------- +@dataclass +class _ReactionPlan: + """What the reaction tables do to the model: existing reactions to update and new ones to add. + + ``new`` pairs a table row with its stoichiometry as ``{(metabolite name, compartment): coefficient}``; + the metabolites are looked up when the reaction is built, after the metabolite table is applied. + """ + + updates: list[tuple[cobra.Reaction, pd.Series]] + new: list[tuple[pd.Series, dict[tuple[str, str], float]]] + + +def _plan_reactions( + model: cobra.Model, rxns_df: pd.DataFrame, coeffs_df: pd.DataFrame, + genes_df: pd.DataFrame | None = None, +) -> _ReactionPlan: + """Decide which reactions of ``rxns_df`` exist in the model, and check them. + + A reaction whose stoichiometry matches a model reaction (compared by metabolite id and + coefficient) is an update; every other one is new, also when it uses a metabolite that + ``mets_df`` adds. Raises :class:`ValueError` for a lower bound above the upper bound (several + rows that update one reaction are checked one after the other, as they are applied) and for + the gene rule of a new reaction that names a gene which neither the model nor ``genes_df`` + has (RAVEN's ``addRxns`` does not add genes either). + """ + by_index = "rxnIdx" in rxns_df.columns and "rxnIdx" in coeffs_df.columns + coeffs_by_key = _group_coeffs(coeffs_df, by_index) + name_index = _name_compartment_index(model) + by_signature = {_stoich_signature(r): r for r in model.reactions} + default_ub = _stored_default_ub(model) + known_genes = {g.id for g in model.genes} + if genes_df is not None: + known_genes |= {_text(g) for g in genes_df["genes"]} + + updates: list[tuple[cobra.Reaction, pd.Series]] = [] + new: list[tuple[pd.Series, dict[tuple[str, str], float]]] = [] + bounds: dict[str, tuple[float, float]] = {} + for _, row in rxns_df.iterrows(): + stoich = coeffs_by_key[_reaction_key(row, by_index)] + existing = None + if all(key in name_index for key in stoich): + existing = by_signature.get( + frozenset((name_index[key].id, coef) for key, coef in stoich.items())) + if existing is not None: + stored = bounds.get(existing.id, (existing.lower_bound, existing.upper_bound)) + lb = float(row["lb"]) if _has_value(row.get("lb")) else stored[0] + ub = float(row["ub"]) if _has_value(row.get("ub")) else stored[1] + bounds[existing.id] = (lb, ub) + updates.append((existing, row)) + else: + lb = float(row["lb"]) if _has_value(row.get("lb")) else -math.inf + ub = float(row["ub"]) if _has_value(row.get("ub")) else default_ub + _check_new_reaction_genes(row, known_genes) + new.append((row, stoich)) + if lb > ub: + raise ValueError( + f"The lower bound {lb:g} of reaction {_text(row['rxnNames'])!r} is larger than its " + f"upper bound {ub:g}") + return _ReactionPlan(updates, new) + + +def _check_new_reaction_genes(row: pd.Series, known_genes: set[str]) -> None: + """Raise :class:`ValueError` when the gene rule of a new reaction names an unknown gene.""" + unknown = sorted(set(_row_rule_genes(row)) - known_genes) + if unknown: + raise ValueError( + f"The gene rule of the new reaction {_text(row['rxnNames'])!r} names genes that are " + f"not in the model or in genes_df: {unknown}") + + def _apply_reactions( model: cobra.Model, + plan: _ReactionPlan, rxns_df: pd.DataFrame, - coeffs_df: pd.DataFrame, id_prefix: str, result: CurationResult, + overwrite_with_empty: bool, ) -> None: - miriam_cols = [c for c in rxns_df.columns if c not in DEFAULT_CORE_RXN_COLUMNS] - _validate_rxns_coeffs(coeffs_df) - coeffs_by_name = _group_coeffs_by_rxn_name(model, coeffs_df) - - # Build {stoichiometry_signature: existing_reaction} for the lookup - # of "is this reaction already in the model?". - by_signature = { - _stoich_signature(r): r for r in model.reactions - } + miriam_cols = [c for c in rxns_df.columns if c not in DEFAULT_CORE_RXN_COLUMNS and c != "rxnIdx"] + any_score = ( + "rxnConfidenceScores" in rxns_df.columns + and bool(rxns_df["rxnConfidenceScores"].map(_has_value).any()) + ) - new_rows: list[tuple[pd.Series, dict[cobra.Metabolite, float]]] = [] - for _, row in rxns_df.iterrows(): - rxn_name = row["rxnNames"] - if rxn_name not in coeffs_by_name: - raise ValueError( - f"Reaction {rxn_name!r} in rxns_df has no matching " - "row(s) in rxns_coeffs_df." - ) - stoich = coeffs_by_name[rxn_name] - sig = _stoich_signature_from_dict(stoich) - existing = by_signature.get(sig) - if existing is not None: - _update_reaction(existing, row, miriam_cols) - result.updated_reactions.append(existing.id) - else: - new_rows.append((row, stoich)) + active = _miriam_with_values(miriam_cols, [r for _, r in plan.updates]) + for existing, row in plan.updates: + _update_reaction(existing, row, miriam_cols, active, overwrite_with_empty, any_score) + result.updated_reactions.append(existing.id) - if new_rows: - new_rxns = _build_new_reactions(model, new_rows, miriam_cols, id_prefix) + if plan.new: + name_index = _name_compartment_index(model) + rows = [(row, {name_index[key]: coef for key, coef in stoich.items()}) + for row, stoich in plan.new] + new_rxns = _build_new_reactions(model, rows, miriam_cols, id_prefix, any_score) model.add_reactions(new_rxns) result.added_reactions.extend(r.id for r in new_rxns) @@ -308,6 +818,48 @@ def _apply_reactions( ) +def _check_reaction_tables(model, mets_df, rxns_df, coeffs_df) -> None: + """Raise :class:`ValueError` before anything is applied when a coefficient row names a + metabolite that neither the model nor ``mets_df`` has, or when a reaction of one table has + no row in the other.""" + by_index = "rxnIdx" in rxns_df.columns and "rxnIdx" in coeffs_df.columns + known = set(_name_compartment_index(model)) + if mets_df is not None: + known |= {(_text(n), _text(c)) for n, c in zip(mets_df["metNames"], mets_df["comps"], + strict=False)} + coeff_keys = set() + for _, row in coeffs_df.iterrows(): + if (_text(row["metNames"]), _text(row["comps"])) not in known: + raise ValueError( + f"Reaction {_text(row['rxnNames'])!r} references metabolite " + f"{_text(row['metNames'])!r}[{_text(row['comps'])}] which is not in the model. " + "Add it via mets_df first, or include it in the same " + "batch_curate call." + ) + coeff_keys.add(_reaction_key(row, by_index)) + rxn_keys = {_reaction_key(row, by_index) for _, row in rxns_df.iterrows()} + for _, row in rxns_df.iterrows(): + if _reaction_key(row, by_index) not in coeff_keys: + raise ValueError( + f"Reaction {_reaction_label(row, by_index)} in rxns_df has no matching " + "row(s) in rxns_coeffs_df." + ) + for _, row in coeffs_df.iterrows(): + if _reaction_key(row, by_index) not in rxn_keys: + raise ValueError( + f"Reaction {_reaction_label(row, by_index)} in rxns_coeffs_df has no matching " + "row in rxns_df." + ) + + +def _reaction_label(row: pd.Series, by_index: bool) -> str: + """The identifiers of a reaction row for an error message: its name and, when the tables are + numbered, its index.""" + if by_index: + return f"{_text(row['rxnNames'])!r} (rxnIdx {_reaction_key(row, by_index)[0]})" + return repr(_text(row["rxnNames"])) + + def _validate_rxns_coeffs(df: pd.DataFrame) -> None: missing = set(DEFAULT_CORE_RXN_COEFFS_COLUMNS) - set(df.columns) if missing: @@ -316,34 +868,29 @@ def _validate_rxns_coeffs(df: pd.DataFrame) -> None: ) -def _group_coeffs_by_rxn_name( - model: cobra.Model, - coeffs_df: pd.DataFrame, -) -> dict[str, dict[cobra.Metabolite, float]]: - """Group coefficient rows by rxn name and resolve the metabolites. - - Each value is a ``{metabolite_object: coefficient}`` dict ready for - ``cobra.Reaction.add_metabolites``. Metabolites are looked up by - ``(name, comp)``; if a coefficient references an unknown - metabolite, that's an error (the caller should add the metabolite - via ``mets_df`` first). +def _reaction_key(row: pd.Series, by_index: bool): + """Key linking a reaction row to its coefficient rows: ``(rxnIdx, name)`` or the name. + + The index is compared by value, so ``1``, ``1.0`` and ``1e0`` are the same index. """ - name_index = _name_compartment_index(model) - grouped: dict[str, dict[cobra.Metabolite, float]] = {} + if by_index: + return (int(row["rxnIdx"]), _text(row["rxnNames"])) + return _text(row["rxnNames"]) + + +def _group_coeffs(coeffs_df: pd.DataFrame, by_index: bool) -> dict: + """Group coefficient rows by reaction. + + The key is ``(rxnIdx, name)`` when ``by_index`` and the reaction name otherwise, so + reactions that share a name stay apart when the tables number them. Each value maps + ``(metabolite name, compartment)`` to the coefficient; two rows of one reaction for the + same metabolite add their coefficients, as RAVEN's ``addRxns`` does. + """ + grouped: dict = {} for _, row in coeffs_df.iterrows(): - rxn_name = str(row["rxnNames"]) - met_name = str(row["metNames"]) - comp = str(row["comps"]) - coef = float(row["coefficient"]) - met = name_index.get((met_name, comp)) - if met is None: - raise ValueError( - f"Reaction {rxn_name!r} references metabolite " - f"{met_name!r}[{comp}] which is not in the model. " - "Add it via mets_df first, or include it in the same " - "batch_curate call." - ) - grouped.setdefault(rxn_name, {})[met] = coef + stoich = grouped.setdefault(_reaction_key(row, by_index), {}) + key = (_text(row["metNames"]), _text(row["comps"])) + stoich[key] = stoich.get(key, 0.0) + float(row["coefficient"]) return grouped @@ -351,25 +898,52 @@ def _stoich_signature(rxn: cobra.Reaction) -> frozenset: return frozenset((m.id, c) for m, c in rxn.metabolites.items()) -def _stoich_signature_from_dict(stoich: dict[cobra.Metabolite, float]) -> frozenset: - return frozenset((m.id, c) for m, c in stoich.items()) +def _update_reaction( + rxn: cobra.Reaction, + row: pd.Series, + miriam_cols: list[str], + active_miriam: set[str], + overwrite_with_empty: bool, + any_score: bool, +) -> None: + if _is_written(row, "rxnNames", overwrite_with_empty): + rxn.name = _text(row["rxnNames"]) + if _has_value(row.get("grRules")): + rxn.gene_reaction_rule = _stored_gene_rule(row["grRules"]) + if _has_value(row.get("lb")) or _has_value(row.get("ub")): + # Both bounds are set together, so that a pair that is valid as a whole is accepted + # whatever the stored bounds are. + rxn.bounds = ( + float(row["lb"]) if _has_value(row.get("lb")) else rxn.lower_bound, + float(row["ub"]) if _has_value(row.get("ub")) else rxn.upper_bound, + ) + if _is_written(row, "subSystems", overwrite_with_empty): + rxn.subsystem = subsystem_to_str(row["subSystems"]) if _has_value(row["subSystems"]) else "" + if _is_written(row, "eccodes", overwrite_with_empty): + if _has_value(row["eccodes"]): + rxn.annotation["ec-code"] = str(row["eccodes"]) + else: + rxn.annotation.pop("ec-code", None) + _apply_rxn_notes(rxn, row, overwrite_with_empty, any_score) + _apply_miriam(rxn, row, miriam_cols, active_miriam, overwrite_with_empty) + +def _stored_default_ub(model: cobra.Model) -> float: + """The ``defaultUB`` stored in the model's ``notes["metaData"]``, else in + ``model.annotation``, else +inf. -def _update_reaction(rxn: cobra.Reaction, row: pd.Series, miriam_cols: list[str]) -> None: - if _has_value(row.get("rxnNames")): - rxn.name = str(row["rxnNames"]) - if _has_value(row.get("grRules")): - rxn.gene_reaction_rule = str(row["grRules"]) - if _has_value(row.get("lb")): - rxn.lower_bound = float(row["lb"]) - if _has_value(row.get("ub")): - rxn.upper_bound = float(row["ub"]) - if _has_value(row.get("subSystems")): - rxn.subsystem = subsystem_to_str(row["subSystems"]) - if _has_value(row.get("eccodes")): - rxn.annotation["ec-code"] = str(row["eccodes"]) - _apply_rxn_notes(rxn, row) - _apply_miriam(rxn, row, miriam_cols) + A reaction without ``ub`` takes this bound, as RAVEN's ``addRxns`` takes + ``model.annotation.defaultUB``. A reaction without ``lb`` keeps ``-inf``: ``addRxns`` + sets the lower bound of reactions given by ``mets`` and ``stoichCoeffs`` to ``-inf`` + before it reads ``defaultLB``, and ``curateModelFromTables`` adds reactions that way. + """ + meta = model.notes.get("metaData") + if isinstance(meta, dict) and meta.get("defaultUB") not in (None, ""): + return float(meta["defaultUB"]) + annotation = getattr(model, "annotation", None) + if isinstance(annotation, dict) and annotation.get("defaultUB") not in (None, ""): + return float(annotation["defaultUB"]) + return math.inf def _build_new_reactions( @@ -377,24 +951,28 @@ def _build_new_reactions( rows: list[tuple[pd.Series, dict[cobra.Metabolite, float]]], miriam_cols: list[str], id_prefix: str, + any_score: bool, ) -> list[cobra.Reaction]: ids = _generate_ids(model.reactions, id_prefix, len(rows)) + default_ub = _stored_default_ub(model) new_rxns: list[cobra.Reaction] = [] for new_id, (row, stoich) in zip(ids, rows, strict=True): rxn = cobra.Reaction( id=new_id, - name=str(row["rxnNames"]), - lower_bound=float(row["lb"]) if _has_value(row.get("lb")) else -1000.0, - upper_bound=float(row["ub"]) if _has_value(row.get("ub")) else 1000.0, + name=_text(row["rxnNames"]), + lower_bound=float(row["lb"]) if _has_value(row.get("lb")) else -math.inf, + upper_bound=float(row["ub"]) if _has_value(row.get("ub")) else default_ub, ) rxn.add_metabolites(stoich) if _has_value(row.get("grRules")): - rxn.gene_reaction_rule = str(row["grRules"]) + rxn.gene_reaction_rule = _stored_gene_rule(row["grRules"]) if _has_value(row.get("subSystems")): rxn.subsystem = subsystem_to_str(row["subSystems"]) if _has_value(row.get("eccodes")): rxn.annotation["ec-code"] = str(row["eccodes"]) _apply_rxn_notes(rxn, row) + if any_score and not _has_value(row.get("rxnConfidenceScores")): + rxn.notes["confidence_score"] = 0 _apply_miriam(rxn, row, miriam_cols) new_rxns.append(rxn) return new_rxns @@ -409,62 +987,208 @@ def _build_new_reactions( _RXN_NOTE_COLUMNS = (("rxnNotes", "note"), ("rxnReferences", "references")) -def _apply_met_notes(met: cobra.Metabolite, row: pd.Series) -> None: +def _apply_met_notes(met: cobra.Metabolite, row: pd.Series, overwrite_with_empty: bool = False) -> None: for col, key in _MET_NOTE_COLUMNS: - if _has_value(row.get(col)): - met.notes[key] = str(row[col]) + if _is_written(row, col, overwrite_with_empty): + _set_note(met.notes, key, row[col]) # An InChI already kept as a MIRIAM annotation is updated as well, so the # metabolite does not carry two different structures. - if _has_value(row.get("inchi")) and "inchi" in met.annotation: - met.annotation["inchi"] = str(row["inchi"]) + if _is_written(row, "inchi", overwrite_with_empty) and "inchi" in met.annotation: + if _has_value(row["inchi"]): + met.annotation["inchi"] = str(row["inchi"]) + else: + del met.annotation["inchi"] -def _apply_rxn_notes(rxn: cobra.Reaction, row: pd.Series) -> None: +def _apply_rxn_notes( + rxn: cobra.Reaction, + row: pd.Series, + overwrite_with_empty: bool = False, + any_score: bool = True, +) -> None: for col, key in _RXN_NOTE_COLUMNS: - if _has_value(row.get(col)): - rxn.notes[key] = str(row[col]) - if _has_value(row.get("rxnConfidenceScores")): - score = float(row["rxnConfidenceScores"]) + if _is_written(row, col, overwrite_with_empty): + _set_note(rxn.notes, key, row[col]) + if _is_written(row, "rxnConfidenceScores", overwrite_with_empty) and any_score: + score = float(row["rxnConfidenceScores"]) if _has_value(row["rxnConfidenceScores"]) else 0.0 rxn.notes["confidence_score"] = int(score) if score.is_integer() else score -def _apply_miriam(entity, row: pd.Series, miriam_cols: list[str]) -> None: +def _miriam_with_values(miriam_cols: list[str], rows: list[pd.Series]) -> set[str]: + """The MIRIAM columns with a value in at least one of ``rows`` (the existing entities of + the table); an empty cell of such a column removes the namespace from its entity.""" + return {c for c in miriam_cols if any(_has_value(r.get(c)) for r in rows)} + + +def _apply_miriam( + entity, + row: pd.Series, + miriam_cols: list[str], + active_miriam: frozenset[str] | set[str] = frozenset(), + overwrite_with_empty: bool = False, +) -> None: for col in miriam_cols: value = row.get(col) if _has_value(value): - entity.annotation[col] = str(value).strip() + entity.annotation[col] = _trim(value) + elif overwrite_with_empty and col in active_miriam: + entity.annotation.pop(col, None) + + +def _is_written(row: pd.Series, col: str, overwrite_with_empty: bool) -> bool: + """Whether the cell ``col`` of ``row`` is applied to an existing entity: it has a value, + or the column is present and empty cells overwrite.""" + if col not in row.index: + return False + return overwrite_with_empty or _has_value(row[col]) + + +def _text(v) -> str: + return str(v) if _has_value(v) else "" + + +def _set_note(notes: dict, key: str, value) -> None: + if _has_value(value): + notes[key] = str(value) + else: + notes.pop(key, None) + + +_DECIMAL = re.compile(r"[+-]?(\d+\.?\d*|\.\d+)([eE][+-]?\d+)?", re.ASCII) +_NAN = re.compile(r"[+-]?nan") +_INF = re.compile(r"[+-]?inf") + +#: The characters MATLAB's ``strtrim`` removes from both ends of a cell: the code points 9 to 13 +#: and 28 to 32 (tab, line feed, vertical tab, form feed, carriage return, the four separators +#: and the space), 5760, 8192 to 8198, 8200 to 8202, 8232, 8233, 8287 and 12288. NUL (0), the +#: non-breaking space (160), the next line character (133), the figure space (8199) and the +#: narrow no-break space (8239) are not removed. +_BLANKS = "".join(chr(c) for c in ( + *range(9, 14), *range(28, 33), 5760, *range(8192, 8199), *range(8200, 8203), 8232, 8233, + 8287, 12288)) + + +def _trim(text) -> str: + """A cell as RAVEN reads it: ``strtrim`` of the text, after MATLAB's ``textscan`` has dropped + the backspace characters at its start.""" + text = str(text).strip(_BLANKS) + if text.startswith("\b"): + text = text.lstrip("\b").strip(_BLANKS) + return text + + +def _number_cell(value, field: str, label: str, allow_boolean: bool = False) -> float: + """Number in a table cell, as RAVEN reads it: NaN for an empty cell or a cell saying NaN. + + ``true`` and ``false`` in any case are 1 and 0 where ``allow_boolean`` is set (the bound and + reversibility columns) and are not numbers elsewhere. Other text must be a plain decimal + number (digits with an optional sign, decimal point and exponent), ``Inf`` or ``NaN`` with + an optional sign; a comma, a hexadecimal number and any expression raise + :class:`ValueError` and nothing is evaluated. + """ + if not _has_value(value): + return math.nan + if isinstance(value, (bool, np.bool_)): + if allow_boolean: + return float(value) + raise ValueError(f'The {field} value "{value}" in the {label} table is not a number') + if isinstance(value, (int, float, np.integer, np.floating)): + return float(value) + text = _trim(value) + low = text.lower() + if low in ("true", "false") and allow_boolean: + return 1.0 if low == "true" else 0.0 + if _NAN.fullmatch(low): + return math.nan + if _INF.fullmatch(low): + return float(low) + if _DECIMAL.fullmatch(text): + return float(text) + raise ValueError(f'The {field} value "{text}" in the {label} table is not a number') + + +def _check_whole_charge(number: float, text: str, field: str, label: str) -> None: + """A charge is a finite whole number: ``0.5``, ``2.5`` and ``Inf`` raise, nothing is rounded.""" + if not math.isnan(number) and (not math.isfinite(number) or number != math.floor(number)): + raise ValueError( + f'The {field} value "{text}" in the {label} table is not a finite whole number') + + +def _check_boolean_rev(number: float, text: str, field: str, label: str) -> None: + """A reversibility cell is 0, 1, ``true`` or ``false``; other numbers raise.""" + if not math.isnan(number) and number not in (0.0, 1.0): + raise ValueError( + f'The {field} value "{text}" in the {label} table must be 0, 1, true or false') + + +def _check_rxn_index(number: float, text: str, field: str, label: str) -> None: + """A reaction index is a finite whole number of at least 1; an empty or NaN cell is not given.""" + if not math.isnan(number) and (not math.isfinite(number) or number != math.floor(number) + or number < 1): + raise ValueError( + f'The {field} value "{text}" in the {label} table is not a whole number of at least 1') + + +def _check_coefficient(number: float, text: str, field: str, label: str) -> None: + """A coefficient is a finite number other than zero; an empty or NaN cell is refused later.""" + if not math.isnan(number) and (not math.isfinite(number) or number == 0): + raise ValueError( + f'The {field} value "{text}" in the {label} table must be a finite number other than 0') + + +def _parse_numeric_columns( + df: pd.DataFrame, columns, label: str, checks=None, boolean_columns=() +) -> pd.DataFrame: + """Copy of ``df`` whose numeric columns hold floats, NaN where the cell is empty. + + ``checks`` maps a column to a function ``(number, text, column, label)`` that raises + :class:`ValueError` for a number the column does not accept. Only the columns in + ``boolean_columns`` read ``true`` and ``false``. + """ + present = [c for c in columns if c in df.columns] + if not present: + return df + checks = checks or {} + df = df.copy() + for column in present: + numbers = [] + for v in df[column]: + number = _number_cell(v, column, label, column in boolean_columns) + if column in checks: + checks[column](number, _trim(v) if _has_value(v) else "", column, label) + numbers.append(number) + df[column] = pd.Series(numbers, index=df.index, dtype=float) + return df def _has_value(v) -> bool: """Return True if ``v`` is a non-empty, non-NaN cell value.""" - if v is None: + if v is None or v is pd.NA or v is pd.NaT: return False - if isinstance(v, float) and v != v: # NaN + if isinstance(v, (float, np.floating)) and v != v: # NaN return False - if isinstance(v, str) and v.strip() == "": + if isinstance(v, str) and _trim(v) == "": return False return True def _coerce_int(v) -> int: - """Pandas reads integer-like strings as floats. Recover the int.""" - return int(round(float(v))) + """Pandas reads integer-like strings as floats. Recover the int. + + The charge column has been checked to hold finite whole numbers; any other value raises + :class:`ValueError` here too, and nothing is rounded. + """ + value = float(v) + if not math.isfinite(value) or value != math.floor(value): + raise ValueError(f"The charge {v!r} is not a finite whole number") + return int(value) def _generate_ids(existing, prefix: str, count: int) -> list[str]: - """Generate ``count`` fresh ids by incrementing from the largest - existing zero-padded suffix matching ``prefix``. + """Generate ``count`` fresh ids after the largest existing number behind ``prefix``. - e.g. with prefix ``s_`` and existing ``s_4100``, the next ids are - ``s_4101``, ``s_4102``, … Width is preserved. + Same rule as :func:`raven_toolbox.utils.ids.generate_new_ids`: with prefix ``s_`` and + existing ``s_4100``, the next ids are ``s_4101``, ``s_4102``, … The suffix has the width of + the largest existing one, or 4 digits when no id has the prefix. """ - max_n = 0 - width = 1 - for entity in existing: - if not entity.id.startswith(prefix): - continue - suffix = entity.id[len(prefix):] - if suffix.isdigit(): - max_n = max(max_n, int(suffix)) - width = max(width, len(suffix)) - return [f"{prefix}{(max_n + i):0{width}d}" for i in range(1, count + 1)] + return next_ids([e.id for e in existing], prefix, count) diff --git a/src/raven_toolbox/gapfilling/fast_lp.py b/src/raven_toolbox/gapfilling/fast_lp.py index 58822196..292b7bd8 100644 --- a/src/raven_toolbox/gapfilling/fast_lp.py +++ b/src/raven_toolbox/gapfilling/fast_lp.py @@ -30,6 +30,7 @@ from cobra.flux_analysis import find_blocked_reactions, flux_variability_analysis from cobra.util import ProcessPool +from raven_toolbox.gapfilling.templates import unique_templates from raven_toolbox.manipulation.transfer import add_reactions_from_model @@ -264,10 +265,10 @@ def fill_gaps_fast_lp( ------- FastLPResult """ - templates = _as_models(templates) + templates = unique_templates(model, _as_models(templates)) # ---- Find blocked reactions in draft ---- - blocked = set(find_blocked_reactions(model)) + blocked = set(find_blocked_reactions(model, processes=1)) if not blocked: if verbose: print("fill_gaps_fast_lp: no blocked reactions found.") @@ -287,7 +288,7 @@ def fill_gaps_fast_lp( blocked_in_model = [r for r in blocked if r in working.reactions] if blocked_in_model: fva = flux_variability_analysis( - working, reaction_list=blocked_in_model, fraction_of_optimum=0.0 + working, reaction_list=blocked_in_model, fraction_of_optimum=0.0, processes=1 ) rescuable = [ r diff --git a/src/raven_toolbox/gapfilling/fill.py b/src/raven_toolbox/gapfilling/fill.py index bdb849d5..ef885f77 100644 --- a/src/raven_toolbox/gapfilling/fill.py +++ b/src/raven_toolbox/gapfilling/fill.py @@ -2,7 +2,7 @@ *blocked* in a draft can carry flux. For the other gap-filling flavour (add the fewest template reactions until the model's -own objective becomes feasible) use ``cobra.flux_analysis.gapfill`` — just align the +own objective becomes feasible) use ``cobra.flux_analysis.gapfill``, after aligning the template's metabolite ids to the draft first, since cobra matches by id. It solves an MILP: pick the minimum-penalty subset of template reactions such that the @@ -11,7 +11,9 @@ :func:`add_reactions_from_model`), so templates in a different identifier namespace than the model still work. Per-reaction ``scores`` (higher = prefer to include) map to RAVEN's ``rxnScores``; the MILP minimises the penalty ``-score`` (default penalty -``1.0``, i.e. minimise the number of reactions added). +``1.0``, i.e. minimise the number of reactions added). As in RAVEN, a score that is not +negative is replaced by the least negative score among the template reactions (or by 0 +when none is negative), so a reaction is never rewarded for being added. """ from __future__ import annotations @@ -20,9 +22,10 @@ from dataclasses import dataclass import cobra -from cobra.flux_analysis import find_blocked_reactions, flux_variability_analysis +from cobra.flux_analysis import flux_variability_analysis from optlang.interface import Variable # untyped (resolves to Any); used for container hints +from raven_toolbox.gapfilling.templates import unique_templates from raven_toolbox.manipulation.transfer import add_reactions_from_model @@ -61,6 +64,57 @@ def _merge_templates(model: cobra.Model, templates: list[cobra.Model]) -> tuple[ return working, template_ids +def _carrying_flux( + model: cobra.Model, reaction_ids: list[str], *, threshold: float, net_production: bool +) -> set[str]: + """Ids among ``reaction_ids`` that can carry a flux larger than ``threshold`` in either direction. + + Infinite bounds count as +-10000, and with ``net_production`` the steady state is + relaxed to ``S @ v >= 0``, as in RAVEN's ``haveFlux`` on a model whose ``b`` upper + bound is infinite. + """ + if not reaction_ids: + return set() + # The rollback of ``with model`` covers reaction bounds but not constraint bounds, so + # the steady-state bounds are restored explicitly. + steady_state_ub = {met.id: model.constraints[met.id].ub for met in model.metabolites} + try: + with model as m: + for rxn in m.reactions: + if rxn.lower_bound == -math.inf: + rxn.lower_bound = -10000.0 + if rxn.upper_bound == math.inf: + rxn.upper_bound = 10000.0 + if net_production: + for met in m.metabolites: + m.constraints[met.id].ub = None + fva = flux_variability_analysis( + m, reaction_list=reaction_ids, fraction_of_optimum=0.0, processes=1 + ) + finally: + for mid, ub in steady_state_ub.items(): + model.constraints[mid].ub = ub + return { + rid for rid in reaction_ids + if rid in fva.index + and (fva.at[rid, "maximum"] > threshold or fva.at[rid, "minimum"] < -threshold) + } + + +def _penalties( + template_ids: list[str], scores: dict[str, float] | None, penalty: float +) -> dict[str, float]: + """Penalty per template reaction: ``-score`` for negative scores. + + A score that is not negative takes the penalty of the least negative score, or 0 when + no score is negative (RAVEN's ``getMinNrFluxes``). + """ + raw = {rid: scores[rid] if scores and rid in scores else -penalty for rid in template_ids} + negatives = [v for v in raw.values() if v < 0] + fallback = -max(negatives) if negatives else 0.0 + return {rid: -v if v < 0 else fallback for rid, v in raw.items()} + + def _solve_min_templates( working: cobra.Model, template_ids: list[str], @@ -103,11 +157,9 @@ def _solve_min_templates( for met in working.metabolites: working.constraints[met.id].ub = None - def pen(rid: str) -> float: - return -scores[rid] if scores and rid in scores else penalty - + penalties = _penalties(template_ids, scores, penalty) working.objective = prob.Objective( - sum(pen(rid) * indicators[rid] for rid in template_ids), direction="min" + sum(penalties[rid] * indicators[rid] for rid in template_ids), direction="min" ) working.slim_optimize() if working.solver.status != "optimal": @@ -133,55 +185,67 @@ def connect_blocked_reactions( scores: dict[str, float] | None = None, penalty: float = 1.0, allow_net_production: bool = False, - eps: float = 1.0, + threshold: float = 1.0, + min_flux: float = 0.1, ) -> GapFillResult: """Add template reactions so blocked draft reactions can carry flux. - Finds reactions that - cannot carry flux in ``model`` under its current bounds, then adds the minimum-penalty set of template - reactions that lets the (irreversible) ones carry flux, and returns the filled - model. Like RAVEN, only irreversible blocked reactions are forced — reversible - ones can carry flux trivially in the split formulation, so forcing them is - uninformative. + Finds reactions that cannot carry a flux larger than ``threshold`` in ``model`` under its current bounds, + then adds the minimum-penalty set of template reactions that lets the (irreversible) + ones carry flux, and returns the filled model. Like RAVEN, only irreversible blocked + reactions are forced; reversible ones can carry flux trivially in the split + formulation, so forcing them is uninformative. - For the *other* gap-filling flavour — adding reactions to make the model's - objective feasible — use ``cobra.flux_analysis.gapfill`` after aligning the + For the *other* gap-filling flavour (adding reactions to make the model's + objective feasible), use ``cobra.flux_analysis.gapfill`` after aligning the template's metabolite ids to the draft. The draft is expected to have exchange reactions for its nutrients (otherwise most reactions are trivially blocked). - ``eps`` is the minimum flux a blocked reaction must be able to carry to count as - rescuable, and the flux forced through it while solving. ``allow_net_production`` + ``scores`` maps a template reaction id to its score; a negative score is the penalty + of adding that reaction, and the default penalty applies to a reaction without one. + A template reaction whose id is already taken by the draft or by an earlier template is + renamed ``_