diff --git a/docs/migrating_from_gecko_matlab.md b/docs/migrating_from_gecko_matlab.md index 552b4af..ed7048b 100644 --- a/docs/migrating_from_gecko_matlab.md +++ b/docs/migrating_from_gecko_matlab.md @@ -250,8 +250,11 @@ Each row maps a MATLAB function to its geckopy equivalent. The | `kcatsFromBrenda` | `kcats_from_brenda` | Fuzzy BRENDA matching, kept as `data/kcats/brenda.tsv`. Loads the BRENDA data and KEGG distances from the adapter unless passed in. | | `kcatsFromDLKcat` | `kcats_from_dlkcat` | DLKcat (run in Docker with `run=True`, or an earlier output), kept as `data/kcats/dlkcat.tsv`. | | `kcatsFromOKP` | `kcats_from_okp` | One OpenKineticsPredictor method, kept as `data/kcats/okp_.tsv`; submits a job when there is none, writes the file once it is done. Returns `(kcat_list, done)`. | -| `saveKcatSource`, `readKcatSource`, `loadKcatSource` | `save_kcat_source`, `read_kcat_source`, `load_kcat_source` | Kcat source files: one per source, keyed by conventional reaction and enzyme set, the same file format in both tools. `load_kcat_source` maps a file to the current model and reports rows that fit no reaction. | +| `saveKcatSource`, `loadKcatSource` | `save_kcat_source`, `load_kcat_source`, `read_kcat_source` | Kcat source files: one per source, keyed by conventional reaction and enzyme set, the same file format in both tools. `load_kcat_source` maps a file to the current model and reports rows that fit no reaction. **`[Py]`** `read_kcat_source` reads a file without a model; GECKO gives the same in `loadKcatSource`'s `report.header` and `report.records`. | | `getEcRxnKeys` | `ec_rxn_keys` | Conventional reaction and enzyme set per `ec.rxns` entry. | +| `collectKcats` | `collect_kcats` | Decides one kcat per reaction from the source files, as `[kcat_assignment]` in model_adapter.toml sets (preset `gecko_median`, the default, or `gecko_classic`, or stages), without changing the model. Writes `output/_kcatDecisions.tsv`. | +| `assignKcats` | `assign_kcats` | `collect_kcats`, then sets `ec.kcat` and `ec.source` and applies the kcat constraints; replaces `merge_kcats` → `apply_kcat_list` → `apply_custom_kcats` → `fill_kcats_from_isozymes` → `assign_standard_kcat` → `apply_kcat_constraints`. In place. | +| `updateKcatSources` | `update_kcat_sources` | Gathers the source files the configuration needs that are missing, or outdated against `[kcat_sources.]`; `dry_run` only reports. | ### 5.3 Get enzyme data (`src/geckomat/get_enzyme_data`) diff --git a/src/geckopy/__init__.py b/src/geckopy/__init__.py index cbe4447..302ae3d 100644 --- a/src/geckopy/__init__.py +++ b/src/geckopy/__init__.py @@ -74,7 +74,9 @@ # Gather-kcats from geckopy.gather_kcats import ( apply_kcat_list, + assign_kcats, assign_standard_kcat, + collect_kcats, fetch_open_kinetics_predictor, fuzzy_kcat_matching, get_standard_kcat, # deprecated; alias of assign_standard_kcat @@ -91,6 +93,7 @@ save_kcat_source, select_kcat_value, # deprecated; alias of apply_kcat_list submit_open_kinetics_predictor, + update_kcat_sources, write_dlkcat_input, ) @@ -195,7 +198,9 @@ "set_prot_pool_size", # Gather-kcats "apply_kcat_list", + "assign_kcats", "assign_standard_kcat", + "collect_kcats", "fetch_open_kinetics_predictor", "fuzzy_kcat_matching", "kcats_from_brenda", @@ -210,6 +215,7 @@ "run_dlkcat", "save_kcat_source", "submit_open_kinetics_predictor", + "update_kcat_sources", "write_dlkcat_input", # get_enzyme_data "copy_ec_to_gem", diff --git a/src/geckopy/adapter/adapterParameters.tsv b/src/geckopy/adapter/adapterParameters.tsv index c4c45f8..7ad87b1 100644 --- a/src/geckopy/adapter/adapterParameters.tsv +++ b/src/geckopy/adapter/adapterParameters.tsv @@ -32,5 +32,10 @@ okp.method okp.method "CataPro" OpenKineticsPredictor method okp.handle_long_sequences okp.handleLongSequences "truncate" How OpenKineticsPredictor handles sequences that are too long okp.include_similarity_columns okp.includeSimilarityColumns true Return similarity to the method's training data okp.canonicalize_substrates okp.canonicalizeSubstrates true Canonicalize substrate SMILES before prediction +kcat_assignment.preset Kcat assignment preset: gecko_median (used when neither preset nor stages are given) or gecko_classic +kcat_assignment.isozymes true Fill reactions without a kcat from their isozymes after the stages +kcat_assignment.standard true Give the remaining reactions the standard kcat +kcat_assignment.stage Stages instead of a preset: action, sources, criteria and BRENDA limits ([[kcat_assignment.stage]]) +kcat_sources Settings per kcat source file, e.g. [kcat_sources.brenda] aggregation matlab.make_model_anaerobic matlab.makeModelAnaerobic MATLAB function that constrains the ecModel to anaerobic growth matlab.change_protein_biomass matlab.changeProteinBiomass MATLAB function that changes the protein content of the biomass diff --git a/src/geckopy/adapter/params.py b/src/geckopy/adapter/params.py index 04bdabd..f1a86d3 100644 --- a/src/geckopy/adapter/params.py +++ b/src/geckopy/adapter/params.py @@ -100,6 +100,54 @@ class MatlabParams(BaseModel): ) +class KcatStageParams(BaseModel): + """One ``[[kcat_assignment.stage]]``: decides each reaction that is still + open by the first listed source with a value.""" + model_config = ConfigDict(extra="forbid") + + action: Literal["select"] = Field(default="select", description="What the stage does.") + sources: list[str] = Field( + min_length=1, + description=( + "Source file names (custom, dlkcat, okp_catapro, ...) and the BRENDA tiers " + "database_exact, database_top and database_bottom, which read brenda.tsv; " + "first listed wins." + ), + ) + criteria: Literal["max", "min", "median", "mean"] = Field( + default="max", description="Among several values of one source for a reaction." + ) + top_origin_limit: int = Field(default=4, description="Origin limit of database_top.") + bottom_origin_limit: int = Field(default=6, description="Origin limit of database_bottom.") + wildcard_limit: int = Field(default=3, description="Wildcard limit of database_bottom.") + + +class KcatAssignmentParams(BaseModel): + """How kcats are assigned from the kcat source files (``collect_kcats``).""" + model_config = ConfigDict(extra="forbid") + + preset: Optional[Literal["gecko_classic", "gecko_median"]] = Field( + default=None, + description=( + "Kcat assignment preset; gecko_median is used when neither a preset nor " + "stages are given." + ), + ) + isozymes: bool = Field( + default=True, description="Fill reactions without a kcat from their isozymes after the stages." + ) + standard: bool = Field(default=True, description="Give the remaining reactions the standard kcat.") + stage: list[KcatStageParams] = Field( + default_factory=list, description="Stages instead of a preset." + ) + + @model_validator(mode="after") + def _preset_or_stages(self) -> "KcatAssignmentParams": + if self.preset is not None and self.stage: + raise ValueError("[kcat_assignment] gives both a preset and stages; give one of the two.") + return self + + class SourceGroupRule(BaseModel): """One trust-tier entry in ``EvotuneParams.source_groups``. @@ -323,4 +371,9 @@ class ModelParameters(BaseModel): evotune: EvotuneParams = Field(default_factory=EvotuneParams) okp: OkpParams = Field(default_factory=OkpParams) spontaneous: SpontaneousParams = Field(default_factory=SpontaneousParams) + kcat_assignment: KcatAssignmentParams = Field(default_factory=KcatAssignmentParams) + kcat_sources: dict[str, dict[str, Any]] = Field( + default_factory=dict, + description="Settings per kcat source file, e.g. kcat_sources.brenda.aggregation.", + ) matlab: MatlabParams = Field(default_factory=MatlabParams) diff --git a/src/geckopy/gather_kcats/__init__.py b/src/geckopy/gather_kcats/__init__.py index f6eeb14..3871f3c 100644 --- a/src/geckopy/gather_kcats/__init__.py +++ b/src/geckopy/gather_kcats/__init__.py @@ -3,6 +3,7 @@ Mirrors GECKO MATLAB's `src/geckomat/gather_kcats/` directory. """ from .fuzzy_kcat_matching import fuzzy_kcat_matching +from .kcat_assignment import assign_kcats, collect_kcats, update_kcat_sources from .kcat_sources import ec_rxn_keys, load_kcat_source, read_kcat_source, save_kcat_source from .kcats_from import kcats_from_brenda, kcats_from_dlkcat, kcats_from_okp from .get_standard_kcat import assign_standard_kcat, get_standard_kcat @@ -28,7 +29,9 @@ "OKPClient", "OKPError", "apply_kcat_list", + "assign_kcats", "assign_standard_kcat", + "collect_kcats", "ec_rxn_keys", "extract_enzyme_substrate_pairs", "fetch_open_kinetics_predictor", @@ -49,5 +52,6 @@ "save_kcat_source", "select_kcat_value", # deprecated; alias of apply_kcat_list "submit_open_kinetics_predictor", + "update_kcat_sources", "write_dlkcat_input", ] diff --git a/src/geckopy/gather_kcats/kcat_assignment.py b/src/geckopy/gather_kcats/kcat_assignment.py new file mode 100644 index 0000000..d4f7e09 --- /dev/null +++ b/src/geckopy/gather_kcats/kcat_assignment.py @@ -0,0 +1,397 @@ +"""Assign kcats from the kcat source files, as model_adapter.toml sets. + +``collect_kcats`` decides one kcat per reaction from the source files in +``data/kcats`` without changing the model; ``assign_kcats`` applies the +decisions; ``update_kcat_sources`` gathers the source files the +configuration needs. GECKO (MATLAB) has the same functions as +``collectKcats``, ``assignKcats`` and ``updateKcatSources``. +""" +from __future__ import annotations + +import json +import logging +import re +import shutil +from dataclasses import dataclass, field +from pathlib import Path +from typing import TYPE_CHECKING, Any, Optional + +import numpy as np +import pandas as pd + +from ..ec_model.pipeline.apply_custom_kcats import apply_custom_kcats +from ..ec_model.pipeline.apply_kcat import apply_kcat_constraints +from ..ec_model.pipeline.fill_kcats import fill_kcats_from_isozymes +from .get_standard_kcat import assign_standard_kcat +from .kcat_sources import load_kcat_source, read_kcat_source +from .kcats_from import _brenda_release, kcats_from_brenda, kcats_from_dlkcat, kcats_from_okp +from .merge_kcats import _strip_okp_prefix, merge_kcats, normalize_source +from .select_kcat_value import apply_kcat_list + +if TYPE_CHECKING: + from ..adapter import ModelAdapter + from ..databases import UniprotDB + from ..ec_model.ec_model import EcModel + +logger = logging.getLogger(__name__) + +PRESETS = ("gecko_classic", "gecko_median") +_DATABASE_TOKENS = {"brenda", "sabio_rk"} +OKP_METHODS = ("CataPro", "CatPred", "CatRange", "DLKcat", "EITLEM", "KinForm-H", "KinForm-L", + "OmniESI", "OmniESI-O2DENet", "TurNuP", "UniKP") + + +@dataclass +class Stage: + sources: list[str] + action: str = "select" + criteria: str = "max" + top_origin_limit: int = 4 + bottom_origin_limit: int = 6 + wildcard_limit: int = 3 + + +@dataclass +class KcatAssignmentConfig: + """The kcat assignment settings of an adapter, with a preset expanded.""" + preset: str + isozymes: bool + standard: bool + stages: list[Stage] + sources: dict[str, dict[str, Any]] = field(default_factory=dict) + + +def _preset_stages(preset: str) -> tuple[list[Stage], dict[str, dict[str, Any]]]: + """The presets of kcat assignment; GECKO defines the same.""" + order = ["custom", "database_top", "dlkcat", "database_bottom"] + if preset == "gecko_classic": + return [Stage(sources=order, top_origin_limit=6)], {"brenda": {"aggregation": "max"}} + if preset == "gecko_median": + return [Stage(sources=order, top_origin_limit=4)], {"brenda": {"aggregation": "median"}} + raise ValueError(f"Unknown kcat assignment preset {preset!r}. Presets: {', '.join(PRESETS)}.") + + +def kcat_assignment_config(adapter: "ModelAdapter") -> KcatAssignmentConfig: + """Read ``[kcat_assignment]`` and ``[kcat_sources]`` of the adapter. + + Without a preset or stages the preset is ``gecko_median``. + """ + a = adapter.params.kcat_assignment + if a.stage: + stages = [Stage(sources=list(s.sources), action=s.action, criteria=s.criteria, + top_origin_limit=s.top_origin_limit, bottom_origin_limit=s.bottom_origin_limit, + wildcard_limit=s.wildcard_limit) for s in a.stage] + preset, sources = "", {} + else: + preset = a.preset or "gecko_median" + stages, sources = _preset_stages(preset) + sources = {**sources, **{k: dict(v) for k, v in adapter.params.kcat_sources.items()}} + return KcatAssignmentConfig(preset=preset, isozymes=a.isozymes, standard=a.standard, + stages=stages, sources=sources) + + +def _source_file(token: str) -> str: + """The source file a stage token reads: the BRENDA tiers read brenda.tsv.""" + return "brenda" if token.startswith("database_") else token + + +def _adapter_of(model: "EcModel", adapter: Optional["ModelAdapter"]) -> "ModelAdapter": + from ..adapter import resolve_adapter + + return resolve_adapter(model, adapter, purpose="kcat assignment needs the adapter") + + +def _decide(model: "EcModel", adapter: "ModelAdapter", uniprot_db: Optional["UniprotDB"]): + """Assign kcats to ``model`` in place; return the decisions and the report.""" + config = kcat_assignment_config(adapter) + path = Path(adapter.params.path) + kcat_dir = path / "data" / "kcats" + ec = model.ec + n = len(ec.rxns) + ec.kcat = np.zeros(n, dtype=float) + ec.source = [""] * n + stage = [""] * n + candidates = np.zeros(n, dtype=int) + spread = np.full(n, np.nan) + report: dict[str, Any] = {"preset": config.preset, "missing": [], "sources": {}} + loaded: dict[str, Optional[pd.DataFrame]] = {} + legacy_custom = False + + for s, st in enumerate(config.stages, start=1): + lists: list[pd.DataFrame] = [] + used: set[str] = set() + priority: list[str] = [] + for token in st.sources: + name = _source_file(token) + if name not in loaded: + file = kcat_dir / f"{name}.tsv" + if file.is_file(): + loaded[name], report["sources"][name] = load_kcat_source(model, file) + else: + loaded[name] = None + if name == "custom" and (path / "data" / "customKcats.tsv").is_file(): + legacy_custom = True + else: + report["missing"].append(name) + logger.warning("No kcat source file %s; it is left out. update_kcat_sources gathers it.", file) + kl = loaded[name] + if kl is None or kl.empty: + continue + if name not in used: + lists.append(kl) + used.add(name) + if token.startswith("database_"): + priority.append(token) + else: + # The labels of the file's rows, at the file's place in the order. + labels = dict.fromkeys(_strip_okp_prefix(normalize_source(x)) for x in kl["source"]) + priority.extend(x for x in labels if x not in _DATABASE_TOKENS) + if not lists or not priority: + continue + merged = merge_kcats(*lists, source_priority=priority, top_origin_limit=st.top_origin_limit, + bottom_origin_limit=st.bottom_origin_limit, wildcard_limit=st.wildcard_limit) + if merged.empty: + continue + open_before = ec.kcat == 0 + updated = apply_kcat_list(model, merged, criteria=st.criteria, overwrite=False) + index = {r: i for i, r in enumerate(ec.rxns)} + valid = merged[merged["kcat"] > 0] + groups = valid.groupby("rxn_id")["kcat"] + for rxn in updated: + i = index[rxn] + if not open_before[i]: + continue + stage[i] = str(s) + values = groups.get_group(rxn).to_numpy(dtype=float) + candidates[i] = len(values) + spread[i] = np.log10(values.max()) - np.log10(values.min()) + + if legacy_custom: + before = ec.kcat.copy() + apply_custom_kcats(model, path=path / "data" / "customKcats.tsv") + changed = ec.kcat != before + for i in np.flatnonzero(changed): + stage[i], candidates[i], spread[i] = "custom", 1, 0.0 + if config.isozymes and not ec.gecko_light: + before = ec.kcat == 0 + fill_kcats_from_isozymes(model) + for i in np.flatnonzero(before & (ec.kcat > 0)): + stage[i] = "isozymes" + if config.standard: + if uniprot_db is None: + from ..databases import load_uniprot_tsv + uniprot_db = load_uniprot_tsv(path / "data" / "uniprot.tsv") + assign_standard_kcat(model, uniprot_db) + extra = len(ec.rxns) - len(stage) + stage += [""] * extra + candidates = np.concatenate([candidates, np.zeros(extra, dtype=int)]) + spread = np.concatenate([spread, np.full(extra, np.nan)]) + for i, src in enumerate(ec.source): + if src == "standard" and ec.kcat[i] > 0: + stage[i] = "standard" + + decisions = pd.DataFrame({ + "ec_rxn": list(ec.rxns), "stage": stage, "source": list(ec.source), + "kcat": np.asarray(ec.kcat, dtype=float), "candidates": candidates, "spread_log10": spread, + }) + logger.info("Kcats decided for %d of %d reactions.", int((decisions["stage"] != "").sum()), len(decisions)) + return decisions, report + + +def _write_decisions(model: "EcModel", adapter: "ModelAdapter", decisions: pd.DataFrame) -> None: + out = Path(adapter.params.path) / "output" + out.mkdir(parents=True, exist_ok=True) + decisions.to_csv(out / f"{model.id or 'ecModel'}_kcatDecisions.tsv", sep="\t", index=False, + lineterminator="\n") + + +def collect_kcats( + model: "EcModel", + *, + adapter: Optional["ModelAdapter"] = None, + uniprot_db: Optional["UniprotDB"] = None, + write_report: bool = True, +) -> tuple[pd.DataFrame, dict[str, Any], "EcModel"]: + """Decide one kcat per reaction from the kcat source files. + + Combines the kcat source files in ``data/kcats`` as ``[kcat_assignment]`` + in model_adapter.toml sets, without changing the model given; + :func:`assign_kcats` applies the decisions. + + The order is a preset (``gecko_median``, the default, or + ``gecko_classic``) or ``[[kcat_assignment.stage]]`` tables. A select + stage decides each reaction that is still open by the first listed + source with a value (see :func:`merge_kcats`); sources are source file + names and the BRENDA tiers ``database_exact``, ``database_top`` and + ``database_bottom``, which read ``brenda.tsv``. Both presets order + custom, database_top, dlkcat and database_bottom: + + - ``gecko_median`` uses the BRENDA median, and database_top holds the + matches without EC wildcards from the organism or other organisms, + with or without the substrate (origin 1 to 4). Matches with EC + wildcards and specific-activity matches (origin 5 and 6) rank below + DLKcat. + - ``gecko_classic`` uses the BRENDA maximum and ranks the + specific-activity matches above DLKcat (origin up to 6), as the GECKO + 3 tutorial. + + Against measured kcats hidden from the search, the BRENDA median of each + origin from 1 to 4 was as close as DLKcat or closer, while the maximum + was further off than the median at every origin; specific-activity + matches could not be tested. The aggregation is + ``[kcat_sources.brenda] aggregation`` (see :func:`update_kcat_sources`). + + After the stages: custom kcats of ``data/customKcats.tsv`` when a stage + lists custom but there is no ``data/kcats/custom.tsv``; with + ``isozymes``, kcats from isozymes; with ``standard``, the standard kcat + for the rest. A source file a stage names that does not exist is + reported and left out. + + Parameters + ---------- + model + The ecModel. + adapter + Model adapter (default: ``model.adapter``). + uniprot_db + UniProt data for the standard kcat (default: ``data/uniprot.tsv``). + write_report + Write the decisions to ``output/_kcatDecisions.tsv``. + + Returns + ------- + (decisions, report, ec_model) + ``decisions`` has one row per ``model.ec`` reaction: ``ec_rxn``, + ``stage`` (the stage number, ``custom``, ``isozymes``, ``standard``, + or empty without a kcat), ``source``, ``kcat``, ``candidates`` (values + the deciding source had) and ``spread_log10``. ``report`` has the + ``preset``, the ``missing`` source files and, per source file, its + load report. ``ec_model`` is a copy of the model with the decided + kcats, before ``apply_kcat_constraints``. + """ + adapter = _adapter_of(model, adapter) + ec_model = model.copy() + decisions, report = _decide(ec_model, adapter, uniprot_db) + if write_report: + _write_decisions(model, adapter, decisions) + return decisions, report, ec_model + + +def assign_kcats( + model: "EcModel", + *, + adapter: Optional["ModelAdapter"] = None, + uniprot_db: Optional["UniprotDB"] = None, + write_report: bool = True, +) -> tuple[pd.DataFrame, dict[str, Any]]: + """Assign kcats from the kcat source files and apply them, in place. + + Decides one kcat per reaction as :func:`collect_kcats` does, writes + ``model.ec.kcat`` and ``model.ec.source``, and applies the kcat + constraints. Any kcats the model had are replaced: the source files are + the input, and nothing in the model is written back to them. + + Returns + ------- + (decisions, report) + As returned by :func:`collect_kcats`. + """ + adapter = _adapter_of(model, adapter) + decisions, report = _decide(model, adapter, uniprot_db) + if write_report: + _write_decisions(model, adapter, decisions) + apply_kcat_constraints(model) + return decisions, report + + +def _okp_method(name: str, wanted: dict[str, Any]) -> str: + if "method" in wanted: + return str(wanted["method"]) + for m in OKP_METHODS: + if "okp_" + re.sub(r"[^A-Za-z0-9]+", "_", m).lower() == name: + return m + raise ValueError(f'Set the OKP method of {name} in model_adapter.toml: [kcat_sources.{name}] method = "..."') + + +def _same(a: Any, b: Any) -> bool: + norm = (lambda v: v.lower() if isinstance(v, str) else (int(v) if isinstance(v, bool) else v)) + return norm(a) == norm(b) + + +def update_kcat_sources( + model: "EcModel", + *, + dry_run: bool = False, + adapter: Optional["ModelAdapter"] = None, +) -> pd.DataFrame: + """Gather the kcat sources that the configuration needs. + + Lists the source files that ``[kcat_assignment]`` uses and gathers those + that are missing or outdated with :func:`kcats_from_brenda`, + :func:`kcats_from_dlkcat` and :func:`kcats_from_okp`. A file is outdated + when its header records other settings than ``[kcat_sources.]`` + asks for (the preset sets these for brenda) or, for brenda, another + BRENDA release than geckopy's BRENDA files. DLKcat reads + ``data/DLKcat.tsv`` if it exists and runs DLKcat (Docker) otherwise; OKP + submits or fetches a job and reports it as pending while it runs. + ``custom.tsv`` and files it cannot gather are left as they are. Calling + it again is safe. + + Returns + ------- + pandas.DataFrame + One row per source: ``source``, ``status`` (current, gathered, + pending, outdated, missing, failed or skipped) and ``note``. + """ + adapter = _adapter_of(model, adapter) + config = kcat_assignment_config(adapter) + kcat_dir = Path(adapter.params.path) / "data" / "kcats" + names = list(dict.fromkeys(_source_file(t) for st in config.stages for t in st.sources)) + rows = [] + for name in names: + wanted = config.sources.get(name, {}) + status, note = "current", "" + file = kcat_dir / f"{name}.tsv" + if name == "custom": + status, note = "skipped", "curated by hand" + elif name not in ("brenda", "dlkcat") and not name.startswith("okp_"): + status, note = "skipped", "not a source geckopy gathers" + elif not file.is_file(): + status = "missing" + else: + _, meta = read_kcat_source(file) + recorded = json.loads(meta.get("parameters") or "{}") + for key, value in wanted.items(): + if key not in recorded or not _same(recorded[key], value): + status, note = "outdated", f"{key} differs from model_adapter.toml" + break + if status == "current" and name == "brenda": + release = _brenda_release(Path(adapter.get_brenda_db_folder())) + if release and meta.get("data") and meta["data"] != release: + status, note = "outdated", f"made from {meta['data']}; geckopy has {release}" + if status in ("missing", "outdated") and not dry_run: + status, note = _gather(model, name, wanted, adapter) + rows.append({"source": name, "status": status, "note": note}) + report = pd.DataFrame(rows, columns=["source", "status", "note"]) + logger.info("Kcat sources:\n%s", report.to_string(index=False)) + return report + + +def _gather(model: "EcModel", name: str, wanted: dict[str, Any], adapter: "ModelAdapter") -> tuple[str, str]: + try: + if name == "brenda": + kcats_from_brenda(model, aggregation=wanted.get("aggregation"), + force_wildcard_level=int(wanted.get("force_wildcard_level", 0)), adapter=adapter) + elif name == "dlkcat": + run = not (Path(adapter.params.path) / "data" / "DLKcat.tsv").is_file() + if run and shutil.which("docker") is None: + return "failed", "DLKcat needs Docker to run" + kcats_from_dlkcat(model, run=run, adapter=adapter) + else: + method = _okp_method(name, wanted) + _, done = kcats_from_okp(model, method=method, save_as=name, adapter=adapter) + if not done: + return "pending", f"OKP job for {method} still running; call again later" + except Exception as err: # reported per source, as GECKO does + return "failed", str(err) + return "gathered", "" diff --git a/src/geckopy/gather_kcats/kcat_sources.py b/src/geckopy/gather_kcats/kcat_sources.py index 8741f29..1609415 100644 --- a/src/geckopy/gather_kcats/kcat_sources.py +++ b/src/geckopy/gather_kcats/kcat_sources.py @@ -11,8 +11,7 @@ tab-separated table with the columns ``rxn``, ``enzymes``, ``substrate``, ``kcat``, ``log10_sd``, ``source``, ``detail``, ``ec_code`` and ``note``, sorted by reaction, enzymes, substrate and source. GECKO (MATLAB) reads and -writes the same files with ``saveKcatSource``, ``readKcatSource`` and -``loadKcatSource``. +writes the same files with ``saveKcatSource`` and ``loadKcatSource``. """ from __future__ import annotations diff --git a/tests/test_kcat_assignment.py b/tests/test_kcat_assignment.py new file mode 100644 index 0000000..92a9863 --- /dev/null +++ b/tests/test_kcat_assignment.py @@ -0,0 +1,179 @@ +"""Tests for collect_kcats, assign_kcats and update_kcat_sources, on ecTestGEM. + +GECKO's kcatAssignmentTests (tc0079-tc0082) test the same. +""" +import shutil +from pathlib import Path + +import cobra +import numpy as np +import pydantic +import pytest + +from geckopy import ModelAdapter, make_ec_model +from geckopy.databases import load_brenda_data, load_phyl_dist, load_uniprot_tsv +from geckopy.ec_model.pipeline import apply_custom_kcats, apply_kcat_constraints, fill_kcats_from_isozymes +from geckopy.gather_kcats import ( + apply_kcat_list, + assign_kcats, + assign_standard_kcat, + collect_kcats, + ec_rxn_keys, + fuzzy_kcat_matching, + kcats_from_brenda, + kcats_from_dlkcat, + merge_kcats, + read_dlkcat_output, + read_kcat_source, + update_kcat_sources, +) + +EXAMPLE_DIR = Path(__file__).parents[1] / "examples" / "ecTestGEM" + +DLKCAT = ( + "rxns\tgenes\tsubstrates\tsmiles\tsequence\tkcat\n" + "R2_EXP_1\tG1\tm1\tC\tMRAL\t12.5\n" + "R2_EXP_1\tG2\tm1\tC\tMNTD\t6.5\n" + "R3\tG4\tm1\tC\tMDFM\t3.25\n" +) + + +def _project(tmp_path: Path, toml_extra: str = ""): + """ecTestGEM in tmp_path, with its BRENDA files as the project's BRENDA + data and DLKcat output for the test model.""" + project = tmp_path / "ecTestGEM" + shutil.copytree(EXAMPLE_DIR, project) + brenda = project / "data" / "brenda" + brenda.mkdir() + for name in ("kcat.tsv", "sa.tsv", "mw.tsv"): + shutil.copy(project / "data" / name, brenda / name) + (project / "data" / "DLKcat.tsv").write_text(DLKCAT, encoding="utf-8") + if toml_extra: + toml = project / "model_adapter.toml" + toml.write_text(toml.read_text(encoding="utf-8") + "\n" + toml_extra, encoding="utf-8") + adapter = ModelAdapter.from_folder(project) + model = make_ec_model(cobra.io.read_sbml_model(str(adapter.params.conv_gem)), adapter) + return model, adapter, project + + +def test_assign_kcats_reproduces_chain(tmp_path): + """gecko_classic gives the same model as the steps it replaces.""" + base, adapter, project = _project(tmp_path, '[kcat_assignment]\npreset = "gecko_classic"\n') + data = project / "data" + brenda = load_brenda_data(data / "brenda") + phyl = load_phyl_dist(data / "PhylDist.mat") + uniprot = load_uniprot_tsv(data / "uniprot.tsv") + + a = base.copy() + fuzzy = fuzzy_kcat_matching(a, brenda, phyl, adapter=adapter) + dlkcat = read_dlkcat_output(a, data / "DLKcat.tsv") + merged = merge_kcats(fuzzy, dlkcat, source_priority=["database_top", "dlkcat", "database_bottom"], + top_origin_limit=6) + apply_kcat_list(a, merged) + apply_custom_kcats(a, path=data / "customKcats.tsv") + fill_kcats_from_isozymes(a) + assign_standard_kcat(a, uniprot) + apply_kcat_constraints(a) + + b = base.copy() + kcats_from_brenda(b, brenda=brenda, phyl_dist=phyl, adapter=adapter) + kcats_from_dlkcat(b, adapter=adapter) + decisions, report = assign_kcats(b, adapter=adapter, uniprot_db=uniprot) + assert list(b.ec.rxns) == list(a.ec.rxns) + np.testing.assert_array_equal(b.ec.kcat, a.ec.kcat) + assert list(b.ec.source) == list(a.ec.source) + stoich = lambda m: {r.id: {x.id: c for x, c in r.metabolites.items()} for r in m.reactions} # noqa: E731 + assert stoich(b) == stoich(a) + assert report["preset"] == "gecko_classic" + assert list(decisions["ec_rxn"]) == list(b.ec.rxns) + assert (project / "output" / f"{base.id}_kcatDecisions.tsv").is_file() + assert {"custom", "1"} <= set(decisions["stage"]) + + +def test_default_preset_is_gecko_median(tmp_path): + """Without [kcat_assignment] the preset is gecko_median: the BRENDA + median, and BRENDA matches beyond origin 4 rank below DLKcat.""" + base, adapter, project = _project(tmp_path) + data = project / "data" + uniprot = load_uniprot_tsv(data / "uniprot.tsv") + a = base.copy() + fuzzy = kcats_from_brenda(a, aggregation="median", adapter=adapter) + dlkcat = kcats_from_dlkcat(a, adapter=adapter) + merged = merge_kcats(fuzzy, dlkcat, source_priority=["database_top", "dlkcat", "database_bottom"], + top_origin_limit=4) + apply_kcat_list(a, merged) + apply_custom_kcats(a, path=data / "customKcats.tsv") + fill_kcats_from_isozymes(a) + assign_standard_kcat(a, uniprot) + apply_kcat_constraints(a) + b = base.copy() + _, report = assign_kcats(b, adapter=adapter, uniprot_db=uniprot, write_report=False) + assert report["preset"] == "gecko_median" + np.testing.assert_array_equal(b.ec.kcat, a.ec.kcat) + assert list(b.ec.source) == list(a.ec.source) + + +def test_collect_kcats_leaves_model(tmp_path): + """collect_kcats decides on a copy; a missing source is reported; + custom.tsv decides first.""" + base, adapter, project = _project(tmp_path) + kcats_from_brenda(base, adapter=adapter) + decisions, report, decided = collect_kcats(base, adapter=adapter, write_report=False) + assert not np.any(base.ec.kcat) + assert report["missing"] == ["dlkcat"] + assert np.any(decided.ec.kcat > 0) + assert not (project / "output").exists() + + rxns, enzymes = ec_rxn_keys(base) + i = base.ec.rxns.index("R3") + (project / "data" / "kcats" / "custom.tsv").write_text( + "# source: custom\n# kind: curated\nrxn\tenzymes\tkcat\tsource\n" + f"{rxns[i]}\t{enzymes[i]}\t777\tcustom\n", encoding="utf-8") + decisions, _, _ = collect_kcats(base, adapter=adapter, write_report=False) + row = decisions.iloc[i] + assert (row["kcat"], row["source"], row["stage"]) == (777, "custom", "1") + + +def test_stages_from_toml(tmp_path): + """Stages from model_adapter.toml, and their errors.""" + extra = ('[kcat_assignment]\nisozymes = false\nstandard = false\n' + '[[kcat_assignment.stage]]\nsources = ["dlkcat"]\ncriteria = "min"\n' + '[kcat_sources.brenda]\naggregation = "median"\n') + base, adapter, project = _project(tmp_path, extra) + kcats_from_dlkcat(base, adapter=adapter) + decisions, report, _ = collect_kcats(base, adapter=adapter, write_report=False) + assert report["preset"] == "" + assert set(decisions.loc[decisions["kcat"] > 0, "source"]) == {"dlkcat"} + assert decisions.loc[decisions["ec_rxn"] == "R2_EXP_1", "kcat"].item() == 6.5 + + params = ModelAdapter.from_folder(project).params + assert params.kcat_assignment.standard is False + assert params.kcat_sources["brenda"]["aggregation"] == "median" + toml = project / "model_adapter.toml" + text = toml.read_text(encoding="utf-8") + toml.write_text(text.replace("isozymes = false", 'preset = "gecko_median"'), encoding="utf-8") + with pytest.raises(pydantic.ValidationError, match="both a preset and stages"): + ModelAdapter.from_folder(project) + toml.write_text(text.replace("isozymes = false", 'preset = "gecko_fast"') + .split("[[kcat_assignment.stage]]")[0], encoding="utf-8") + with pytest.raises(pydantic.ValidationError): + ModelAdapter.from_folder(project) + + +def test_update_kcat_sources(tmp_path): + """Missing and outdated sources are gathered; current ones are not.""" + base, adapter, project = _project(tmp_path, '[kcat_assignment]\npreset = "gecko_median"\n') + report = update_kcat_sources(base, dry_run=True, adapter=adapter) + assert list(report["source"]) == ["custom", "brenda", "dlkcat"] + assert list(report["status"]) == ["skipped", "missing", "missing"] + assert not (project / "data" / "kcats").exists() + report = update_kcat_sources(base, adapter=adapter) + assert list(report["status"]) == ["skipped", "gathered", "gathered"] + _, meta = read_kcat_source(project / "data" / "kcats" / "brenda.tsv") + assert meta["parameters"] == '{"aggregation":"median","force_wildcard_level":0}' + report = update_kcat_sources(base, adapter=adapter) + assert list(report["status"]) == ["skipped", "current", "current"] + # gecko_classic asks for the BRENDA maximum: brenda.tsv is outdated. + adapter.params.kcat_assignment.preset = "gecko_classic" + report = update_kcat_sources(base, dry_run=True, adapter=adapter) + assert report["status"].iloc[1] == "outdated" diff --git a/tutorials/full_ecModel/model_adapter.toml b/tutorials/full_ecModel/model_adapter.toml index 8e852b8..63bc8e9 100644 --- a/tutorials/full_ecModel/model_adapter.toml +++ b/tutorials/full_ecModel/model_adapter.toml @@ -27,6 +27,13 @@ reviewed = true [complex] taxonomic_id = 559292 +# Kcat assignment (assign_kcats): the stages of the gecko_median preset, +# with the BRENDA aggregation that kcats_from_brenda takes from +# kcat_aggregate_brenda (the maximum). +[[kcat_assignment.stage]] +sources = ["custom", "database_top", "dlkcat", "database_bottom"] +top_origin_limit = 4 + [evotune] sigma0_log_default = 0.3 rmse_threshold = 0.2 diff --git a/tutorials/full_ecModel/protocol.py b/tutorials/full_ecModel/protocol.py index cf53ee3..fad03af 100644 --- a/tutorials/full_ecModel/protocol.py +++ b/tutorials/full_ecModel/protocol.py @@ -25,14 +25,11 @@ from geckopy import ( ModelAdapter, apply_complex_data, - apply_custom_kcats, - apply_kcat_constraints, - apply_kcat_list, - assign_standard_kcat, + assign_kcats, fill_eccodes_from_database, - fill_kcats_from_isozymes, find_met_smiles, - fuzzy_kcat_matching, + kcats_from_brenda, + kcats_from_dlkcat, load_brenda_data, load_conventional_gem, load_dlkcat_ignore_lists, @@ -40,8 +37,6 @@ load_phyl_dist, load_uniprot_tsv, make_ec_model, - merge_kcats, - read_dlkcat_output, save_ec_model, set_kcat_for_reactions, set_prot_pool_size, @@ -125,7 +120,11 @@ # %% [markdown] # **STEP 18-19** Gather kcat values from BRENDA via fuzzy matching. -# The BRENDA data ships with geckopy. When BRENDA has no value for +# Each kcat source is kept as a file in `data/kcats` (`brenda.tsv`, +# `dlkcat.tsv`, `okp_.tsv`, `custom.tsv`), which `assign_kcats` +# combines in STEP 26-31; `[kcat_assignment]` in model_adapter.toml sets +# how. `update_kcat_sources(ec_model)` gathers all sources that the +# configuration needs in one call. The BRENDA data ships with geckopy. When BRENDA has no value for # S. cerevisiae itself, the value of the phylogenetically closest # organism is used, based on the KEGG phylogenetic distances in # `data/PhylDist.mat`. That file (~5 MB) is downloaded from RAVEN on @@ -137,8 +136,8 @@ if not phyl_dist_path.is_file(): download_phyl_dist(phyl_dist_path) phyl_dist = load_phyl_dist(phyl_dist_path) -kcat_list_fuzzy = fuzzy_kcat_matching(ec_model, brenda, phyl_dist) -print(f"Fuzzy BRENDA matches: {len(kcat_list_fuzzy)} rows") +kcat_list_fuzzy = kcats_from_brenda(ec_model, brenda=brenda, phyl_dist=phyl_dist) +print(f"Fuzzy BRENDA matches: {len(kcat_list_fuzzy)} rows, kept as data/kcats/brenda.tsv") # %% [markdown] # **STEP 20-22** Gather metabolite SMILES. `find_met_smiles` queries @@ -166,66 +165,44 @@ # shipped `data/DLKcat.tsv` already has output from a previous run.) # %% [markdown] -# **STEP 25** Load DLKcat output. +# **STEP 25** Load DLKcat output from `data/DLKcat.tsv`, kept as +# `data/kcats/dlkcat.tsv`. # %% -kcat_list_dlkcat = read_dlkcat_output( - ec_model, params.path / "data" / "DLKcat.tsv", -) +kcat_list_dlkcat = kcats_from_dlkcat(ec_model) print(f"DLKcat predictions: {len(kcat_list_dlkcat)} rows") # %% [markdown] -# **STEP 26** Combine kcat from BRENDA and DLKcat. For each reaction, -# the first source in the priority list that has a kcat is used: -# BRENDA matches without EC wildcard and with origin up to 6 -# (`database_top`), then DLKcat, then any remaining BRENDA match -# (`database_bottom`). See `help(merge_kcats)` for the tier -# definitions. - -# %% -kcat_list_merged = merge_kcats( - kcat_list_fuzzy, kcat_list_dlkcat, - source_priority=["database_top", "dlkcat", "database_bottom"], -) - -# %% [markdown] -# **STEP 27** Populate `ec_model.ec.kcat` from the merged list. - -# %% -apply_kcat_list(ec_model, kcat_list_merged) - -# %% [markdown] -# **STEP 28** Apply manually-curated kcat values from -# `data/customKcats.tsv`. - -# %% -apply_custom_kcats(ec_model, path=params.path / "data" / "customKcats.tsv") - -# %% [markdown] -# **STEP 29** Propagate kcat values across isozymes (sibling -# enzymes that catalyse the same reaction). - -# %% -fill_kcats_from_isozymes(ec_model) - -# %% [markdown] -# **STEP 30** Get standard kcat. Assigns a protein cost to -# reactions without a gene association, except for exchange / -# transport / pseudoreactions and any reactions listed in -# `data/pseudoRxns.tsv`. - -# %% -assign_standard_kcat(ec_model, uniprot_db) - -# %% [markdown] -# **STEP 31** Apply kcat constraints to the stoichiometric matrix. -# `apply_kcat_constraints` translates `ec.kcat`, `ec.mw`, and -# `ec.rxn_enz_mat` into the protein pseudo-substrate stoichiometry -# in `model.S`. Re-run any time `ec.kcat`, `ec.rxn_enz_mat`, or -# `ec.mw` changes. +# **STEP 26-31** Assign kcat values. `assign_kcats` combines the kcat +# source files and applies the result to `ec_model`, in the stages that +# model_adapter.toml sets: +# +# - **STEP 26** For each reaction, the first source in the priority list +# that has a kcat is used: custom kcats, then BRENDA matches without EC +# wildcard and with origin up to 4 (`database_top`), then DLKcat, then +# any remaining BRENDA match (`database_bottom`). See +# `help(merge_kcats)` for the tier definitions. (The GECKO tutorial +# uses the `gecko_classic` preset, which takes origin up to 6.) +# - **STEP 27** The kcats are written to `ec_model.ec.kcat`. +# - **STEP 28** Manually curated kcat values from `data/customKcats.tsv`, +# or from `data/kcats/custom.tsv` when that exists. +# - **STEP 29** Reactions without a kcat take one from their isozymes. +# - **STEP 30** Get standard kcat. Assigns a protein cost to reactions +# without a gene association, except for exchange / transport / +# pseudoreactions and any reactions listed in `data/pseudoRxns.tsv`. +# - **STEP 31** Apply kcat constraints to the stoichiometric matrix: +# `ec.kcat`, `ec.mw` and `ec.rxn_enz_mat` become the protein +# pseudo-substrate stoichiometry in `model.S`. Re-run +# `apply_kcat_constraints` any time `ec.kcat`, `ec.rxn_enz_mat`, or +# `ec.mw` changes. +# +# The decision per reaction is in `output/_kcatDecisions.tsv`: +# which stage and source gave its kcat, and how many values there were. +# `collect_kcats` gives the same decisions without changing the model. # %% -apply_kcat_constraints(ec_model) +kcat_decisions, kcat_report = assign_kcats(ec_model, uniprot_db=uniprot_db) +print(kcat_decisions["stage"].value_counts().to_string()) # %% [markdown] # **STEP 32** Set the upper bound of the protein pool exchange.