Kcat strategies, BRENDA from toolbox-assets and GECKO 3 project upgrades

.gitignore

@@ -17,6 +17,8 @@ geckopy/geckopy/data_files/*.xml
 *.ppt*
 # Only allowed .mat file
 !phylDist.mat
+# BRENDA files, downloaded by fetchBrendaAsset
+databases/brenda-*/
 # Ignore intermediate model files
 *_stage*.yml
 

databases/README.md

@@ -8,21 +8,17 @@
 
 `DLKcatIgnoreMets.tsv` is a table of small metabolites/ions that `writeDLKcatInput` filters out as DLKcat does not predict kcat values for such substrates. This can either be in this folder (and a pull request to the GitHub repository will make this more widely available to other users), or by keeping a copy of this file in the `data` subfolder of the model adapter folder.
 
-### k<sub>cat</sub> values from BRENDA
+### BRENDA
 
-`kcat.tsv` collates kcat values per organism, reaction and substrate, as gathered from the BRENDA database. Both a per-triple maximum (`kcat_max`) and median (`kcat_median`) are recorded. Which one is used is set by `params.kcatAggregation` on the model adapter, or per call with `fuzzyKcatMatching(model,'kcatAggregation','median')`; the default is `max`. Produced by the `geckopy brenda-refresh` CLI in the [geckopy](https://github.com/CompMetEng/geckopy) Python package, which maintains the actively updated BRENDA-parsing tooling for both toolboxes; `mw.tsv` and `sa.tsv` below come from the same tool and share its `#`-prefixed release-version header line.
+GECKO does not ship BRENDA data. `fetchBrendaAsset` downloads the files of the [toolbox-assets](https://github.com/CompMetEng/toolbox-assets) release `brenda-2026.1-1` ([BRENDA](https://www.brenda-enzymes.org) 2026.1, CC BY 4.0) on first use into `databases/brenda-2026.1-1`, checked against the manifest checksums. geckopy uses the same release. A project can use its own BRENDA files instead, in `data/brenda`, written by the `geckopy brenda-refresh` CLI of [geckopy](https://github.com/CompMetEng/geckopy).
 
-### Molecular weights
+- `kcat.tsv`: kcat values per organism, reaction and substrate. Both a per-triple maximum (`kcat_max`) and median (`kcat_median`) are recorded. Which one is used is set by `params.kcatAggregation` on the model adapter, or per call with `fuzzyKcatMatching(model,'kcatAggregation','median')`; the default is `max`.
+- `sa.tsv`: specific activities per organism, reaction and substrate, with a maximum (`sa_max`) and median (`sa_median`); `loadBRENDAdata` uses `sa_max`.
+- `mw.tsv`: molecular weights per organism and reaction (without explicitly referring to a protein identifier).
+- `measurements.tsv`: one row per wild-type k<sub>cat</sub> measurement with an EC number: EC number, substrate (name and InChIKey), organism, UniProt accession and sequence id, k<sub>cat</sub> (1/s; the geometric midpoint of a range), temperature, pH and PubMed id. Values more than four orders of magnitude below their EC class, and clusters of values below 0.01 /s that sit three orders of magnitude under the other values of their EC number, are left out as unit errors. `kcatsFromECTree` reads it (see `loadBrendaMeasurements`).
+- `sequences.fasta`, `substrates.tsv` and `kcatCombineParams.tsv`: the sequences and substrates of the measurements, and the parameters of the `combine` stage.
 
-`mw.tsv` is a collation of molecular weights per organism and reaction (without explicitly referring to a protein identifier), as gathered from the BRENDA database.
-
-### Specific activities from BRENDA
-
-`sa.tsv` collates specific activities per organism, reaction and substrate, as gathered from the BRENDA database. Like `kcat.tsv`, both a maximum (`sa_max`) and median (`sa_median`) are recorded; `loadBRENDAdata` currently uses `sa_max` only.
-
-### BRENDA measurements
-
-`measurements.tsv` holds one row per wild-type k<sub>cat</sub> measurement in [BRENDA](https://www.brenda-enzymes.org) with an EC number: EC number, substrate (name and InChIKey), organism, UniProt accession and sequence id, k<sub>cat</sub> (1/s; the geometric midpoint of a range), temperature, pH and PubMed id. Values more than four orders of magnitude below their EC class, and clusters of values below 0.01 /s that sit three orders of magnitude under the other values of their EC number, are left out as unit errors. `kcatsFromECTree` reads it (see `loadBrendaMeasurements`); geckopy ships the same file. The first line gives the BRENDA release.
+Each file starts with a `#` line naming the BRENDA release.
 
 ### Taxonomic tree
 

src/change_model/upgradeEcModel.m

@@ -0,0 +1,57 @@
+function [model, changed] = upgradeEcModel(model)
+% upgradeEcModel  Bring an ecModel made with GECKO 3 to the GECKO 4 format.
+%
+% GECKO 3 defined the protein usage reactions (usage_prot_<id>, prot_<id>
+% -> prot_pool) and the protein pool exchange (prot_pool_exchange,
+% prot_pool ->) in the reverse direction, used with negative flux. GECKO 4
+% defines them forward (prot_pool -> prot_<id>, -> prot_pool), so that an
+% enzyme's usage is the positive flux of its usage reaction. This turns
+% each reaction that is still reverse: its column of S, its bounds and its
+% objective coefficient change sign. Reactions already forward are left as
+% they are, so calling it again changes nothing. Flux values of the model
+% (fluxes of usage reactions and the pool exchange) change sign with it.
+% geckopy has the same (upgrade_ec_model).
+%
+% Parameters
+% ----------
+% model : struct
+%     an ecModel.
+%
+% Returns
+% -------
+% model : struct
+%     the ecModel in the GECKO 4 format.
+% changed : cell
+%     the reactions that were turned.
+%
+% Examples
+% --------
+%     ecModel = upgradeEcModel(loadEcModel('ecYeastGEM_GECKO3.yml'));
+%
+% See also
+% --------
+% upgradeGECKOproject, makeEcModel, setProtPoolSize
+
+pool = find(strcmp(model.mets, 'prot_pool'));
+changed = {};
+if isempty(pool)
+    return
+end
+usage = find(startsWith(model.rxns, 'usage_prot_'));
+reverse = usage(full(model.S(pool, usage)) > 0);
+exchange = find(strcmp(model.rxns, 'prot_pool_exchange'));
+if ~isempty(exchange) && full(model.S(pool, exchange)) < 0
+    reverse(end+1, 1) = exchange;
+end
+if isempty(reverse)
+    return
+end
+model.S(:, reverse) = -model.S(:, reverse);
+lb = model.lb(reverse);
+model.lb(reverse) = -model.ub(reverse);
+model.ub(reverse) = -lb;
+if isfield(model, 'c')
+    model.c(reverse) = -model.c(reverse);
+end
+changed = model.rxns(reverse);
+end

src/download_databases/fetchBrendaAsset.m

@@ -0,0 +1,62 @@
+function folder = fetchBrendaAsset(files)
+% fetchBrendaAsset  Download (once) GECKO's BRENDA files and return their folder.
+%
+% GECKO does not ship BRENDA data. The files of the toolbox-assets release
+% brenda-2026.1-1 (BRENDA 2026.1, CC BY 4.0) are downloaded on first use
+% with RAVEN's fetchToolboxAsset, which checks them against the checksums
+% of the toolbox-assets manifest, and are unpacked into
+% databases/brenda-2026.1-1 in the GECKO folder. geckopy reads the same
+% release (fetch_brenda_asset). A project can use its own BRENDA files
+% instead, in data/brenda (see ModelAdapter.getBrendaDBFolder).
+%
+% Parameters
+% ----------
+% files : cell
+%     file names (default: kcat.tsv, sa.tsv, mw.tsv and measurements.tsv);
+%     the release also has sequences.fasta, substrates.tsv and
+%     kcatCombineParams.tsv.
+%
+% Returns
+% -------
+% folder : char
+%     the folder with the files, unpacked.
+%
+% Examples
+% --------
+%     folder = fetchBrendaAsset();
+%     folder = fetchBrendaAsset({'kcatCombineParams.tsv'});
+%
+% See also
+% --------
+% loadBRENDAdata, loadBrendaMeasurements, ModelAdapter
+
+release = 'brenda-2026.1-1';
+if nargin < 1 || isempty(files)
+    files = {'kcat.tsv', 'sa.tsv', 'mw.tsv', 'measurements.tsv'};
+end
+files = cellstr(files);
+packed = struct('kcat_tsv', 1, 'sa_tsv', 1, 'mw_tsv', 1, 'measurements_tsv', 1, ...
+    'sequences_fasta', 1, 'substrates_tsv', 1, 'kcatCombineParams_tsv', 0);
+folder = fullfile(findGECKOroot(), 'databases', release);
+if exist('fetchToolboxAsset', 'file') ~= 2
+    error('fetchBrendaAsset:raven', ['GECKO downloads its BRENDA files with RAVEN''s ' ...
+        'fetchToolboxAsset, which this RAVEN does not have. Update RAVEN to version 3.'])
+end
+for i = 1:numel(files)
+    name = files{i};
+    key = matlab.lang.makeValidName(name);
+    if ~isfield(packed, key)
+        error('fetchBrendaAsset:file', '%s is not a file of %s.', name, release)
+    end
+    target = fullfile(folder, name);
+    if isfile(target)
+        continue
+    end
+    if packed.(key)
+        gz = fetchToolboxAsset(folder, release, [name '.gz']);
+        gunzip(gz, folder);
+    else
+        fetchToolboxAsset(folder, release, name);
+    end
+end
+end