Estimate kcats from the EC tree; make brenda_ec_tree the default preset

databases/README.md

@@ -20,6 +20,10 @@
 
 `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.
+
 ### Taxonomic tree
 
 `PhylDist.mat` is a taxonomic tree of KEGG organisms, as generated by RAVEN Toolbox.

src/gather_kcats/collectKcats.m

@@ -9,7 +9,7 @@
 % The order is a preset or a list of stages:
 %
 %     [kcat_assignment]
-%     preset = "gecko_median"      # the default, or gecko_classic
+%     preset = "brenda_ec_tree"    # the default; or gecko_median, gecko_classic
 %     isozymes = true              # fill from isozymes after the stages
 %     standard = true              # standard kcat for the rest
 %
@@ -26,20 +26,24 @@
 % listed source with a value (see mergeKcats). Sources are the names of
 % source files (custom, dlkcat, okp_catapro, ...) and the BRENDA quality
 % tiers database_exact, database_top and database_bottom, which read
-% brenda.tsv. Both presets order custom, database_top, dlkcat and
-% database_bottom:
-% - gecko_median (the default) 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 updateKcatSources).
+% brenda.tsv. The presets:
+% - brenda_ec_tree (the default): custom kcats, the BRENDA median of
+%   measurements of the organism on the reaction's substrate (origin 1,
+%   without EC wildcards), the EC tree (kcatsFromECTree), then DLKcat for
+%   reactions without an EC number.
+% - gecko_median: custom kcats, the BRENDA median of matches without EC
+%   wildcards from the organism or other organisms, with or without the
+%   substrate (origin 1 to 4), DLKcat, then the other BRENDA matches (EC
+%   wildcards, and specific activities, origin 5 and 6).
+% - gecko_classic: as gecko_median with the BRENDA maximum and the
+%   specific-activity matches above DLKcat (origin up to 6), as the GECKO 3
+%   tutorial.
+% Against measured kcats hidden from the search, the EC tree was closer
+% than the BRENDA median at origins 2 and 4 and with EC wildcards and as
+% close at origins 1 and 3; the BRENDA median was as close as DLKcat or
+% closer at every origin, and the maximum further off than the median. The
+% BRENDA aggregation is [kcat_sources.brenda] aggregation (see
+% updateKcatSources).
 %
 % After the stages:
 % - custom kcats of data/customKcats.tsv are applied (applyCustomKcats),

src/gather_kcats/ecTreeEstimate.m

@@ -0,0 +1,122 @@
+function [log10kcat, n, depth, ecUsed] = ecTreeEstimate(eccodes, measurements, k)
+% ecTreeEstimate  EC-tree kcat estimates for EC numbers.
+%
+% For each EC number, the median log10 kcat of its measurements is pulled
+% toward the estimate of its parent class with weight n / (n + k), where n
+% is the number of measurements: the global median, then the 2-digit
+% class, the 3-digit class and the full EC number, in that order. A
+% well-measured EC keeps its own median; a sparsely measured one leans on
+% its class. geckopy has the same (geckopy.gather_kcats.ec_tree.ECTree).
+%
+% Parameters
+% ----------
+% eccodes : cell
+%     per query one or more EC numbers separated by ';'. The one with most
+%     measurements of its full EC is used (the first on a tie).
+% measurements : struct
+%     BRENDA measurement table with fields ec and log10kcat (see
+%     loadBrendaMeasurements).
+% k : double
+%     shrinkage constant (default 5): the number of measurements at which
+%     an EC's own median and its class count equally.
+%
+% Returns
+% -------
+% log10kcat : double
+%     the estimates; the global median for queries without measurements.
+% n : double
+%     number of measurements of the full EC number used.
+% depth : double
+%     deepest EC level with measurements: 4, 3, 2, or 0 when there are
+%     none (or no EC number).
+% ecUsed : cell
+%     the EC number used.
+%
+% Examples
+% --------
+%     m = loadBrendaMeasurements();
+%     [v, n, depth] = ecTreeEstimate({'1.1.1.1'; '2.7.1.1;2.7.1.2'}, m);
+%     kcat = 10.^v;
+%
+% See also
+% --------
+% kcatsFromECTree, loadBrendaMeasurements
+
+if nargin < 3 || isempty(k)
+    k = 5;
+end
+if ~(isnumeric(k) && isscalar(k) && k > 0)
+    error('ecTreeEstimate:k', 'k must be a positive number.')
+end
+if ischar(eccodes) || isstring(eccodes)
+    eccodes = cellstr(eccodes);
+end
+ec = measurements.ec(:);
+v = measurements.log10kcat(:);
+globalMedian = median(v);
+lv = cell(4, 1);
+for d = [2 3 4]
+    if d == 4
+        key = ec;
+    else
+        key = cellfun(@(e) parentEC(e, d), ec, 'UniformOutput', false);
+    end
+    [keys, ~, g] = unique(key);
+    lv{d} = struct('keys', {keys}, 'median', accumarray(g, v, [], @median), ...
+        'n', accumarray(g, 1));
+end
+
+nq = numel(eccodes);
+log10kcat = repmat(globalMedian, nq, 1);
+n = zeros(nq, 1);
+depth = zeros(nq, 1);
+ecUsed = repmat({''}, nq, 1);
+for q = 1:nq
+    ecs = strtrim(strsplit(char(eccodes{q}), ';'));
+    ecs = ecs(~cellfun(@isempty, ecs));
+    if isempty(ecs)
+        continue
+    end
+    best = ecs{1};
+    bestN = -1;
+    for e = 1:numel(ecs)
+        c = countOf(lv{4}, ecs{e});
+        if c > bestN
+            best = ecs{e};
+            bestN = c;
+        end
+    end
+    value = globalMedian;
+    for d = [2 3 4]
+        if d == 4
+            key = best;
+        else
+            key = parentEC(best, d);
+        end
+        [c, med] = countOf(lv{d}, key);
+        if c > 0
+            w = c / (c + k);
+            value = w * med + (1 - w) * value;
+            depth(q) = d;
+        end
+    end
+    log10kcat(q) = value;
+    n(q) = bestN;
+    ecUsed{q} = best;
+end
+end
+
+function p = parentEC(ec, d)
+parts = strsplit(ec, '.');
+p = strjoin(parts(1:min(d, numel(parts))), '.');
+end
+
+function [c, med] = countOf(level, key)
+[found, i] = ismember(key, level.keys);
+c = 0;
+med = NaN;
+if found
+    c = level.n(i);
+    med = level.median(i);
+end
+end

src/gather_kcats/kcatAssignmentConfig.m

@@ -4,7 +4,7 @@
 % Reads [kcat_assignment] and [kcat_sources.*] of model_adapter.toml (or
 % params.kcat_assignment and params.kcat_sources of an adapter class), and
 % expands a preset into its stages, as collectKcats and updateKcatSources
-% use them. Without a preset or stages the preset is gecko_median.
+% use them. Without a preset or stages the preset is brenda_ec_tree.
 % geckopy's kcat_assignment_config gives the same settings.
 %
 % Parameters
@@ -48,7 +48,7 @@
 if hasStages
     config.stages = arrayfun(@readStage, a.stage(:)');
 else
-    preset = 'gecko_median';
+    preset = 'brenda_ec_tree';
     if hasPreset
         preset = char(a.preset);
     end
@@ -75,9 +75,13 @@
     case 'gecko_median'
         stages = select(order, 4);
         sources.brenda = struct('aggregation', 'median');
+    case 'brenda_ec_tree'
+        stages = select({'custom', 'database_top', 'ec_tree', 'dlkcat'}, 1);
+        sources.brenda = struct('aggregation', 'median');
+        sources.ec_tree = struct('k', 5);
     otherwise
         error('kcatAssignment:config', 'Unknown kcat assignment preset "%s". Presets: %s.', ...
-            preset, 'gecko_classic, gecko_median')
+            preset, 'gecko_classic, gecko_median, brenda_ec_tree')
 end
 end
 

src/gather_kcats/kcatsFromECTree.m

@@ -0,0 +1,136 @@
+function kcatList = kcatsFromECTree(model, varargin)
+% kcatsFromECTree  Estimate kcats from the EC tree and keep them as a source file.
+%
+% Each reaction gets the median kcat of all wild-type BRENDA measurements
+% of its EC number, from any organism and substrate, shrunk toward the
+% median of the EC's 3-digit and 2-digit classes with weight n / (n + k)
+% for n measurements (see ecTreeEstimate). A reaction with several EC
+% numbers uses the one with most measurements. Reactions without an EC
+% number, or whose EC classes have no measurements, get no value. The
+% result is written to data/kcats/ec_tree.tsv in the project folder (see
+% saveKcatSource).
+%
+% Parameters
+% ----------
+% model : struct
+%     an ecModel in GECKO 3 format (with ecModel.ec structure), with EC
+%     numbers (getECfromDatabase or getECfromGEM).
+%
+% Name-Value Arguments
+% --------------------
+% k : double
+%     shrinkage constant (default 5): the number of measurements at which
+%     an EC's own median and its class count equally.
+% ecRxns : logical
+%     for which model.ec.rxns kcats are estimated (default: all).
+% measurements : struct
+%     BRENDA measurement table (default: loadBrendaMeasurements with the
+%     model adapter).
+% save : logical
+%     write the source file (default true).
+% saveAs : char
+%     name of the source file, without .tsv (default 'ec_tree'). With
+%     ecRxns, the file then only holds those reactions.
+% modelAdapter : ModelAdapter
+%     a loaded model adapter (default: the current default model adapter).
+%
+% Returns
+% -------
+% kcatList : struct
+%     the kcat values, with fields source, rxns, eccodes, kcats and detail
+%     (the number of measurements of the EC, n, and the deepest EC level
+%     with measurements, depth).
+%
+% Examples
+% --------
+%     kcatsFromECTree(ecModel);
+%     kcatList = kcatsFromECTree(ecModel, 'k', 10, 'save', false);
+%
+% See also
+% --------
+% ecTreeEstimate, loadBrendaMeasurements, kcatsFromBrenda, saveKcatSource, assignKcats
+
+checkOptions(varargin, {'k','ecRxns','measurements','save','saveAs','modelAdapter'}, 'kcatsFromECTree');
+p = parseGECKOargs(varargin, { ...
+    'k',            5; ...
+    'ecRxns',       []; ...
+    'measurements', []; ...
+    'save',         true; ...
+    'saveAs',       'ec_tree'; ...
+    'modelAdapter', []});
+modelAdapter = p.modelAdapter;
+if isempty(modelAdapter)
+    modelAdapter = ModelAdapterManager.getDefault();
+    if isempty(modelAdapter)
+        error('Either send in a modelAdapter or set the default model adapter in the ModelAdapterManager.')
+    end
+end
+if ~isfield(model.ec, 'eccodes')
+    error('No EC codes defined in model.ec.eccodes. First run getECfromGEM() and/or getECfromDatabase().')
+end
+if p.save
+    file = fullfile(modelAdapter.params.path, 'data', 'kcats', [char(p.saveAs) '.tsv']);
+    warnSubset(p.ecRxns, numel(model.ec.rxns), file);
+end
+measurements = p.measurements;
+if isempty(measurements)
+    measurements = loadBrendaMeasurements('modelAdapter', modelAdapter);
+end
+
+idx = (1:numel(model.ec.rxns))';
+if ~isempty(p.ecRxns)
+    idx = idx(p.ecRxns(:));
+end
+eccodes = model.ec.eccodes(idx);
+eccodes(cellfun(@(e) ~ischar(e) && ~isstring(e), eccodes)) = {''};
+[values, n, depth, ecUsed] = ecTreeEstimate(eccodes, measurements, p.k);
+found = depth > 0;
+idx = idx(found);
+kcatList.source  = 'ec_tree';
+kcatList.rxns    = model.ec.rxns(idx);
+kcatList.eccodes = ecUsed(found);
+kcatList.kcats   = 10 .^ values(found);
+kcatList.detail  = cellstr("n=" + string(n(found)) + "; depth=" + string(depth(found)));
+
+if p.save
+    parameters = struct('k', p.k);
+    if ~isempty(p.ecRxns)
+        parameters.ec_rxns = countRxns(p.ecRxns);
+    end
+    saveKcatSource(model, kcatList, p.saveAs, 'modelAdapter', modelAdapter, ...
+        'generator', 'kcatsFromECTree', 'parameters', parameters, ...
+        'data', measurements.release);
+    fprintf('EC-tree kcats for %d reactions written to %s\n', numel(kcatList.rxns), file);
+end
+end
+
+function n = countRxns(ecRxns)
+if islogical(ecRxns)
+    n = sum(ecRxns);
+else
+    n = numel(unique(ecRxns));
+end
+end
+
+function warnSubset(ecRxns, nRxns, file)
+% A subset written under the default name replaces the full file.
+if ~isempty(ecRxns) && countRxns(ecRxns) < nRxns && isfile(file)
+    printOrange(sprintf(['%s is replaced by kcats for %d reactions only. Give saveAs ' ...
+        'to keep it.\n'], file, countRxns(ecRxns)));
+end
+end
+
+function checkOptions(args, names, functionName)
+% Options are name-value pairs with known names, so a misspelled option is
+% an error instead of being taken as a positional argument.
+id = [functionName ':arguments'];
+given = args(1:2:end);
+if mod(numel(args), 2) ~= 0 || ~all(cellfun(@(x) ischar(x) || (isstring(x) && isscalar(x)), given))
+    error(id, 'Options of %s are name-value pairs. Options: %s.', functionName, strjoin(names, ', '))
+end
+unknown = setdiff(cellfun(@char, given, 'UniformOutput', false), names);
+if ~isempty(unknown)
+    error(id, 'Unknown option of %s: %s. Options: %s.', functionName, strjoin(unknown, ', '), ...
+        strjoin(names, ', '))
+end
+end