Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions .gitattributes
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@

# A fixture that must keep its byte-order mark and CRLF line endings.
test/unit_tests/ecTestGEM/data/kcatSourceFixture.tsv -text
11 changes: 9 additions & 2 deletions src/gather_kcats/fetchOpenKineticsPredictor.m
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,9 @@
% return (default false).
% pollInterval : double
% seconds between polls when wait=true (default 30).
% outFile : char
% where the result is stored, and read from with useStored=true
% (default: data/OKP_output.csv).
% modelAdapter : ModelAdapter
% a loaded model adapter (default: the current default model adapter).
%
Expand Down Expand Up @@ -56,7 +59,8 @@
'jobId', []; ...
'wait', []; ...
'pollInterval', []; ...
'modelAdapter', []});
'modelAdapter', []; ...
'outFile', []});
useStored = p.useStored;
jobId = p.jobId;
wait = p.wait;
Expand All @@ -74,7 +78,10 @@
if isempty(wait); wait = false; end
if isempty(pollInterval); pollInterval = 30; end

outFile = fullfile(params.path,'data','OKP_output.csv');
outFile = p.outFile;
if isempty(outFile)
outFile = fullfile(params.path,'data','OKP_output.csv');
end

%% Stored-result mode: parse the saved file, no API call
if useStored
Expand Down
51 changes: 51 additions & 0 deletions src/gather_kcats/getEcRxnKeys.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,51 @@
function [baseRxns, enzymeSets] = getEcRxnKeys(model)
% getEcRxnKeys Build-independent keys of the model.ec reactions.
%
% A kcat source file (see saveKcatSource) identifies a reaction by its
% conventional reaction and its enzyme set, not by model.ec.rxns, whose
% identifiers depend on how the ecModel was built. This function gives
% both keys for every model.ec row.
%
% Parameters
% ----------
% model : struct
% an ecModel in GECKO 3 format (with ecModel.ec structure).
%
% Returns
% -------
% baseRxns : cell
% per model.ec.rxns entry, the conventional reaction: without a
% trailing _EXP_<n> (full model) or a leading <nnn>_ (light model). A
% reverse direction keeps its _REV.
% enzymeSets : cell
% per model.ec.rxns entry, the UniProt accessions of its enzymes
% (model.ec.rxnEnzMat), sorted and joined by ' + '; empty without
% enzymes.
%
% Examples
% --------
% [baseRxns, enzymeSets] = getEcRxnKeys(ecModel);
%
% See also
% --------
% saveKcatSource, loadKcatSource

if model.ec.geckoLight
baseRxns = regexprep(model.ec.rxns, '^\d+_', '');
else
baseRxns = regexprep(model.ec.rxns, '_EXP_\d+$', '');
end
baseRxns = baseRxns(:);

n = numel(model.ec.rxns);
enzymeSets = repmat({''}, n, 1);
[rows, cols] = find(model.ec.rxnEnzMat);
if ~isempty(rows)
% Enzymes per reaction in one pass, rather than a search per reaction.
perRxn = accumarray(rows(:), cols(:), [n 1], @(c) {c});
for i = find(~cellfun(@isempty, perRxn))'
enzymes = unique(model.ec.enzymes(perRxn{i}));
enzymeSets{i} = strjoin(enzymes(:)', ' + ');
end
end
end
128 changes: 128 additions & 0 deletions src/gather_kcats/kcatsFromBrenda.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,128 @@
function kcatList = kcatsFromBrenda(model, varargin)
% kcatsFromBrenda Gather kcats from BRENDA and keep them as a source file.
%
% Matches the model's EC numbers and substrates to BRENDA by fuzzy
% matching (see fuzzyKcatMatching), and writes the result to
% data/kcats/brenda.tsv in the project folder (see saveKcatSource). The
% matching runs locally on the BRENDA files of GECKO.
%
% Parameters
% ----------
% model : struct
% an ecModel in GECKO 3 format (with ecModel.ec structure), with EC
% numbers (getECfromDatabase or getECfromGEM).
%
% Name-Value Arguments
% --------------------
% aggregation : char
% which BRENDA value per EC number, substrate and organism is used:
% 'max' or 'median' (default: params.kcatAggregation of the model
% adapter, or 'max').
% forceWClvl : double
% force a minimum wildcard level in the EC numbers (default 0).
% ecRxns : logical
% for which model.ec.rxns kcats are gathered (default: all).
% save : logical
% write the source file (default true).
% saveAs : char
% name of the source file, without .tsv (default 'brenda'). 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, as returned by fuzzyKcatMatching.
%
% Examples
% --------
% kcatsFromBrenda(ecModel);
% kcatList = kcatsFromBrenda(ecModel, 'aggregation', 'median');
%
% See also
% --------
% fuzzyKcatMatching, kcatsFromDLKcat, kcatsFromOKP, saveKcatSource, loadKcatSource

checkOptions(varargin, {'aggregation','forceWClvl','ecRxns','save','saveAs','modelAdapter'}, 'kcatsFromBrenda');
p = parseGECKOargs(varargin, { ...
'aggregation', []; ...
'forceWClvl', 0; ...
'ecRxns', []; ...
'save', true; ...
'saveAs', 'brenda'; ...
'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
aggregation = resolveKcatAggregation(p.aggregation, modelAdapter);
if p.save
file = fullfile(modelAdapter.params.path, 'data', 'kcats', [char(p.saveAs) '.tsv']);
warnSubset(p.ecRxns, numel(model.ec.rxns), file);
end

kcatList = fuzzyKcatMatching(model, 'ecRxns', p.ecRxns, 'modelAdapter', modelAdapter, ...
'forceWClvl', p.forceWClvl, 'kcatAggregation', aggregation);

if p.save
parameters = struct('aggregation', aggregation, 'force_wildcard_level', p.forceWClvl);
if ~isempty(p.ecRxns)
parameters.ec_rxns = countRxns(p.ecRxns);
end
saveKcatSource(model, kcatList, p.saveAs, 'modelAdapter', modelAdapter, ...
'generator', 'kcatsFromBrenda', 'parameters', parameters, ...
'data', brendaRelease(modelAdapter));
fprintf('BRENDA kcats for %d reactions written to %s\n', sum(kcatList.kcats > 0), file);
end
end

function release = brendaRelease(modelAdapter)
% The release line at the top of kcat.tsv, e.g. 'BRENDA release 2026.1'.
release = '';
file = fullfile(modelAdapter.getBrendaDBFolder(), 'kcat.tsv');
fid = fopen(file, 'r');
if fid < 0
return
end
line = fgetl(fid);
fclose(fid);
tok = regexp(line, '^#\s*(BRENDA release [^ ]+)', 'tokens', 'once');
if ~isempty(tok)
release = tok{1};
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
145 changes: 145 additions & 0 deletions src/gather_kcats/kcatsFromDLKcat.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,145 @@
function kcatList = kcatsFromDLKcat(model, varargin)
% kcatsFromDLKcat Gather kcats from DLKcat and keep them as a source file.
%
% Predicts kcats with DLKcat and writes them to data/kcats/dlkcat.tsv in
% the project folder (see saveKcatSource). With 'run', true it writes the
% DLKcat input (writeDLKcatInput) and runs DLKcat through its Docker image
% (runDLKcat), which needs Docker; otherwise it reads the output of an
% earlier run (readDLKcatOutput). 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 : struct
% an ecModel in GECKO 3 format (with ecModel.ec structure), with
% metabolite SMILES (findMetSmiles) and protein sequences.
%
% Name-Value Arguments
% --------------------
% run : logical
% write the input and run DLKcat; false reads the existing output in
% data/DLKcat.tsv (default false).
% ecRxns : logical
% for which model.ec.rxns DLKcat is run (default: all). Only used with
% 'run', true.
% save : logical
% write the source file (default true).
% saveAs : char
% name of the source file, without .tsv (default 'dlkcat').
% modelAdapter : ModelAdapter
% a loaded model adapter (default: the current default model adapter).
%
% Returns
% -------
% kcatList : struct
% the kcat values, as returned by readDLKcatOutput.
%
% Examples
% --------
% kcatsFromDLKcat(ecModel, 'run', true);
% kcatList = kcatsFromDLKcat(ecModel);
%
% See also
% --------
% writeDLKcatInput, runDLKcat, readDLKcatOutput, kcatsFromBrenda, saveKcatSource

checkOptions(varargin, {'run','ecRxns','save','saveAs','modelAdapter'}, 'kcatsFromDLKcat');
p = parseGECKOargs(varargin, { ...
'run', false; ...
'ecRxns', []; ...
'save', true; ...
'saveAs', 'dlkcat'; ...
'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 ~p.run && ~isempty(p.ecRxns)
error('kcatsFromDLKcat:ecRxns', ['ecRxns selects the reactions DLKcat is run for, so ' ...
'it needs ''run'', true.'])
end
if p.save
file = fullfile(modelAdapter.params.path, 'data', 'kcats', [char(p.saveAs) '.tsv']);
if p.run
warnSubset(p.ecRxns, numel(model.ec.rxns), file);
end
end

runFile = fullfile(modelAdapter.params.path, 'data', 'DLKcat_run.txt');
if p.run
% As runDLKcat: on macOS, MATLAB started from the launcher may not
% have Docker on its path.
if ismac
setenv('PATH', strcat('/usr/local/bin', ':', getenv("PATH")));
end
[status, ~] = system('docker --version');
if status ~= 0
error('kcatsFromDLKcat:noDocker', ['DLKcat runs in a Docker image, but Docker was ' ...
'not found. Install it from https://www.docker.com, or call kcatsFromDLKcat ' ...
'without ''run'' to read an earlier DLKcat output.'])
end
writeDLKcatInput(model, 'ecRxns', p.ecRxns, 'modelAdapter', modelAdapter, 'overwrite', true);
runDLKcat('modelAdapter', modelAdapter);
fid = fopen(runFile, 'w');
if ~isempty(p.ecRxns)
fprintf(fid, 'ecRxns: %d\n', countRxns(p.ecRxns));
end
fclose(fid);
else
outFile = fullfile(modelAdapter.params.path, 'data', 'DLKcat.tsv');
if ~isfile(outFile)
error('kcatsFromDLKcat:noOutput', ['No DLKcat output at %s. Run ' ...
'kcatsFromDLKcat(model, ''run'', true) first.'], outFile)
end
end
kcatList = readDLKcatOutput(model, 'modelAdapter', modelAdapter);

if p.save
parameters = struct('run', logical(p.run));
if isfile(runFile)
tok = regexp(fileread(runFile), 'ecRxns:\s*(\d+)', 'tokens', 'once');
if ~isempty(tok)
parameters.ec_rxns = str2double(tok{1});
end
end
saveKcatSource(model, kcatList, p.saveAs, 'modelAdapter', modelAdapter, ...
'generator', 'kcatsFromDLKcat', 'parameters', parameters);
fprintf('DLKcat kcats for %d reactions written to %s\n', ...
numel(unique(kcatList.rxns(kcatList.kcats > 0))), 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
Loading
Loading