diff --git a/.gitattributes b/.gitattributes index 9e5a235..5f2601a 100644 --- a/.gitattributes +++ b/.gitattributes @@ -13,3 +13,6 @@ *.txt text eol=lf *.xml text eol=lf *.md text eol=lf + +# A fixture that must keep its byte-order mark and CRLF line endings. +tests/data/kcatSourceFixture.tsv -text diff --git a/docs/migrating_from_gecko_matlab.md b/docs/migrating_from_gecko_matlab.md index 28fbaca..035ddf9 100644 --- a/docs/migrating_from_gecko_matlab.md +++ b/docs/migrating_from_gecko_matlab.md @@ -247,6 +247,14 @@ Each row maps a MATLAB function to its geckopy equivalent. The | `fetchOpenKineticsPredictor` | `fetch_open_kinetics_predictor` | **`[3→4]`** New REST-based fetch (replaced `readOpenKineticsPredictorOutput`); per-row `kcatSource` reports actual provenance (CataPro / BRENDA / Sabio-RK). | | `writeOpenKineticsPredictorInput` | *(removed)* | **`[3→4]`** Replaced by `submit_open_kinetics_predictor`. | | `readOpenKineticsPredictorOutput` | *(removed)* | **`[3→4]`** Replaced by `fetch_open_kinetics_predictor`. | +| `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. | +| `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 90d3e8a..62da514 100644 --- a/src/geckopy/__init__.py +++ b/src/geckopy/__init__.py @@ -74,17 +74,26 @@ # 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 + kcats_from_brenda, + kcats_from_dlkcat, + kcats_from_okp, + load_kcat_source, merge_dlkcat_and_fuzzy_kcats, # deprecated; alias of merge_kcats merge_kcats, read_dlkcat_output, + read_kcat_source, remove_standard_kcat, run_dlkcat, + save_kcat_source, select_kcat_value, # deprecated; alias of apply_kcat_list submit_open_kinetics_predictor, + update_kcat_sources, write_dlkcat_input, ) @@ -183,15 +192,24 @@ "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", + "kcats_from_dlkcat", + "kcats_from_okp", + "load_kcat_source", "merge_dlkcat_and_fuzzy_kcats", # deprecated; alias of merge_kcats "merge_kcats", "read_dlkcat_output", + "read_kcat_source", "remove_standard_kcat", "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/adapter.py b/src/geckopy/adapter/adapter.py index d9042d5..ac6b90d 100644 --- a/src/geckopy/adapter/adapter.py +++ b/src/geckopy/adapter/adapter.py @@ -87,9 +87,11 @@ def from_folder(cls, folder: str | Path) -> "ModelAdapter": data["path"] = folder params = ModelParameters.model_validate(data) - # Resolve conv_gem relative to the folder if it is not absolute. + # Resolve project paths relative to the folder if they are not absolute. if not params.conv_gem.is_absolute(): params.conv_gem = folder / params.conv_gem + if params.spontaneous.file is not None and not params.spontaneous.file.is_absolute(): + params.spontaneous.file = folder / params.spontaneous.file # Note: a subclass in the project's `adapter.py` is not auto-discovered; # import and instantiate your ModelAdapter subclass explicitly. @@ -101,10 +103,29 @@ def from_folder(cls, folder: str | Path) -> "ModelAdapter": def get_spontaneous_reactions(self, model: "cobra.Model") -> list[str]: """Return IDs of reactions that proceed without enzyme catalysis. - The default returns an empty list. Override in a subclass to - identify spontaneous reactions for your organism. + The default reads the ``[spontaneous]`` section of the adapter: + reactions whose name contains ``name_contains``, and those listed + in ``file``. Without that section it returns an empty list. + Override in a subclass for other rules. """ - return [] + rules = self.params.spontaneous + if rules.name_contains is None and rules.file is None: + return [] + listed: set[str] = set() + if rules.file is not None: + with open(rules.file, "r", encoding="utf-8") as f: + next(f) # header + for line in f: + parts = line.rstrip("\n").split("\t") + if not parts[0]: + continue + if len(parts) == 1 or parts[1].strip().lower() in ("1", "true"): + listed.add(parts[0].strip('"')) + return [ + r.id for r in model.reactions + if r.id in listed + or (rules.name_contains is not None and rules.name_contains in (r.name or "")) + ] def get_uniprot_compatible_genes(self, in_genes: list[str]) -> list[str]: """Map model gene IDs to IDs usable in UniProt queries. diff --git a/src/geckopy/adapter/adapterParameters.tsv b/src/geckopy/adapter/adapterParameters.tsv new file mode 100644 index 0000000..7ad87b1 --- /dev/null +++ b/src/geckopy/adapter/adapterParameters.tsv @@ -0,0 +1,41 @@ +toml_key matlab_field default description +conv_gem convGEM Path to the conventional GEM, relative to the project folder +org_name org_name Scientific name of the organism +sigma sigma 0.5 Average enzyme saturation factor +p_tot Ptot 0.5 Total protein content [g protein/gDW] +f f 0.5 Fraction of enzymes in the model [g enzyme/g protein] +gr_exp gR_exp 0.41 Growth rate the model should reach when not limited by nutrient uptake [1/h] +c_source c_source "" Exchange reaction of the main carbon source +bio_rxn bioRxn "" Biomass pseudoreaction +enzyme_comp enzyme_comp "cytoplasm" Compartment of the protein pseudometabolites +kcat_aggregate_brenda kcatAggregation "max" BRENDA aggregate read by fuzzyKcatMatching: max or median +kcat_aggregate_candidates "max" Aggregation of several candidate kcats for one reaction (geckopy) +kcat_aggregate_isozymes "mean" Aggregation of isozyme kcats when filling a missing value (geckopy) +kegg.id kegg.ID "sce" KEGG organism code +kegg.gene_id kegg.geneID "kegg" Field in the KEGG entry that matches model gene identifiers +uniprot.type uniprot.type "taxonomy" Whether uniprot.id is a taxonomy or a proteome identifier +uniprot.id uniprot.ID "559292" Taxonomy or proteome identifier +uniprot.gene_id_field uniprot.geneIDfield "gene_oln" UniProt field that matches model gene identifiers +uniprot.reviewed uniprot.reviewed false Use only reviewed (Swiss-Prot) entries +complex.taxonomic_id complex.taxonomicID Taxonomy identifier for Complex Portal +spontaneous.name_contains spontaneous.name_contains Reactions whose name contains this text are spontaneous +spontaneous.file spontaneous.file TSV of spontaneous reactions, relative to the project folder: reaction identifiers, optionally a 0/1 column +evotune.sigma0_log_default evotune.sigma0logDefault 0.5 Prior standard deviation (log) for kcats whose source is in no group +evotune.source_groups evotune.kcatSources Trust tiers: group name to the ec.source values it covers +evotune.sigma0_log_source evotune.sigma0logSource Prior standard deviation (log) per group +evotune.rmse_threshold evotune.rmseThreshold 0.2 Stop when the best RMSE reaches this +evotune.max_generations evotune.maxGenerations 150 Maximum number of CMA-ES generations +evotune.max_growth_weight evotune.maxGrowthWeight 1.0 Weight of the maximum-growth RMSE against the flux RMSE +evotune.prior_penalty_weight evotune.priorPenaltyWeight 0.03 Weight of the prior term in the search objective +evotune.tie_isozymes evotune.tieIsozymes true Give isozyme copies of one reaction a single kcat +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 b709593..0db00a3 100644 --- a/src/geckopy/adapter/params.py +++ b/src/geckopy/adapter/params.py @@ -8,7 +8,7 @@ Ported from GECKO MATLAB: src/geckomat/model_adapter/adapterTemplate.m. """ from pathlib import Path -from typing import Literal, Optional +from typing import Any, Literal, Optional from pydantic import BaseModel, ConfigDict, Field, model_validator @@ -56,6 +56,98 @@ class ComplexParams(BaseModel): ) +class SpontaneousParams(BaseModel): + """Which reactions proceed without an enzyme. + + Both rules may be given; a reaction is spontaneous if either matches. + GECKO (MATLAB) reads the same keys, so a project needs no code for this. + """ + model_config = ConfigDict(extra="forbid") + + name_contains: Optional[str] = Field( + default=None, + description="Reactions whose name contains this text are spontaneous.", + ) + file: Optional[Path] = Field( + default=None, + description=( + "TSV of spontaneous reactions, relative to the project folder: " + "a header row, reaction identifiers in the first column and, " + "optionally, a second column where 1 or true marks the " + "spontaneous ones." + ), + ) + + +class MatlabParams(BaseModel): + """Settings that only GECKO (MATLAB) reads; geckopy ignores them. + + They let one model_adapter.toml configure both tools. + """ + model_config = ConfigDict(extra="forbid") + + make_model_anaerobic: Optional[str] = Field( + default=None, + description="MATLAB function that constrains the ecModel to anaerobic growth.", + ) + change_protein_biomass: Optional[str] = Field( + default=None, + description="MATLAB function that changes the protein content of the biomass.", + ) + params: dict[str, Any] = Field( + default_factory=dict, + description="Further GECKO parameters, kept under their MATLAB names.", + ) + + +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``. @@ -277,3 +369,10 @@ class ModelParameters(BaseModel): complex: ComplexParams = Field(default_factory=ComplexParams) 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 d1e0628..3871f3c 100644 --- a/src/geckopy/gather_kcats/__init__.py +++ b/src/geckopy/gather_kcats/__init__.py @@ -3,6 +3,9 @@ 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 from .merge_dlkcat_and_fuzzy_kcats import merge_dlkcat_and_fuzzy_kcats from .merge_kcats import merge_kcats, normalize_source @@ -26,19 +29,29 @@ "OKPClient", "OKPError", "apply_kcat_list", + "assign_kcats", "assign_standard_kcat", + "collect_kcats", + "ec_rxn_keys", "extract_enzyme_substrate_pairs", "fetch_open_kinetics_predictor", "format_kcat_source", "fuzzy_kcat_matching", "get_standard_kcat", # deprecated; alias of assign_standard_kcat + "kcats_from_brenda", + "kcats_from_dlkcat", + "kcats_from_okp", + "load_kcat_source", "merge_dlkcat_and_fuzzy_kcats", # deprecated; use merge_kcats "merge_kcats", "normalize_source", + "read_kcat_source", "read_dlkcat_output", "remove_standard_kcat", "run_dlkcat", + "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/fuzzy_kcat_matching/__init__.py b/src/geckopy/gather_kcats/fuzzy_kcat_matching/__init__.py index 73ee43c..573bbc2 100644 --- a/src/geckopy/gather_kcats/fuzzy_kcat_matching/__init__.py +++ b/src/geckopy/gather_kcats/fuzzy_kcat_matching/__init__.py @@ -36,6 +36,7 @@ from ._organism import filter_by_organism, resolve_organism_index if TYPE_CHECKING: + from ...adapter import ModelAdapter from ...databases import BrendaData, PhylDist from ...ec_model.ec_model import EcModel @@ -81,6 +82,7 @@ def fuzzy_kcat_matching( ec_rxns: Optional[Iterable[str]] = None, force_wildcard_level: int = 0, aggregate: Optional[str] = None, + adapter: Optional["ModelAdapter"] = None, ) -> pd.DataFrame: """Match each ec reaction to BRENDA kcat data with progressive relaxation of EC, substrate, and organism specificity. @@ -161,6 +163,9 @@ def fuzzy_kcat_matching( across whichever rows survive the EC + organism + substrate gates. See ``docs/kcat_aggregation.md`` for the empirical rationale. + adapter + Model adapter for ``org_name`` and the aggregation default + (default: ``model.adapter``). Returns ------- @@ -182,7 +187,7 @@ def fuzzy_kcat_matching( """ from ...adapter import resolve_adapter adapter = resolve_adapter( - model, + model, adapter, purpose="fuzzy_kcat_matching reads params.org_name from the adapter", ) if force_wildcard_level < 0: 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 new file mode 100644 index 0000000..8741f29 --- /dev/null +++ b/src/geckopy/gather_kcats/kcat_sources.py @@ -0,0 +1,539 @@ +"""Kcat source files: gathered kcats kept in the project, one file per source. + +A kcat source file, ``data/kcats/.tsv`` in the project folder, holds +the kcats of one source (BRENDA, DLKcat, an OpenKineticsPredictor method, +or curated values). Rows are keyed by their conventional reaction and +enzyme set rather than by ``model.ec.rxns``, so a file stays valid when the +ecModel is rebuilt, and applies to full and light models alike. + +The file is UTF-8 text with LF line endings. It starts with ``#`` header +lines (source, generator, parameters, data, kind, created), then a +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``. +""" +from __future__ import annotations + +import json +import logging +import math +import re +from collections import defaultdict +from datetime import datetime, timezone +from importlib.metadata import PackageNotFoundError, version +from pathlib import Path +from typing import TYPE_CHECKING, Any, Optional + +import numpy as np +import pandas as pd + +if TYPE_CHECKING: + from ..adapter import ModelAdapter + from ..ec_model.ec_model import EcModel + +logger = logging.getLogger(__name__) + +COLUMNS = ["rxn", "enzymes", "substrate", "kcat", "log10_sd", "source", "detail", "ec_code", "note"] +_NAME_RE = re.compile(r"^[A-Za-z0-9_.-]+$") +_EXP_RE = re.compile(r"_EXP_\d+$") +_LIGHT_RE = re.compile(r"^\d+_") +_NUMBER_RE = re.compile(r"^[-+]?([0-9]+\.?[0-9]*|\.[0-9]+)([eE][-+]?[0-9]+)?$|^[-+]?Inf$") +_META_RE = re.compile(r"^\s*#\s*([A-Za-z][A-Za-z0-9_]{0,62})\s*:\s?(.*)$", re.ASCII) + + +def ec_rxn_keys(model: "EcModel") -> tuple[list[str], list[str]]: + """Build-independent keys of the ``model.ec`` reactions. + + Returns + ------- + (base_rxns, enzyme_sets) + Per ``model.ec.rxns`` entry: the conventional reaction (without a + trailing ``_EXP_`` in a full model, or a leading ``_`` in a + light model; a reverse direction keeps its ``_REV``), and the UniProt + accessions of its enzymes, sorted and joined by ``" + "`` (empty + without enzymes). + """ + ec = model.ec + pattern = _LIGHT_RE if ec.gecko_light else _EXP_RE + base = [pattern.sub("", r) for r in ec.rxns] + mat = ec.rxn_enz_mat.tocsr() + enzymes = [] + for i in range(len(ec.rxns)): + if mat.shape[1]: + row = mat[i] + cols = row.indices[row.data != 0] + else: + cols = [] + enzymes.append(" + ".join(sorted({ec.enzymes[j] for j in cols}))) + return base, enzymes + + +def format_number(x: Any) -> str: + """Whole numbers as integers, others in the shortest decimal form that + reads back as the same float; ``"Inf"`` for infinity, ``""`` for NaN. + Matches GECKO's ``saveKcatSource``.""" + if x is None or pd.isna(x): + return "" + x = float(x) + if math.isinf(x): + return "Inf" if x > 0 else "-Inf" + if x == round(x) and abs(x) < 1e15: + return str(int(round(x))) + for digits in range(1, 18): + text = f"{x:.{digits}g}" + if float(text) == x: + return text + return repr(x) + + +def _json_string(text: str) -> str: + out = [] + for ch in text: + if ch == '"': + out.append('\\"') + elif ch == "\\": + out.append("\\\\") + elif ch in "\b\t\n\f\r": + out.append({"\b": "\\b", "\t": "\\t", "\n": "\\n", "\f": "\\f", "\r": "\\r"}[ch]) + elif ord(ch) < 32: + out.append(f"\\u{ord(ch):04x}") + else: + out.append(ch) + return '"' + "".join(out) + '"' + + +def _parameters_json(value: Any) -> str: + """The parameters line as compact JSON with keys sorted at every level. + + Numbers are written with :func:`format_number`, NaN and infinity as + ``null``; text as it is, non-ASCII included. GECKO's ``saveKcatSource`` + writes the same text. + """ + if value is None: + return "{}" + if isinstance(value, dict): + return "{" + ",".join(f"{_json_string(str(k))}:{_parameters_json_value(value[k])}" + for k in sorted(value, key=str)) + "}" + return _parameters_json_value(value) + + +def _parameters_json_value(value: Any) -> str: + if isinstance(value, dict): + return _parameters_json(value) + if isinstance(value, np.ndarray): + value = value.tolist() + if isinstance(value, (list, tuple)): + return "[" + ",".join(_parameters_json_value(v) for v in value) + "]" + if value is None: + return "null" + if isinstance(value, (bool, np.bool_)): + return "true" if value else "false" + if isinstance(value, (int, float, np.integer, np.floating)): + return "null" if not math.isfinite(float(value)) else format_number(value) + if isinstance(value, str): + return _json_string(value) + raise TypeError(f"Parameters of type {type(value).__name__} cannot be written.") + + +def _project_path(model: "EcModel", adapter: Optional["ModelAdapter"]) -> Path: + from ..adapter import resolve_adapter + + adapter = resolve_adapter(model, adapter, purpose="kcat source files need params.path") + return Path(adapter.params.path) + + +def _header_line(text: str) -> Optional[tuple[str, str]]: + """Key (in lower case) and value of a ``# key: value`` header line.""" + m = _META_RE.match(_unquote(text)) + return (m.group(1).lower(), m.group(2).strip()) if m else None + + +def _read_header_lines(path: Path) -> dict[str, str]: + """Only the ``#`` lines at the top.""" + meta: dict[str, str] = {} + with open(path, encoding="utf-8-sig", errors="replace", newline="") as f: + for line in f: + text = line.rstrip("\r\n") + if not text.strip(): + continue + if not _unquote(text).startswith("#"): + break + kv = _header_line(text) + if kv: + meta[kv[0]] = kv[1] + return meta + + +def kcat_source_file( + name: str, + *, + model: Optional["EcModel"] = None, + adapter: Optional["ModelAdapter"] = None, + folder: Optional[str | Path] = None, +) -> Path: + """Path of a kcat source file that may be written. + + Checks the name, resolves the folder (``data/kcats`` in the project + folder by default), and refuses a file curated by hand: named + ``custom.tsv``, or whose header says ``kind: curated`` or + ``source: custom`` (in any case). + """ + if not _NAME_RE.match(name): + raise ValueError(f"The source name {name!r} may only contain letters, digits, '_', '-' and '.'.") + folder = Path(folder) if folder is not None else _project_path(model, adapter) / "data" / "kcats" + path = folder / f"{name}.tsv" + if name.lower() == "custom": + raise PermissionError("custom.tsv holds curated kcats and is not written by geckopy; edit it by hand.") + if path.is_file(): + meta = _read_header_lines(path) + if meta.get("kind", "").lower() == "curated" or meta.get("source", "").lower() == "custom": + raise PermissionError(f"{path} is curated by hand and is not overwritten.") + return path + + +def _geckopy_version() -> str: + try: + return version("geckopy") + except PackageNotFoundError: + return "" + + +def _clean(value) -> str: + if isinstance(value, (list, tuple, np.ndarray)): + value = "; ".join(str(v) for v in value) + elif value is None or pd.isna(value): + return "" + return re.sub(r"[\t\r\n]", " ", str(value)) + + +def save_kcat_source( + model: "EcModel", + kcat_list: pd.DataFrame, + name: str, + *, + adapter: Optional["ModelAdapter"] = None, + folder: Optional[str | Path] = None, + generator: str = "save_kcat_source", + parameters: Optional[dict] = None, + data: str = "", + kind: str = "gathered", +) -> Path: + """Write a kcat list to the kcat source file ``.tsv``. + + Left out are rows with kcat 0, NaN or below (no value found), and rows + of ``model.ec`` reactions without enzymes, which have no enzyme + constraint and could not be told apart from their isozymes. + + Parameters + ---------- + model + The ecModel the kcat list refers to. + kcat_list + A kcat list as returned by the ``kcats_from_*`` functions, + ``fuzzy_kcat_matching``, ``read_dlkcat_output`` or ``merge_kcats``: + columns ``rxn_id`` and ``kcat``, optionally ``substrates``, + ``eccode``, ``wildcard_level``, ``origin``, ``source`` and + ``log10_sd``. + name + Name of the source; the file is ``.tsv``. + adapter + Model adapter (default: ``model.adapter``); gives the project folder. + folder + Folder to write to (default: ``data/kcats`` in the project folder). + generator, parameters, data + What produced the values, its settings (written as JSON with sorted + keys) and the version of the data they come from. + kind + ``"gathered"`` (default) or ``"curated"``. + + A file named ``custom.tsv``, or whose header says ``kind: curated`` or + ``source: custom``, is never written. + + Returns + ------- + pathlib.Path + The file written. + """ + if kind not in ("gathered", "curated"): + raise ValueError("kind must be 'gathered' or 'curated'") + path = kcat_source_file(name, model=model, adapter=adapter, folder=folder) + + base, enzymes = ec_rxn_keys(model) + index = {r: i for i, r in enumerate(model.ec.rxns)} + unknown = sorted({r for r in kcat_list["rxn_id"] if r not in index}) + if unknown: + raise ValueError(f"These kcat list reactions are not in model.ec.rxns: {unknown}") + + def col(column): + if column in kcat_list.columns: + return kcat_list[column] + default = name if column == "source" else None + return pd.Series([default] * len(kcat_list), index=kcat_list.index) + + rows = [] + no_enzymes = 0 + for rxn, kcat, subs, src, wc, origin, ec_code, sd in zip( + kcat_list["rxn_id"], kcat_list["kcat"], col("substrates"), col("source"), + col("wildcard_level"), col("origin"), col("eccode"), col("log10_sd"), + ): + if kcat is None or pd.isna(kcat) or not (math.isfinite(float(kcat)) and float(kcat) > 0): + continue + i = index[rxn] + if not enzymes[i]: + no_enzymes += 1 + continue + detail = [] + if wc is not None and pd.notna(wc): + detail.append(f"wc={int(wc)}") + if origin is not None and pd.notna(origin): + detail.append(f"origin={int(origin)}") + rows.append([ + _clean(base[i]), _clean(enzymes[i]), _clean(subs), format_number(kcat), + format_number(sd), _clean(src), "; ".join(detail), _clean(ec_code), "", + ]) + if no_enzymes: + logger.info("%d kcats of reactions without enzymes are not written to %s.", no_enzymes, path) + rows.sort(key=lambda r: (r[0], r[1], r[2], r[5], r[3], r[4], r[6], r[7], r[8])) + + path.parent.mkdir(parents=True, exist_ok=True) + created = datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ") + header = [ + f"# source: {name}", + f"# generator: {_clean(generator)} (geckopy {_geckopy_version()})", + f"# parameters: {_parameters_json(parameters)}", + f"# data: {_clean(data)}", + f"# kind: {kind}", + f"# created: {created}", + "\t".join(COLUMNS), + ] + path.write_text("\n".join(header + ["\t".join(r) for r in rows]) + "\n", + encoding="utf-8", newline="\n") + return path + + +def _unquote(text: str) -> str: + """A field without surrounding spaces and double quotes; a field wrapped + in double quotes, as a spreadsheet may save it, is unwrapped with ``""`` + for a quote inside, and spaces inside the quotes are trimmed too.""" + text = text.strip() + if len(text) >= 2 and text[0] == '"' and text[-1] == '"': + text = text[1:-1].replace('""', '"').strip() + return text + + +def _parse_number(text: str) -> float: + return float(text) if _NUMBER_RE.match(text) else math.nan + + +def read_kcat_source(path: str | Path) -> tuple[pd.DataFrame, dict[str, str]]: + """Read a kcat source file as it is, without a model. + + The file is read as UTF-8 (a leading byte-order mark is ignored). The + ``#`` lines before the column header are the header; ``#`` lines and + blank lines further down are skipped. Fields wrapped in double quotes, + as a spreadsheet may save them, are unwrapped. ``kcat`` and + ``log10_sd`` must be plain numbers (``12``, ``0.5``, ``1e-05``); other + text, such as ``5,5`` with a decimal comma, is read as NaN, which + :func:`load_kcat_source` reports. A row with more filled fields than + there are columns is an error. + + Returns + ------- + (records, meta) + ``records`` has one row per data row: ``kcat`` and ``log10_sd`` as + floats (NaN when empty or not a number), every other column as text; + columns other than the standard ones are kept, and ``line`` is the + line number of each row in the file. ``meta`` holds the ``#`` header + lines as ``{key: value}``. + """ + path = Path(path) + with open(path, encoding="utf-8-sig", errors="replace", newline="") as f: + lines = re.split(r"\r?\n", f.read()) + meta: dict[str, str] = {} + k = 0 + while k < len(lines) and (not lines[k].strip() or _unquote(lines[k]).startswith("#")): + kv = _header_line(lines[k]) + if kv: + meta[kv[0]] = kv[1] + k += 1 + if k >= len(lines): + raise ValueError(f"{path} has no column header.") + columns = [_unquote(c) for c in lines[k].split("\t")] + if not {"kcat", "source"} <= set(columns) or not {"rxn", "enzymes"} & set(columns): + raise ValueError(f"{path} needs the columns kcat and source, and rxn and/or enzymes.") + named = [c for c in columns if c] + if len(set(named)) < len(named) or "line" in columns: + raise ValueError(f"{path} has a column name twice, or a column named line.") + rows, line_numbers = [], [] + for j in range(k + 1, len(lines)): + line = lines[j] + if not line.strip() or _unquote(line).startswith("#"): + continue + fields = line.split("\t") + if len(fields) > len(columns): + if any(f.strip() for f in fields[len(columns):]): + raise ValueError(f"{path}, line {j + 1} has more fields than the {len(columns)} columns.") + fields = fields[: len(columns)] + fields += [""] * (len(columns) - len(fields)) + if not any(f.strip() for f in fields): + continue # a row of tabs only + rows.append([_unquote(f) for f in fields]) + line_numbers.append(j + 1) + keep = [i for i, c in enumerate(columns) if c] + records = pd.DataFrame([[r[i] for i in keep] for r in rows], columns=[columns[i] for i in keep], dtype=object) + for c in COLUMNS: + if c not in records.columns: + records[c] = pd.Series([""] * len(records), dtype=object) + for c in ("kcat", "log10_sd"): + records[c] = pd.Series([_parse_number(v) for v in records[c]], dtype=float) + records["line"] = pd.Series(line_numbers, dtype=int) + return records, meta + + +def _normalize_enzymes(text: str) -> str: + parts = sorted({p.strip() for p in text.split("+") if p.strip()}) + return " + ".join(parts) + + +def load_kcat_source( + model: "EcModel", + name: str | Path, + *, + adapter: Optional["ModelAdapter"] = None, +) -> tuple[pd.DataFrame, dict[str, pd.DataFrame]]: + """Read a kcat source file and map its rows to the ``model.ec`` reactions. + + 1. A row with ``rxn`` and ``enzymes`` applies to the ``model.ec`` rows + with that conventional reaction and that enzyme set. + 2. A row with only ``enzymes`` applies to every row with that enzyme + set; a row with only ``rxn`` to every row of that reaction (curated + rows; :func:`save_kcat_source` always writes both). + 3. A row whose kcat is missing, not a number, infinite, or not above 0 + is reported as invalid and not applied. + 4. A row that matches nothing is reported, never dropped silently. If a + ``model.ec`` row has at least half of the row's enzymes, that partial + match is reported as well, but not applied. + + The enzymes of a row may be listed in any order. + + Parameters + ---------- + model + The ecModel. + name + Name of the source (``data/kcats/.tsv`` in the project folder), + or the path to a kcat source file. + + Returns + ------- + (kcat_list, report) + ``kcat_list`` has one row per applied match, in the kcat-list schema + (``rxn_id``, ``source``, ``eccode``, ``substrates``, ``genes``, + ``kcat``, ``wildcard_level``, ``origin``) plus ``log10_sd`` and + ``line`` (the line of the file). ``report`` has the ``invalid``, + ``unmatched`` and ``partial`` rows. + """ + path = Path(name) + if not path.is_file(): + path = _project_path(model, adapter) / "data" / "kcats" / f"{name}.tsv" + if not path.is_file(): + raise FileNotFoundError( + f"No kcat source file {path}. Make it with the kcats_from_* function of this " + f"source (kcats_from_brenda, kcats_from_dlkcat, kcats_from_okp), or with " + f"save_kcat_source." + ) + records, meta = read_kcat_source(path) + base, enzymes = ec_rxn_keys(model) + + by_pair: dict[tuple[str, str], list[int]] = defaultdict(list) + by_base: dict[str, list[int]] = defaultdict(list) + by_enz: dict[str, list[int]] = defaultdict(list) + by_acc: dict[str, list[int]] = defaultdict(list) + acc_sets = [set(e.split(" + ")) if e else set() for e in enzymes] + for i, (b, e) in enumerate(zip(base, enzymes)): + by_pair[(b, e)].append(i) + by_base[b].append(i) + if e: + by_enz[e].append(i) + for a in acc_sets[i]: + by_acc[a].append(i) + + kcats = records["kcat"].to_numpy(dtype=float) + invalid = [r for r in range(len(records)) if not (math.isfinite(kcats[r]) and kcats[r] > 0)] + invalid_set = set(invalid) + rec_idx, ec_idx, unmatched, partial = [], [], [], [] + for r, (rxn, enz_text) in enumerate(zip(records["rxn"], records["enzymes"])): + if r in invalid_set: + continue + enz = _normalize_enzymes(enz_text) if enz_text else "" + if rxn and enz: + hit = by_pair.get((rxn, enz), []) + elif enz: + hit = by_enz.get(enz, []) + elif rxn: + hit = by_base.get(rxn, []) + else: + hit = [] + if not hit: + unmatched.append(r) + if enz: + wanted = enz.split(" + ") + candidates = sorted({c for a in wanted for c in by_acc.get(a, [])}) + if rxn: + candidates = [c for c in candidates if base[c] == rxn] + for c in candidates: + if np.mean([w in acc_sets[c] for w in wanted]) >= 0.5: + partial.append((r, c)) + else: + rec_idx.extend([r] * len(hit)) + ec_idx.extend(hit) + + sub = records.iloc[rec_idx] + wc, origin = [], [] + for d in sub["detail"]: + m = re.search(r"wc=(\d+)", d or "") + wc.append(float(m.group(1)) if m else np.nan) + m = re.search(r"origin=(\d+)", d or "") + origin.append(float(m.group(1)) if m else np.nan) + kcat_list = pd.DataFrame({ + "rxn_id": [model.ec.rxns[i] for i in ec_idx], + "source": sub["source"].tolist(), + "eccode": sub["ec_code"].tolist(), + "substrates": [[s.strip() for s in t.split("; ") if s.strip()] for t in sub["substrate"]], + "genes": [[] for _ in rec_idx], + "kcat": sub["kcat"].to_numpy(dtype=float), + "wildcard_level": pd.Series(wc, dtype=float), + "origin": pd.Series(origin, dtype=float), + "log10_sd": sub["log10_sd"].to_numpy(dtype=float), + "line": pd.Series(sub["line"].tolist(), dtype=int), + }) + kcat_list.attrs["source_name"] = meta.get("source", "") + + def rows_of(idx): + return pd.DataFrame({ + "line": records["line"].iloc[idx].tolist(), + "rxn": records["rxn"].iloc[idx].tolist(), + "enzymes": records["enzymes"].iloc[idx].tolist(), + "kcat": records["kcat"].iloc[idx].tolist(), + }) + + report = { + "invalid": rows_of(invalid), + "unmatched": rows_of(unmatched), + "partial": pd.DataFrame({ + "line": [int(records["line"].iloc[r]) for r, _ in partial], + "rxn": [records["rxn"].iloc[r] for r, _ in partial], + "enzymes": [records["enzymes"].iloc[r] for r, _ in partial], + "ec_rxn": [model.ec.rxns[c] for _, c in partial], + "ec_enzymes": [enzymes[c] for _, c in partial], + }), + } + if invalid: + logger.warning("%d of %d rows of %s have no valid kcat; see report['invalid'].", + len(invalid), len(records), path) + if unmatched: + logger.warning("%d of %d rows of %s apply to no reaction of this model; see report['unmatched'].", + len(unmatched), len(records), path) + return kcat_list, report diff --git a/src/geckopy/gather_kcats/kcats_from.py b/src/geckopy/gather_kcats/kcats_from.py new file mode 100644 index 0000000..7d29b47 --- /dev/null +++ b/src/geckopy/gather_kcats/kcats_from.py @@ -0,0 +1,365 @@ +"""Gather the kcats of one source and keep them as a kcat source file. + +One function per source, each with its own options, each writing +``data/kcats/.tsv`` in the project folder by default (see +:mod:`geckopy.gather_kcats.kcat_sources`). GECKO (MATLAB) has the same +functions as ``kcatsFromBrenda``, ``kcatsFromDLKcat`` and ``kcatsFromOKP``. +""" +from __future__ import annotations + +import logging +import lzma +import re +import shutil +from datetime import datetime, timezone +from pathlib import Path +from typing import TYPE_CHECKING, Iterable, Optional + +import pandas as pd + +from ..databases import load_brenda_data, load_dlkcat_ignore_lists, load_phyl_dist +from .fuzzy_kcat_matching import fuzzy_kcat_matching +from .kcat_sources import kcat_source_file, save_kcat_source +from .open_kinetics_predictor import fetch_open_kinetics_predictor, submit_open_kinetics_predictor +from .open_kinetics_predictor.parse_output import parse_okp_output +from .read_dlkcat_output import read_dlkcat_output +from .run_dlkcat import run_dlkcat +from .write_dlkcat_input import write_dlkcat_input + +if TYPE_CHECKING: + from ..adapter import ModelAdapter + from ..databases import BrendaData, PhylDist + from ..ec_model.ec_model import EcModel + +logger = logging.getLogger(__name__) + + +def _adapter(model: "EcModel", adapter: Optional["ModelAdapter"], purpose: str) -> "ModelAdapter": + from ..adapter import resolve_adapter + + return resolve_adapter(model, adapter, purpose=purpose) + + +def _rxn_subset(ec_rxns: Optional[Iterable[str]]) -> Optional[list[str]]: + if ec_rxns is None: + return None + ec_rxns = list(ec_rxns) + if not ec_rxns: + raise ValueError("ec_rxns is empty; leave it out to gather for all reactions.") + return ec_rxns + + +def _warn_subset(ec_rxns: Optional[list[str]], model: "EcModel", path: Path) -> None: + if ec_rxns is not None and len(set(ec_rxns)) < len(model.ec.rxns) and path.is_file(): + logger.warning("%s is replaced by kcats for %d reactions only. Give save_as to keep it.", + path, len(ec_rxns)) + + +def _brenda_release(folder: Path) -> str: + """The release line at the top of kcat.tsv, e.g. 'BRENDA release 2026.1'.""" + for path, opener in ((folder / "kcat.tsv.xz", lzma.open), (folder / "kcat.tsv", open)): + if path.is_file(): + with opener(path, "rt", encoding="utf-8") as f: + m = re.match(r"^#\s*(BRENDA release \S+)", f.readline()) + return m.group(1) if m else "" + return "" + + +def kcats_from_brenda( + model: "EcModel", + *, + aggregation: Optional[str] = None, + force_wildcard_level: int = 0, + ec_rxns: Optional[Iterable[str]] = None, + brenda: Optional["BrendaData"] = None, + phyl_dist: Optional["PhylDist"] = None, + save: bool = True, + save_as: str = "brenda", + adapter: Optional["ModelAdapter"] = None, +) -> pd.DataFrame: + """Gather kcats from BRENDA by fuzzy matching and keep them as ``brenda.tsv``. + + Runs locally on the BRENDA files of geckopy (see + :func:`~geckopy.gather_kcats.fuzzy_kcat_matching`). + + Parameters + ---------- + model + EcModel with EC numbers. + aggregation + Which BRENDA value per EC number, substrate and organism is used: + ``"max"`` or ``"median"`` (default: ``params.kcat_aggregate_brenda``). + force_wildcard_level + Force a minimum wildcard level in the EC numbers. + ec_rxns + ``model.ec.rxns`` to gather kcats for (default: all). With the + default ``save_as`` the file then only holds those reactions. + brenda, phyl_dist + Loaded BRENDA data and KEGG phylogenetic distances (default: loaded + from the adapter's folders). + save, save_as + Write the source file, and its name without ``.tsv``. + adapter + Model adapter (default: ``model.adapter``). + + Returns + ------- + pandas.DataFrame + The kcat list. + """ + adapter = _adapter(model, adapter, "kcats_from_brenda needs the BRENDA folder and params") + if aggregation is None: + aggregation = adapter.params.kcat_aggregate_brenda + aggregation = aggregation.lower() + if aggregation not in ("max", "median"): + raise ValueError("aggregation must be 'max' or 'median'") + ec_rxns = _rxn_subset(ec_rxns) + if save: + path = kcat_source_file(save_as, model=model, adapter=adapter) # fails before the matching + _warn_subset(ec_rxns, model, path) + folder = Path(adapter.get_brenda_db_folder()) + if brenda is None: + brenda = load_brenda_data(folder) + if phyl_dist is None: + phyl_dist = load_phyl_dist(adapter.get_phyl_dist_path()) + kcat_list = fuzzy_kcat_matching( + model, brenda, phyl_dist, ec_rxns=ec_rxns, + force_wildcard_level=force_wildcard_level, aggregate=aggregation, adapter=adapter, + ) + if save: + parameters = {"aggregation": aggregation, "force_wildcard_level": force_wildcard_level} + if ec_rxns is not None: + parameters["ec_rxns"] = len(ec_rxns) + path = save_kcat_source( + model, kcat_list, save_as, adapter=adapter, generator="kcats_from_brenda", + parameters=parameters, data=_brenda_release(folder), + ) + logger.info("BRENDA kcats for %d reactions written to %s", + int((kcat_list["kcat"] > 0).sum()), path) + return kcat_list + + +def kcats_from_dlkcat( + model: "EcModel", + *, + run: bool = False, + ec_rxns: Optional[Iterable[str]] = None, + save: bool = True, + save_as: str = "dlkcat", + adapter: Optional["ModelAdapter"] = None, +) -> pd.DataFrame: + """Gather kcats from DLKcat and keep them as ``dlkcat.tsv``. + + With ``run=True`` the DLKcat input is written and DLKcat runs through + its Docker image (which needs Docker); otherwise the output of an + earlier run in ``data/DLKcat.tsv`` is read. A run records how many + reactions it covered in ``data/DLKcat_run.txt``, so that the source file + says so also when the output is read later. + + Parameters + ---------- + model + EcModel with metabolite SMILES and protein sequences. + run + Write the input and run DLKcat (default False). + ec_rxns + ``model.ec.rxns`` to run DLKcat for (default: all); needs ``run``. + save, save_as, adapter + As in :func:`kcats_from_brenda`. + + Returns + ------- + pandas.DataFrame + The kcat list. + """ + adapter = _adapter(model, adapter, "kcats_from_dlkcat needs params.path") + ec_rxns = _rxn_subset(ec_rxns) + if ec_rxns is not None and not run: + raise ValueError("ec_rxns selects the reactions DLKcat is run for, so it needs run=True.") + if save: + path = kcat_source_file(save_as, model=model, adapter=adapter) + if run: + _warn_subset(ec_rxns, model, path) + data_dir = Path(adapter.params.path) / "data" + dlkcat_file = data_dir / "DLKcat.tsv" + run_file = data_dir / "DLKcat_run.txt" + if run: + if shutil.which("docker") is None: + raise RuntimeError( + "DLKcat runs in a Docker image, but Docker was not found. Install it from " + "https://www.docker.com, or call kcats_from_dlkcat without run=True to read " + "an earlier DLKcat output." + ) + write_dlkcat_input(model, dlkcat_file, load_dlkcat_ignore_lists(data_dir), + ec_rxns=ec_rxns, overwrite=True) + run_dlkcat(dlkcat_file) + run_file.write_text(f"ecRxns: {len(set(ec_rxns))}\n" if ec_rxns is not None else "", + encoding="utf-8", newline="\n") + elif not dlkcat_file.is_file(): + raise FileNotFoundError( + f"No DLKcat output at {dlkcat_file}. Run kcats_from_dlkcat(model, run=True) first." + ) + kcat_list = read_dlkcat_output(model, dlkcat_file) + if save: + parameters = {"run": bool(run)} + if run_file.is_file(): + m = re.search(r"ecRxns:\s*(\d+)", run_file.read_text(encoding="utf-8")) + if m: + parameters["ec_rxns"] = int(m.group(1)) + path = save_kcat_source( + model, kcat_list, save_as, adapter=adapter, generator="kcats_from_dlkcat", + parameters=parameters, + ) + logger.info("DLKcat kcats written to %s", path) + return kcat_list + + +_JOB_KEYS = ("jobId", "method", "submittedAt", "ecRxns", "resultOf") +_DATABASE_SOURCES = {"okp-brenda", "okp-sabio-rk", "okp-uniprot"} + + +def _read_job_file(path: Path) -> dict[str, str]: + job: dict[str, str] = {} + if path.is_file(): + for line in path.read_text(encoding="utf-8").splitlines(): + m = re.match(r"^([A-Za-z]\w*):\s*(.*?)\s*$", line) + if m: + job[m.group(1)] = m.group(2) + return job + + +def _write_job_file(path: Path, job: dict[str, str]) -> None: + path.write_text("".join(f"{k}: {job[k]}\n" for k in _JOB_KEYS if k in job), + encoding="utf-8", newline="\n") + + +def _adopt_shared_result(model: "EcModel", data_dir: Path, method: str, job_file: Path, out_file: Path) -> None: + """Take over a result that ``fetch_open_kinetics_predictor`` stored by + itself, in ``data/OKP_output.csv`` with ``data/OKP_job.txt``, only when + its job is for this method and its predictions are labelled with this + method.""" + shared = data_dir / "OKP_output.csv" + job = _read_job_file(data_dir / "OKP_job.txt") + ok = shared.is_file() and "jobId" in job and job.get("method", "").lower() == method.lower() + if ok: + sources = [str(x) for x in parse_okp_output(model, shared)["source"]] + predicted = [x for x in sources if x.lower().startswith("okp-") and x.lower() not in _DATABASE_SOURCES] + ok = bool(predicted) and all(x.lower() == f"okp-{method}".lower() for x in predicted) + if not ok: + raise FileNotFoundError( + f"No stored OKP result for {method} ({out_file}). Call kcats_from_okp without " + f"use_stored to submit or fetch a job." + ) + shutil.copyfile(shared, out_file) + _write_job_file(job_file, {**job, "method": method, "resultOf": job["jobId"]}) + + +def kcats_from_okp( + model: "EcModel", + *, + method: Optional[str] = None, + wait: bool = False, + resubmit: bool = False, + use_stored: bool = False, + ec_rxns: Optional[Iterable[str]] = None, + save: bool = True, + save_as: Optional[str] = None, + adapter: Optional["ModelAdapter"] = None, +) -> tuple[Optional[pd.DataFrame], bool]: + """Gather kcats from OpenKineticsPredictor and keep them as ``okp_.tsv``. + + OKP is a web service with a job queue, so this can be called repeatedly: + without a job for ``method`` it submits one; with a job it fetches the + result, and writes the file once the job is done. Each method keeps its + own job in ``data/OKP_job_.txt`` and its own result in + ``data/OKP_output_.csv``, so jobs of several methods can run at + the same time. The job file records which job the stored result came + from; ``use_stored`` only uses a result of the current job. An API key + is needed. + + Parameters + ---------- + model + EcModel with metabolite SMILES and protein sequences. + method + OKP method, e.g. ``"CataPro"``, ``"CatRange"``, ``"UniKP"`` (default: + ``params.okp.method``). + wait + Wait until the job is done. + resubmit + Submit a new job even if one exists for this method; its earlier + result is removed. + use_stored + Read the stored result of this method's job without contacting OKP. + ec_rxns + ``model.ec.rxns`` to predict for (default: all); only used when a + job is submitted. + save, save_as, adapter + As in :func:`kcats_from_brenda`; ``save_as`` defaults to + ``okp_``, in lower case, with other characters than letters + and digits as ``_``. + + Returns + ------- + (kcat_list, done) + The kcat list (``None`` while the job is not done), and whether the + result was available. + """ + adapter = _adapter(model, adapter, "kcats_from_okp needs params.path") + if not method: + method = adapter.params.okp.method + ec_rxns = _rxn_subset(ec_rxns) + tag = re.sub(r"[^A-Za-z0-9]+", "_", method).lower() + data_dir = Path(adapter.params.path) / "data" + job_file = data_dir / f"OKP_job_{tag}.txt" + out_file = data_dir / f"OKP_output_{tag}.csv" + name = save_as or f"okp_{tag}" + job = _read_job_file(job_file) + submitting = not use_stored and (resubmit or "jobId" not in job) + if ec_rxns is not None and not submitting: + raise ValueError( + f"ec_rxns selects the reactions of a new job, but this call submits none (a job for " + f"{method} exists; use resubmit=True for a new one)." + ) + if save: + path = kcat_source_file(name, model=model, adapter=adapter) # fails before submitting + if submitting: + _warn_subset(ec_rxns, model, path) + + if use_stored: + if "jobId" not in job: + _adopt_shared_result(model, data_dir, method, job_file, out_file) + job = _read_job_file(job_file) + if not out_file.is_file() or job.get("resultOf") != job.get("jobId"): + raise FileNotFoundError( + f"No stored result of the current {method} job ({job_file}). Call kcats_from_okp " + f"without use_stored to fetch it." + ) + kcat_list, done = parse_okp_output(model, out_file), True + else: + if submitting: + out_file.unlink(missing_ok=True) + job_id = submit_open_kinetics_predictor(model, overwrite=True, method=method, + ec_rxns=ec_rxns, adapter=adapter) + job = {"jobId": job_id, "method": method, + "submittedAt": datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ")} + if ec_rxns is not None: + job["ecRxns"] = str(len(set(ec_rxns))) + _write_job_file(job_file, job) + done, kcat_list = fetch_open_kinetics_predictor(model, job_id=job["jobId"], wait=wait, + adapter=adapter, out_file=out_file) + if done: + job["resultOf"] = job["jobId"] + _write_job_file(job_file, job) + + if not done: + logger.info("Call kcats_from_okp(model, method=%r) again later to collect the result.", method) + return None, False + if save: + parameters: dict = {"method": method} + if "ecRxns" in job: + parameters["ec_rxns"] = int(job["ecRxns"]) + path = save_kcat_source(model, kcat_list, name, adapter=adapter, generator="kcats_from_okp", + parameters=parameters, data=f"OKP job {job['jobId']}") + logger.info("OKP (%s) kcats written to %s", method, path) + return kcat_list, True diff --git a/src/geckopy/gather_kcats/open_kinetics_predictor/run.py b/src/geckopy/gather_kcats/open_kinetics_predictor/run.py index 3ecd283..2643d5f 100644 --- a/src/geckopy/gather_kcats/open_kinetics_predictor/run.py +++ b/src/geckopy/gather_kcats/open_kinetics_predictor/run.py @@ -141,6 +141,7 @@ def fetch_open_kinetics_predictor( adapter: Optional["ModelAdapter"] = None, base_url: str = _BASE_URL, session: Optional[requests.Session] = None, + out_file: Optional[str | Path] = None, ) -> tuple[bool, Optional[pd.DataFrame]]: """Check an OKP job and, when finished, download + parse the result. @@ -159,6 +160,9 @@ def fetch_open_kinetics_predictor( Polling cadence and overall deadline (seconds) when ``wait``. api_key, adapter, base_url, session Key resolution / adapter / testing knobs (see submit). + out_file + Where the result is stored, and read from with ``use_stored`` + (default: ``data/OKP_output.csv``). Returns ------- @@ -173,7 +177,7 @@ def fetch_open_kinetics_predictor( model, adapter, purpose="fetch_open_kinetics_predictor needs params.path", ) data_dir = Path(adapter.params.path) / "data" - out_file = data_dir / "OKP_output.csv" + out_file = Path(out_file) if out_file is not None else data_dir / "OKP_output.csv" if use_stored: if not out_file.is_file(): diff --git a/tests/data/kcatSourceFixture.tsv b/tests/data/kcatSourceFixture.tsv new file mode 100644 index 0000000..4cfdc9f --- /dev/null +++ b/tests/data/kcatSourceFixture.tsv @@ -0,0 +1,16 @@ +# source: fixture +# generator: hand-written +"# parameters: {""aggregation"":""median""}" +# kind: curated + +rxn enzymes substrate kcat log10_sd source detail ec_code note +R2 P2 + P1 m1 12.5 0.3 brenda wc=0; origin=1 1.1.1.1 enzymes in another order +# a comment between rows +R3 P4 "2,3-bisphospho-D-glycerate" 1e-05 custom quoted substrate + +R5 P5 m2 5,5 custom decimal comma: invalid +R2_REV m2 7 custom reaction only + P3 3 custom enzymes only +R9 X1 -1 custom negative: invalid +R9 X1 4 custom unmatched + diff --git a/tests/test_adapter.py b/tests/test_adapter.py index 5bbb093..f7244e2 100644 --- a/tests/test_adapter.py +++ b/tests/test_adapter.py @@ -134,3 +134,120 @@ def test_shipped_example_adapter_toml_is_valid(example_dir): """Verify every folder under examples/ has a model_adapter.toml that validates against the current ModelParameters schema.""" ModelAdapter.from_folder(example_dir) + + +SPONTANEOUS_TOML = MINIMAL_TOML + """ +[spontaneous] +name_contains = "spontaneous" +file = "data/spontaneousReactions.tsv" +""" + + +def test_spontaneous_reactions_from_toml(tmp_path): + import cobra + + folder = _write_adapter_folder(tmp_path, SPONTANEOUS_TOML) + (folder / "data").mkdir() + (folder / "data" / "spontaneousReactions.tsv").write_text( + "rxns\tspontaneous\nR1\t0\nR2\t1\nR3\t1\n" + ) + model = cobra.Model("m") + for rid, name in [("R1", "a"), ("R2", "b"), ("R4", "spontaneous c"), ("R5", "d")]: + model.add_reactions([cobra.Reaction(rid, name=name)]) + adapter = ModelAdapter.from_folder(folder) + assert adapter.params.spontaneous.file == folder / "data" / "spontaneousReactions.tsv" + assert adapter.get_spontaneous_reactions(model) == ["R2", "R4"] + + +# A file as GECKO's convertAdapterToToml writes it, MATLAB-only settings included. +CONVERTED_TOML = """\ +# Model adapter for testus testus, converted from TestGEMAdapter. +conv_gem = "models/testModel.xml" +org_name = "testus testus" +sigma = 0.5 +p_tot = 0.5 +f = 4 +gr_exp = 0.41 +c_source = "E1" +bio_rxn = "R4" +enzyme_comp = "c" + +[kegg] +id = "tst" +gene_id = "kegg" + +[uniprot] +type = "taxonomy" +id = "TE00000000" +gene_id_field = "gene_tst" +reviewed = true + +[complex] +taxonomic_id = 123456 + +[evotune.source_groups.dlkcat] +sources = ["dlkcat"] + +[evotune.source_groups.brenda] +sources = ["brenda"] + +[evotune.sigma0_log_source] +dlkcat = 0.4 +brenda = 0.2 + +[matlab] +make_model_anaerobic = "anaerobicModel_GECKO" + +[matlab.params] +c_UptakeExp = 0.641339301 +""" + + +def test_toml_converted_from_matlab_adapter_loads(tmp_path): + folder = _write_adapter_folder(tmp_path, CONVERTED_TOML) + adapter = ModelAdapter.from_folder(folder) + assert adapter.params.f == 4 + assert adapter.params.matlab.make_model_anaerobic == "anaerobicModel_GECKO" + assert adapter.params.matlab.params == {"c_UptakeExp": 0.641339301} + assert set(adapter.params.evotune.source_groups) == {"dlkcat", "brenda"} + + +def _schema_leaves(model_cls, prefix=""): + """Dotted keys of every leaf field of a pydantic model, with defaults.""" + from pydantic import BaseModel + + leaves = {} + instance = model_cls.model_construct() + for name, field in model_cls.model_fields.items(): + key = f"{prefix}{name}" + annotation = field.annotation + if isinstance(annotation, type) and issubclass(annotation, BaseModel): + leaves.update(_schema_leaves(annotation, key + ".")) + else: + leaves[key] = field + return leaves + + +def test_adapter_parameter_table_matches_schema(): + """adapterParameters.tsv, shared with GECKO, lists every key and its default.""" + import csv + import tomllib + + from geckopy.adapter.params import ModelParameters + + table = Path(__file__).resolve().parent.parent / "src" / "geckopy" / "adapter" / "adapterParameters.tsv" + with open(table, newline="", encoding="utf-8") as f: + rows = {r["toml_key"]: r for r in csv.DictReader(f, delimiter="\t", quoting=csv.QUOTE_NONE)} + + leaves = _schema_leaves(ModelParameters) + leaves.pop("path") # set from the project folder, never written + leaves.pop("matlab.params") # free-form + assert set(rows) == set(leaves) + + for key, field in leaves.items(): + default_text = rows[key]["default"] + if field.is_required() or not default_text: + continue + expected = tomllib.loads(f"x = {default_text}")["x"] + default = field.get_default(call_default_factory=True) + assert default == expected, key 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/tests/test_kcat_sources.py b/tests/test_kcat_sources.py new file mode 100644 index 0000000..0cb2a11 --- /dev/null +++ b/tests/test_kcat_sources.py @@ -0,0 +1,406 @@ +"""Tests for kcat source files and the kcats_from_* functions.""" +from pathlib import Path + +import cobra +import numpy as np +import pandas as pd +import pytest +from scipy import sparse + +from geckopy import EcModel, ModelAdapter +from geckopy.ec_model.ec_data import EcData +from geckopy.gather_kcats import ( + ec_rxn_keys, + kcats_from_brenda, + kcats_from_dlkcat, + kcats_from_okp, + load_kcat_source, + read_kcat_source, + save_kcat_source, +) +from geckopy.gather_kcats import kcats_from as kcats_from_module +from geckopy.gather_kcats.kcat_sources import format_number + +# Conventional reactions with their isozymes (each a set of enzymes). +REACTIONS = { + "R1": [["P2", "P1"], ["P3"]], + "R1_REV": [["P2", "P1"], ["P3"]], + "R2": [["P4"]], +} +ENZYMES = ["P1", "P2", "P3", "P4"] + + +def _adapter(tmp_path: Path) -> ModelAdapter: + (tmp_path / "model_adapter.toml").write_text('conv_gem = "dummy.xml"\norg_name = "testus testus"\n', encoding="utf-8") + return ModelAdapter.from_folder(tmp_path) + + +def _model(adapter: ModelAdapter, light: bool) -> EcModel: + """Full model: R1_EXP_1, R1_EXP_2, ...; light model: 001_R1, 002_R1, ...""" + model = EcModel("test", adapter=adapter, gecko_light=light) + met = cobra.Metabolite("m1", compartment="c", name="m1") + model.add_metabolites([met]) + rxns, rows = [], [] + for base, isozymes in REACTIONS.items(): + for n, enzymes in enumerate(isozymes, start=1): + if light: + rxns.append(f"{n:03d}_{base}") + else: + rxns.append(f"{base}_EXP_{n}" if len(isozymes) > 1 else base) + rows.append([ENZYMES.index(e) for e in enzymes]) + if not light: + for n in range(len(isozymes)): + rid = f"{base}_EXP_{n + 1}" if len(isozymes) > 1 else base + model.add_reactions([cobra.Reaction(rid)]) + model.reactions.get_by_id(rid).add_metabolites({met: -1}) + else: + model.add_reactions([cobra.Reaction(base)]) + mat = sparse.lil_matrix((len(rxns), len(ENZYMES))) + for i, cols in enumerate(rows): + for j in cols: + mat[i, j] = 1 + n = len(rxns) + model.ec = EcData( + gecko_light=light, rxns=rxns, kcat=np.zeros(n), source=[""] * n, notes=[""] * n, + eccodes=["1.1.1.1"] * n, genes=[f"G{k}" for k in range(len(ENZYMES))], enzymes=list(ENZYMES), + mw=np.full(len(ENZYMES), 40000.0), sequence=["M"] * len(ENZYMES), + concs=np.full(len(ENZYMES), np.nan), rxn_enz_mat=mat.tocsr(), + ) + return model + + +def _kcat_list(model: EcModel) -> pd.DataFrame: + n = len(model.ec.rxns) + return pd.DataFrame({ + "rxn_id": model.ec.rxns, + "source": ["brenda"] * n, + "eccode": ["1.1.1.1"] * n, + "substrates": [["m1", "ATP"]] * n, + "genes": [[] for _ in range(n)], + "kcat": [10.0 * (k + 1) for k in range(n)], + "wildcard_level": [0] * n, + "origin": [k % 6 + 1 for k in range(n)], + }) + + +def test_ec_rxn_keys_full_and_light(tmp_path): + adapter = _adapter(tmp_path) + base_f, enz_f = ec_rxn_keys(_model(adapter, light=False)) + base_l, enz_l = ec_rxn_keys(_model(adapter, light=True)) + assert base_f[:2] == ["R1", "R1"] and base_f[2] == "R1_REV" + assert enz_f[0] == "P1 + P2" + assert sorted(zip(base_f, enz_f)) == sorted(zip(base_l, enz_l)) + + +def test_save_read_load_round_trip(tmp_path): + adapter = _adapter(tmp_path) + full = _model(adapter, light=False) + kcats = _kcat_list(full) + path = save_kcat_source(full, kcats, "brenda", parameters={"aggregation": "max"}, data="BRENDA release test") + assert path == tmp_path / "data" / "kcats" / "brenda.tsv" + records, meta = read_kcat_source(path) + assert meta["source"] == "brenda" and meta["kind"] == "gathered" + assert meta["parameters"] == '{"aggregation":"max"}' + assert records["substrate"].iloc[0] == "m1; ATP" + assert records["detail"].iloc[0] == "wc=0; origin=1" + loaded, report = load_kcat_source(full, "brenda") + assert report["unmatched"].empty + merged = loaded.merge(kcats, on="rxn_id", suffixes=("", "_in")) + assert len(merged) == len(kcats) + assert (merged["kcat"] == merged["kcat_in"]).all() + assert (merged["origin"] == merged["origin_in"]).all() + + +def test_full_file_applies_to_light_model(tmp_path): + adapter = _adapter(tmp_path) + full, light = _model(adapter, light=False), _model(adapter, light=True) + save_kcat_source(full, _kcat_list(full), "brenda") + loaded, report = load_kcat_source(light, "brenda") + assert report["unmatched"].empty + assert sorted(loaded["rxn_id"]) == sorted(light.ec.rxns) + + +def test_matching_rules_and_curated_file(tmp_path): + adapter = _adapter(tmp_path) + full = _model(adapter, light=False) + folder = tmp_path / "data" / "kcats" + folder.mkdir(parents=True) + (folder / "custom.tsv").write_text( + "# source: custom\n# kind: curated\n" + "rxn\tenzymes\tsubstrate\tkcat\tlog10_sd\tsource\tdetail\tec_code\tnote\n" + "\tP2 + P1\t\t10\t\tcustom\t\t\tenzymes only, other order\n" + "R2\t\t\t20\t\tcustom\t\t\treaction only\n" + "R999\tX1\t\t30\t\tcustom\t\t\tunmatched\n" + "R1\tP1 + P2 + XNEW\t\t40\t\tcustom\t\t\tpartial\n" + , encoding="utf-8") + loaded, report = load_kcat_source(full, "custom") + assert sorted(loaded.loc[loaded["line"] == 4, "rxn_id"]) == ["R1_EXP_1", "R1_REV_EXP_1"] + assert list(loaded.loc[loaded["line"] == 5, "rxn_id"]) == ["R2"] + assert list(report["unmatched"]["line"]) == [6, 7] + assert set(report["partial"]["line"]) == {7} + # Curated files are never written: by name, 'kind: curated' and + # 'source: custom', in any case; the file is left as it was. + before = (folder / "custom.tsv").read_bytes() + with pytest.raises(PermissionError): + save_kcat_source(full, _kcat_list(full), "custom") + assert (folder / "custom.tsv").read_bytes() == before + (folder / "mine.tsv").write_text("# source: Custom\nrxn\tkcat\tsource\n", encoding="utf-8") + (folder / "mine2.tsv").write_text("# kind: Curated\nrxn\tkcat\tsource\n", encoding="utf-8") + for name in ("mine", "mine2"): + with pytest.raises(PermissionError): + save_kcat_source(full, _kcat_list(full), name) + + +@pytest.mark.parametrize("value, text", [ + (10.0, "10"), (0.1, "0.1"), (1.5e-05, "1.5e-05"), (123.456, "123.456"), + (np.nan, ""), (2.0 / 3.0, "0.6666666666666666"), (np.inf, "Inf"), (-np.inf, "-Inf"), + (-3.0, "-3"), (np.int64(7), "7"), (999999999999999.0, "999999999999999"), (1e15, "1e+15"), + (2e20, "2e+20"), (0.30000000000000004, "0.30000000000000004"), +]) +def test_format_number_matches_gecko(value, text): + assert format_number(value) == text + + +def test_kcats_from_dlkcat_and_okp_write_files(tmp_path, monkeypatch): + adapter = _adapter(tmp_path) + full = _model(adapter, light=False) + data = tmp_path / "data" + data.mkdir() + (data / "DLKcat.tsv").write_text("placeholder\n", encoding="utf-8") + dlkcat = _kcat_list(full).assign(source="DLKcat", wildcard_level=np.nan, origin=np.nan) + monkeypatch.setattr(kcats_from_module, "read_dlkcat_output", lambda model, path: dlkcat) + kcats_from_dlkcat(full) + records, meta = read_kcat_source(data / "kcats" / "dlkcat.tsv") + assert meta["generator"].startswith("kcats_from_dlkcat") + assert set(records["source"]) == {"DLKcat"} + + # A result stored by fetch_open_kinetics_predictor alone is taken over + # when its job is for the method and its predictions are of that method. + (data / "OKP_job.txt").write_text("jobId: abc123\nmethod: CataPro\n", encoding="utf-8") + (data / "OKP_output.csv").write_text("placeholder\n", encoding="utf-8") + okp = _kcat_list(full).assign(source="OKP-CataPro", wildcard_level=np.nan, origin=np.nan) + monkeypatch.setattr(kcats_from_module, "parse_okp_output", lambda model, path: okp) + kcat_list, done = kcats_from_okp(full, method="CataPro", use_stored=True) + assert done + assert (data / "OKP_output_catapro.csv").is_file() + records, meta = read_kcat_source(data / "kcats" / "okp_catapro.tsv") + assert meta["data"] == "OKP job abc123" + assert meta["parameters"] == '{"method":"CataPro"}' + # The shared job names UniKP, but the shared result holds CataPro + # predictions: refused, and nothing is written. + (data / "OKP_job.txt").write_text("jobId: uni1\nmethod: UniKP\n", encoding="utf-8") + with pytest.raises(FileNotFoundError): + kcats_from_okp(full, method="UniKP", use_stored=True) + assert not (data / "kcats" / "okp_unikp.tsv").exists() + with pytest.raises(ValueError): + kcats_from_dlkcat(full, ec_rxns=["R2"]) + with pytest.raises(ValueError): + kcats_from_okp(full, ec_rxns=[]) + + +def test_okp_jobs_per_method(tmp_path, monkeypatch): + """The sequence that mixed up methods: two jobs in flight, one fetched, + then use_stored for the other. Each method keeps its own job and result.""" + adapter = _adapter(tmp_path) + full = _model(adapter, light=False) + data = tmp_path / "data" + data.mkdir() + results = {"job-CataPro": "OKP-CataPro", "job-UniKP": "OKP-UniKP"} + finished = set() + calls = [] + + def fake_submit(model, **kw): + calls.append(("submit", kw["method"], kw.get("ec_rxns"))) + (data / "OKP_job.txt").write_text(f"jobId: job-{kw['method']}\nmethod: {kw['method']}\n", encoding="utf-8") + return f"job-{kw['method']}" + + def fake_fetch(model, **kw): + calls.append(("fetch", kw["job_id"])) + if kw["job_id"] not in finished: + return False, None + Path(kw["out_file"]).write_text(results[kw["job_id"]], encoding="utf-8") + (data / "OKP_output.csv").write_text(results[kw["job_id"]], encoding="utf-8") + return True, _kcat_list(model).assign(source=results[kw["job_id"]], wildcard_level=np.nan, origin=np.nan) + + def fake_parse(model, path): + return _kcat_list(model).assign(source=Path(path).read_text(), wildcard_level=np.nan, origin=np.nan) + + monkeypatch.setattr(kcats_from_module, "submit_open_kinetics_predictor", fake_submit) + monkeypatch.setattr(kcats_from_module, "fetch_open_kinetics_predictor", fake_fetch) + monkeypatch.setattr(kcats_from_module, "parse_okp_output", fake_parse) + + assert kcats_from_okp(full, method="CataPro") == (None, False) + assert kcats_from_okp(full, method="UniKP", ec_rxns=["R2"]) == (None, False) + finished.add("job-CataPro") + assert kcats_from_okp(full, method="CataPro")[1] + with pytest.raises(FileNotFoundError): + kcats_from_okp(full, method="UniKP", use_stored=True) + assert not (data / "kcats" / "okp_unikp.tsv").exists() + records, _ = read_kcat_source(data / "kcats" / "okp_catapro.tsv") + assert set(records["source"]) == {"OKP-CataPro"} + + finished.add("job-UniKP") + kcats_from_okp(full, method="UniKP") + records, meta = read_kcat_source(data / "kcats" / "okp_unikp.tsv") + assert set(records["source"]) == {"OKP-UniKP"} + assert meta["parameters"] == '{"ec_rxns":1,"method":"UniKP"}' + assert meta["data"] == "OKP job job-UniKP" + + # ec_rxns without a new job is an error; a resubmit removes the old result. + with pytest.raises(ValueError): + kcats_from_okp(full, method="UniKP", ec_rxns=["R2"]) + finished.discard("job-UniKP") + assert kcats_from_okp(full, method="UniKP", resubmit=True) == (None, False) + assert not (data / "OKP_output_unikp.csv").exists() + with pytest.raises(FileNotFoundError): + kcats_from_okp(full, method="UniKP", use_stored=True) + assert calls[0] == ("submit", "CataPro", None) and calls[2] == ("submit", "UniKP", ["R2"]) + + +def test_kcats_from_brenda_writes_file(tmp_path, monkeypatch): + adapter = _adapter(tmp_path) + full = _model(adapter, light=False) + seen = {} + + def fake_fuzzy(model, brenda, phyl, **kw): + seen.update(kw) + return _kcat_list(model) + + monkeypatch.setattr(kcats_from_module, "fuzzy_kcat_matching", fake_fuzzy) + kcats_from_brenda(full, aggregation="Median", brenda=object(), phyl_dist=object()) + assert seen["adapter"] is adapter and seen["aggregate"] == "median" + records, meta = read_kcat_source(tmp_path / "data" / "kcats" / "brenda.tsv") + assert meta["parameters"] == '{"aggregation":"median","force_wildcard_level":0}' + assert len(records) == len(full.ec.rxns) + + +def test_parameters_json_matches_gecko(): + """The same text as GECKO's testGuardsNumbersAndJson_tc0078.""" + from geckopy.gather_kcats.kcat_sources import _parameters_json + + value = {"z": 1.0, "a": {"y": "café", "b": True}, "n": np.nan, "f": 1 / 3, "g": 1e-7, + "h": 1e15, "k": 1e20, "l": [5], "s": 'tab\tquote"'} + assert _parameters_json(value) == ( + '{"a":{"b":true,"y":"café"},"f":0.3333333333333333,"g":1e-07,"h":1e+15,' + '"k":1e+20,"l":[5],"n":null,"s":"tab\\tquote\\"","z":1}' + ) + assert _parameters_json({"x": np.int64(3), "y": np.array([1.5, 2.0]), "i": np.inf}) == \ + '{"i":null,"x":3,"y":[1.5,2]}' + + +def test_enzymeless_rows_are_not_written(tmp_path): + adapter = _adapter(tmp_path) + full = _model(adapter, light=False) + mat = full.ec.rxn_enz_mat.tolil() + i = full.ec.rxns.index("R1_EXP_2") + mat[i, :] = 0 + full.ec.rxn_enz_mat = mat.tocsr() + save_kcat_source(full, _kcat_list(full), "brenda") + loaded, _ = load_kcat_source(full, "brenda") + assert "R1_EXP_2" not in set(loaded["rxn_id"]) + assert (loaded["rxn_id"] == "R1_EXP_1").sum() == 1 + + +def test_rows_are_sorted(tmp_path): + adapter = _adapter(tmp_path) + full = _model(adapter, light=False) + kl = _kcat_list(full).iloc[::-1] + path = save_kcat_source(full, kl, "brenda") + records, _ = read_kcat_source(path) + keys = list(zip(records["rxn"], records["enzymes"], records["substrate"], records["source"])) + assert keys == sorted(keys) + assert b"\r\n" not in path.read_bytes() + + +def test_shared_fixture(): + """The fixture file, read the same way by GECKO: byte-order mark, CRLF, + a quoted header line and field, comments and blank lines between rows, + a row of tabs, a decimal comma, a negative kcat.""" + path = Path(__file__).resolve().parent / "data" / "kcatSourceFixture.tsv" + records, meta = read_kcat_source(path) + assert meta["source"] == "fixture" + assert meta["parameters"] == '{"aggregation":"median"}' + assert meta["kind"] == "curated" + assert list(records["line"]) == [7, 9, 11, 12, 13, 14, 15] + assert list(records["rxn"]) == ["R2", "R3", "R5", "R2_REV", "", "R9", "R9"] + assert records["substrate"].iloc[1] == "2,3-bisphospho-D-glycerate" + kcats = records["kcat"].tolist() + assert kcats[:2] == [12.5, 1e-05] and np.isnan(kcats[2]) and kcats[3:] == [7, 3, -1, 4] + assert records["log10_sd"].iloc[0] == 0.3 + assert records["kcat"].dtype == float + + +def test_extra_fields_are_an_error(tmp_path): + path = tmp_path / "bad.tsv" + path.write_text("rxn\tkcat\tsource\nR2\t1\tbrenda\textra\n", encoding="utf-8") + with pytest.raises(ValueError, match="more fields"): + read_kcat_source(path) + + +def test_guards_numbers_and_parity_of_reading(tmp_path): + adapter = _adapter(tmp_path) + full = _model(adapter, light=False) + folder = tmp_path / "data" / "kcats" + folder.mkdir(parents=True) + # The curated guard reads quoted header lines and ignores the case of keys. + (folder / "a.tsv").write_text('"# kind: curated"\nrxn\tkcat\tsource\n', encoding="utf-8") + (folder / "b.tsv").write_text("# Kind: Curated\nrxn\tkcat\tsource\n", encoding="utf-8") + (folder / "c.tsv").write_text("# SOURCE: custom\nrxn\tkcat\tsource\n", encoding="utf-8") + for name in ("a", "b", "c"): + with pytest.raises(PermissionError): + save_kcat_source(full, _kcat_list(full), name) + # The source functions check the file before they gather. + with pytest.raises(PermissionError): + kcats_from_brenda(full, save_as="custom", brenda=object(), phyl_dist=object()) + + # Infinite, NaN and non-positive kcats are not written, and not applied. + kl = _kcat_list(full) + kl.loc[0, "kcat"], kl.loc[1, "kcat"], kl.loc[2, "kcat"] = np.inf, np.nan, -2.0 + path = save_kcat_source(full, kl, "numbers", data="two\nlines") + records, meta = read_kcat_source(path) + assert (records["kcat"] > 0).all() and np.isfinite(records["kcat"]).all() + assert meta["data"] == "two lines" + (folder / "inf.tsv").write_text("rxn\tenzymes\tkcat\tsource\nR2\tP4\tInf\tcustom\nR2\tP4\t1e999\tcustom\n" + "R2\tP4\t\u0661\u0662\tcustom\nR2\tP4\t\" 12 \"\tcustom\n", encoding="utf-8") + loaded, report = load_kcat_source(full, folder / "inf.tsv") + assert list(report["invalid"]["line"]) == [2, 3, 4] + assert list(loaded["kcat"]) == [12.0] # quoted with spaces: trimmed, as in GECKO + + # Substrates come back as lists; whole-number kcats as float. + loaded, _ = load_kcat_source(full, path) + assert loaded["substrates"].iloc[0] == ["m1", "ATP"] + assert loaded["kcat"].dtype == float + + # A column name twice, or a column named line, is an error. + (folder / "dup.tsv").write_text("rxn\tkcat\tkcat\tsource\n", encoding="utf-8") + (folder / "line.tsv").write_text("rxn\tkcat\tsource\tline\n", encoding="utf-8") + for name in ("dup", "line"): + with pytest.raises(ValueError): + read_kcat_source(folder / f"{name}.tsv") + + +def test_rows_sorted_by_every_column(tmp_path): + """The same order as GECKO's testGuardsNumbersAndJson_tc0078.""" + adapter = _adapter(tmp_path) + full = _model(adapter, light=False) + folder = tmp_path / "data" / "kcats" + folder.mkdir(parents=True) + (folder / "unsorted.tsv").write_text( + "rxn\tenzymes\tsubstrate\tkcat\tsource\tdetail\n" + "R2\tP4\tm1\t2\tb\t\nR2\tP4\tm1\t2\ta\twc=1\nR2\tP4\t\t5\ta\t\nR1_REV\tP3\tm2\t1\ta\t\n", encoding="utf-8") + loaded, _ = load_kcat_source(full, folder / "unsorted.tsv") + path = save_kcat_source(full, loaded, "sorted") + records, _ = read_kcat_source(path) + assert list(records["rxn"]) == ["R1_REV", "R2", "R2", "R2"] + assert list(records["substrate"]) == ["m2", "", "m1", "m1"] + assert list(records["source"]) == ["a", "a", "a", "b"] + + +def test_subset_warns(tmp_path, monkeypatch, caplog): + adapter = _adapter(tmp_path) + full = _model(adapter, light=False) + monkeypatch.setattr(kcats_from_module, "fuzzy_kcat_matching", lambda model, b, p, **kw: _kcat_list(model)) + kcats_from_brenda(full, brenda=object(), phyl_dist=object()) + with caplog.at_level("WARNING"): + kcats_from_brenda(full, ec_rxns=["R2"], brenda=object(), phyl_dist=object()) + assert "replaced by kcats for 1 reactions only" in caplog.text 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.