From 9f7a12f0bd7fc3bb79de31e9caaa2a3eae015dce Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Fri, 16 Oct 2026 12:43:44 +0200 Subject: [PATCH 1/3] Keep gathered kcats as kcat source files, one per source Gathered kcats are written to data/kcats/.tsv in the project, keyed by conventional reaction and enzyme set rather than model.ec.rxns, so a file stays valid when the ecModel is rebuilt and applies to full and light models alike. geckopy reads and writes the same files. - kcatsFromBrenda, kcatsFromDLKcat and kcatsFromOKP gather one source each, with that source's own options, and write its file. OKP keeps a job and a result per method, and only uses a stored result of the current job. - saveKcatSource writes any kcat list; readKcatSource reads a file; loadKcatSource maps a file to a model and reports rows with an invalid kcat, rows that fit no reaction, and rows that match only part of an enzyme set. - getEcRxnKeys gives the conventional reaction and enzyme set of every model.ec reaction. - fetchOpenKineticsPredictor takes the file to store the result in. Original-commit: https://github.com/SysBioChalmers/GECKO/commit/f3a696cade6ea3129f036723954cff859212dbc1 --- .gitattributes | 3 + src/gather_kcats/fetchOpenKineticsPredictor.m | 11 +- src/gather_kcats/getEcRxnKeys.m | 51 ++++ src/gather_kcats/kcatsFromBrenda.m | 97 ++++++ src/gather_kcats/kcatsFromDLKcat.m | 114 +++++++ src/gather_kcats/kcatsFromOKP.m | 214 ++++++++++++++ src/gather_kcats/loadKcatSource.m | 218 ++++++++++++++ src/gather_kcats/private/checkNameValue.m | 21 ++ src/gather_kcats/private/countRxns.m | 8 + src/gather_kcats/private/formatKcatNumber.m | 24 ++ src/gather_kcats/private/kcatSourceFile.m | 58 ++++ src/gather_kcats/private/kcatSourceJson.m | 61 ++++ src/gather_kcats/private/parseHeaderLine.m | 14 + src/gather_kcats/private/parseKcatNumber.m | 11 + src/gather_kcats/private/unquoteField.m | 11 + src/gather_kcats/private/warnSubset.m | 7 + src/gather_kcats/readKcatSource.m | 137 +++++++++ src/gather_kcats/saveKcatSource.m | 226 ++++++++++++++ test/VerificationMatrix.tsv | 7 + .../ecTestGEM/data/kcatSourceFixture.tsv | 16 + test/unit_tests/kcatSourceTests.m | 279 ++++++++++++++++++ 21 files changed, 1586 insertions(+), 2 deletions(-) create mode 100644 .gitattributes create mode 100644 src/gather_kcats/getEcRxnKeys.m create mode 100644 src/gather_kcats/kcatsFromBrenda.m create mode 100644 src/gather_kcats/kcatsFromDLKcat.m create mode 100644 src/gather_kcats/kcatsFromOKP.m create mode 100644 src/gather_kcats/loadKcatSource.m create mode 100644 src/gather_kcats/private/checkNameValue.m create mode 100644 src/gather_kcats/private/countRxns.m create mode 100644 src/gather_kcats/private/formatKcatNumber.m create mode 100644 src/gather_kcats/private/kcatSourceFile.m create mode 100644 src/gather_kcats/private/kcatSourceJson.m create mode 100644 src/gather_kcats/private/parseHeaderLine.m create mode 100644 src/gather_kcats/private/parseKcatNumber.m create mode 100644 src/gather_kcats/private/unquoteField.m create mode 100644 src/gather_kcats/private/warnSubset.m create mode 100644 src/gather_kcats/readKcatSource.m create mode 100644 src/gather_kcats/saveKcatSource.m create mode 100644 test/unit_tests/ecTestGEM/data/kcatSourceFixture.tsv create mode 100644 test/unit_tests/kcatSourceTests.m diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 00000000..814b1d13 --- /dev/null +++ b/.gitattributes @@ -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 diff --git a/src/gather_kcats/fetchOpenKineticsPredictor.m b/src/gather_kcats/fetchOpenKineticsPredictor.m index a992f856..1666e498 100644 --- a/src/gather_kcats/fetchOpenKineticsPredictor.m +++ b/src/gather_kcats/fetchOpenKineticsPredictor.m @@ -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). % @@ -56,7 +59,8 @@ 'jobId', []; ... 'wait', []; ... 'pollInterval', []; ... - 'modelAdapter', []}); + 'modelAdapter', []; ... + 'outFile', []}); useStored = p.useStored; jobId = p.jobId; wait = p.wait; @@ -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 diff --git a/src/gather_kcats/getEcRxnKeys.m b/src/gather_kcats/getEcRxnKeys.m new file mode 100644 index 00000000..652750ea --- /dev/null +++ b/src/gather_kcats/getEcRxnKeys.m @@ -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_ (full model) or a leading _ (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 diff --git a/src/gather_kcats/kcatsFromBrenda.m b/src/gather_kcats/kcatsFromBrenda.m new file mode 100644 index 00000000..39721861 --- /dev/null +++ b/src/gather_kcats/kcatsFromBrenda.m @@ -0,0 +1,97 @@ +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 + +checkNameValue(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 = kcatSourceFile(p.saveAs, [], modelAdapter); % fails before the matching, not after + 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 diff --git a/src/gather_kcats/kcatsFromDLKcat.m b/src/gather_kcats/kcatsFromDLKcat.m new file mode 100644 index 00000000..fb50a39e --- /dev/null +++ b/src/gather_kcats/kcatsFromDLKcat.m @@ -0,0 +1,114 @@ +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 + +checkNameValue(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 = kcatSourceFile(p.saveAs, [], modelAdapter); + 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 diff --git a/src/gather_kcats/kcatsFromOKP.m b/src/gather_kcats/kcatsFromOKP.m new file mode 100644 index 00000000..4e0ed5d1 --- /dev/null +++ b/src/gather_kcats/kcatsFromOKP.m @@ -0,0 +1,214 @@ +function [kcatList, done] = kcatsFromOKP(model, varargin) +% kcatsFromOKP Gather kcats from OpenKineticsPredictor and keep them as a source file. +% +% Predicts kcats with one OpenKineticsPredictor (OKP) method and writes +% them to data/kcats/okp_.tsv in the project folder (see +% saveKcatSource). OKP is a web service with a job queue, so this function +% can be called repeatedly: +% +% - without a job for this method, it submits one +% (submitOpenKineticsPredictor); +% - with a job, it fetches the result (fetchOpenKineticsPredictor). When +% the job is done the file is written; otherwise it reports the job's +% status and writes nothing. +% +% 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; useStored only uses a result of the current job. An API key +% is needed (see resolveOkpApiKey). +% +% Parameters +% ---------- +% model : struct +% an ecModel in GECKO 3 format (with ecModel.ec structure), with +% metabolite SMILES (findMetSmiles) and protein sequences. +% +% Name-Value Arguments +% -------------------- +% method : char +% OKP method, e.g. 'CataPro', 'CatPred', 'CatRange', 'DLKcat', +% 'EITLEM', 'KinForm-H', 'KinForm-L', 'OmniESI', 'TurNuP', 'UniKP' +% (default: params.okp.method of the model adapter, or 'CataPro'). +% wait : logical +% wait until the job is done (default false). +% resubmit : logical +% submit a new job even if one exists for this method; its earlier +% result is removed (default false). +% useStored : logical +% read the stored result of this method's job, without contacting OKP +% (default false). +% ecRxns : logical +% for which model.ec.rxns kcats are predicted (default: all). Only +% used when a job is submitted. +% save : logical +% write the source file (default true). +% saveAs : char +% name of the source file, without .tsv (default 'okp_', in +% lower case, with other characters than letters and digits as '_'). +% modelAdapter : ModelAdapter +% a loaded model adapter (default: the current default model adapter). +% +% Returns +% ------- +% kcatList : struct +% the kcat values, as returned by readOpenKineticsPredictorOutput; +% empty while the job is not done. +% done : logical +% whether the result was available. +% +% Examples +% -------- +% kcatsFromOKP(ecModel, 'method', 'CataPro'); % submits +% kcatsFromOKP(ecModel, 'method', 'CataPro'); % later: fetches and writes +% kcatsFromOKP(ecModel, 'method', 'UniKP', 'wait', true); +% +% See also +% -------- +% submitOpenKineticsPredictor, fetchOpenKineticsPredictor, kcatsFromBrenda, saveKcatSource + +checkNameValue(varargin, {'method','wait','resubmit','useStored','ecRxns','save','saveAs','modelAdapter'}, 'kcatsFromOKP'); +p = parseGECKOargs(varargin, { ... + 'method', []; ... + 'wait', false; ... + 'resubmit', false; ... + 'useStored', false; ... + 'ecRxns', []; ... + 'save', true; ... + 'saveAs', []; ... + '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 +params = modelAdapter.params; +method = char(p.method); +if isempty(method) + if isfield(params, 'okp') && isfield(params.okp, 'method') && ~isempty(params.okp.method) + method = char(params.okp.method); + else + method = 'CataPro'; + end +end +tag = lower(regexprep(method, '[^A-Za-z0-9]+', '_')); +dataDir = fullfile(params.path, 'data'); +jobFile = fullfile(dataDir, ['OKP_job_' tag '.txt']); +outFile = fullfile(dataDir, ['OKP_output_' tag '.csv']); +name = char(p.saveAs); +if isempty(name) + name = ['okp_' tag]; +end +job = readJobFile(jobFile); +submitting = ~p.useStored && (p.resubmit || ~isfield(job, 'jobId')); +if ~isempty(p.ecRxns) && ~submitting + error('kcatsFromOKP:ecRxns', ['ecRxns selects the reactions of a new job, but this call ' ... + 'submits none (a job for %s exists; use ''resubmit'', true for a new one).'], method) +end +if p.save + file = kcatSourceFile(name, [], modelAdapter); % fails before submitting, not after + if submitting + warnSubset(p.ecRxns, numel(model.ec.rxns), file); + end +end + +if p.useStored + if ~isfield(job, 'jobId') + adoptSharedResult(model, modelAdapter, dataDir, method, jobFile, outFile); + job = readJobFile(jobFile); + end + if ~isfile(outFile) || ~isfield(job, 'resultOf') || ~strcmp(job.resultOf, job.jobId) + error('kcatsFromOKP:noStored', ['No stored result of the current %s job (%s). Call ' ... + 'kcatsFromOKP without ''useStored'' to fetch it.'], method, jobFile) + end + kcatList = readOpenKineticsPredictorOutput(model, 'outFile', outFile, 'modelAdapter', modelAdapter); + done = true; +else + if submitting + if isfile(outFile) + delete(outFile); + end + jobId = submitOpenKineticsPredictor(model, 'method', method, 'ecRxns', p.ecRxns, ... + 'overwrite', true, 'modelAdapter', modelAdapter); + job = struct('jobId', jobId, 'method', method, 'submittedAt', ... + char(datetime('now', 'TimeZone', 'UTC', 'Format', 'yyyy-MM-dd''T''HH:mm:ss''Z'''))); + if ~isempty(p.ecRxns) + job.ecRxns = sprintf('%d', countRxns(p.ecRxns)); + end + writeJobFile(jobFile, job); + end + [done, kcatList] = fetchOpenKineticsPredictor(model, 'jobId', job.jobId, 'wait', p.wait, ... + 'outFile', outFile, 'modelAdapter', modelAdapter); + if done + job.resultOf = job.jobId; + writeJobFile(jobFile, job); + end +end + +if ~done + fprintf('Call kcatsFromOKP(model, ''method'', ''%s'') again later to collect the result.\n', method); + return +end +if p.save + parameters = struct('method', method); + if isfield(job, 'ecRxns') + parameters.ec_rxns = str2double(job.ecRxns); + end + saveKcatSource(model, kcatList, name, 'modelAdapter', modelAdapter, ... + 'generator', 'kcatsFromOKP', 'parameters', parameters, 'data', ['OKP job ' job.jobId]); + fprintf('OKP (%s) kcats for %d reactions written to %s\n', method, ... + numel(unique(kcatList.rxns(kcatList.kcats > 0))), file); +end +end + +function adoptSharedResult(model, modelAdapter, dataDir, method, jobFile, outFile) +% A result that fetchOpenKineticsPredictor stored by itself, in +% data/OKP_output.csv with data/OKP_job.txt, is taken over for this method +% only when its job is for this method and its predictions are labelled +% with this method. +shared = fullfile(dataDir, 'OKP_output.csv'); +sharedJob = readJobFile(fullfile(dataDir, 'OKP_job.txt')); +ok = isfile(shared) && isfield(sharedJob, 'jobId') && isfield(sharedJob, 'method') && ... + strcmpi(sharedJob.method, method); +if ok + list = readOpenKineticsPredictorOutput(model, 'outFile', shared, 'modelAdapter', modelAdapter); + predicted = list.kcatSource(startsWith(list.kcatSource, 'OKP-') & ... + ~ismember(lower(list.kcatSource), {'okp-brenda','okp-sabio-rk','okp-uniprot'})); + ok = ~isempty(predicted) && all(strcmpi(predicted, ['OKP-' method])); +end +if ~ok + error('kcatsFromOKP:noStored', ['No stored OKP result for %s (%s). Call kcatsFromOKP ' ... + 'without ''useStored'' to submit or fetch a job.'], method, outFile) +end +copyfile(shared, outFile); +sharedJob.method = method; +sharedJob.resultOf = sharedJob.jobId; +writeJobFile(jobFile, sharedJob); +end + +function job = readJobFile(file) +% Fields of an OKP job file ('key: value' per line), as a struct. +job = struct(); +if ~isfile(file) + return +end +lines = strsplit(fileread(file), newline); +for i = 1:numel(lines) + tok = regexp(lines{i}, '^([A-Za-z]\w*):\s*(.*?)\s*$', 'tokens', 'once'); + if ~isempty(tok) + job.(tok{1}) = tok{2}; + end +end +end + +function writeJobFile(file, job) +fid = fopen(file, 'w'); +for f = {'jobId','method','submittedAt','ecRxns','resultOf'} + if isfield(job, f{1}) + fprintf(fid, '%s: %s\n', f{1}, job.(f{1})); + end +end +fclose(fid); +end diff --git a/src/gather_kcats/loadKcatSource.m b/src/gather_kcats/loadKcatSource.m new file mode 100644 index 00000000..87d99867 --- /dev/null +++ b/src/gather_kcats/loadKcatSource.m @@ -0,0 +1,218 @@ +function [kcatList, report] = loadKcatSource(model, name, varargin) +% loadKcatSource Read a kcat source file and map its rows to an ecModel. +% +% Reads data/kcats/.tsv (see saveKcatSource) and maps every row to +% the model.ec reactions it applies to, by conventional reaction and +% enzyme set (see getEcRxnKeys): +% +% 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 model.ec row with that +% enzyme set; a row with only rxn to every model.ec row of that +% conventional reaction (curated rows; saveKcatSource 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 : struct +% an ecModel in GECKO 3 format (with ecModel.ec structure). +% name : char +% name of the source (data/kcats/.tsv in the project folder), or +% the path to a kcat source file. +% +% Name-Value Arguments +% -------------------- +% modelAdapter : ModelAdapter +% a loaded model adapter (default: the current default model adapter). +% +% Returns +% ------- +% kcatList : struct +% kcat values per model.ec reaction, as used by mergeKcats and +% assignKcatValues: source (the source name), rxns, kcats, +% substrates, eccodes, wildcardLvl and origin (from the detail +% column), kcatSource (the per-row source), log10sd, and line (the +% line of the file each entry comes from). +% report : struct +% - invalid: table of the file rows with a missing or invalid kcat; +% - unmatched: table of the file rows that apply to no model.ec row; +% - partial: table of rows that only partly match an enzyme set, with +% the model.ec reaction and its enzymes. +% +% Examples +% -------- +% kcatList = loadKcatSource(ecModel, 'brenda'); +% [kcatList, report] = loadKcatSource(ecModel, 'custom'); +% +% See also +% -------- +% saveKcatSource, readKcatSource, kcatsFromBrenda, getEcRxnKeys + +checkNameValue(varargin, {'modelAdapter'}, 'loadKcatSource'); +p = parseGECKOargs(varargin, {'modelAdapter', []}); + +if isfile(name) + filename = name; +else + 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 + filename = fullfile(modelAdapter.params.path, 'data', 'kcats', [name '.tsv']); + if ~isfile(filename) + error('loadKcatSource:missing', ['No kcat source file %s. Make it with the ' ... + 'kcatsFrom function of this source (kcatsFromBrenda, kcatsFromDLKcat, ' ... + 'kcatsFromOKP), or with saveKcatSource.'], filename) + end +end +[records, meta] = readKcatSource(filename); +sourceName = ''; +if isfield(meta, 'source') + sourceName = meta.source; +end + +[baseRxns, enzymeSets] = getEcRxnKeys(model); +nEc = numel(baseRxns); +nRec = numel(records.kcat); + +% Lookup tables: key -> model.ec rows. +byPair = groupIndex(strcat(baseRxns, '|', enzymeSets)); +byBase = groupIndex(baseRxns); +byEnz = groupIndex(enzymeSets); +% accession -> model.ec rows containing it, for partial matches. +accRows = cell(nEc, 1); +for i = 1:nEc + if ~isempty(enzymeSets{i}) + accRows{i} = strsplit(enzymeSets{i}, ' + '); + end +end +accIdx = repelem((1:nEc)', cellfun(@numel, accRows)); +byAcc = groupIndex([accRows{:}]', accIdx); + +invalid = ~(records.kcat > 0) | ~isfinite(records.kcat); +matchRec = cell(nRec, 1); +matchEc = cell(nRec, 1); +unmatched = false(nRec, 1); +partialRec = zeros(0, 1); partialEc = zeros(0, 1); +for r = find(~invalid)' + rxn = records.rxn{r}; + enz = normalizeEnzymes(records.enzymes{r}); + if ~isempty(rxn) && ~isempty(enz) + hit = lookup(byPair, [rxn '|' enz]); + elseif ~isempty(enz) + hit = lookup(byEnz, enz); + elseif ~isempty(rxn) + hit = lookup(byBase, rxn); + else + hit = []; + end + if isempty(hit) + unmatched(r) = true; + if ~isempty(enz) + wanted = strsplit(enz, ' + '); + candidates = unique(cell2mat(cellfun(@(a) lookup(byAcc, a), wanted(:), 'UniformOutput', false))); + if ~isempty(rxn) + candidates = intersect(candidates, lookup(byBase, rxn)); + end + for c = candidates(:)' + if mean(ismember(wanted, accRows{c})) >= 0.5 + partialRec(end+1, 1) = r; %#ok + partialEc(end+1, 1) = c; %#ok + end + end + end + else + matchRec{r} = repmat(r, numel(hit), 1); + matchEc{r} = hit(:); + end +end +rec = vertcat(matchRec{:}, zeros(0, 1)); +ec = vertcat(matchEc{:}, zeros(0, 1)); + +kcatList.source = sourceName; +kcatList.rxns = reshape(model.ec.rxns(ec), [], 1); +kcatList.kcats = records.kcat(rec); +kcatList.substrates = records.substrate(rec); +kcatList.eccodes = records.ec_code(rec); +kcatList.kcatSource = records.source(rec); +kcatList.log10sd = records.log10_sd(rec); +[kcatList.wildcardLvl, kcatList.origin] = parseDetail(records.detail(rec)); +kcatList.line = records.line(rec); + +report.invalid = table(records.line(invalid), records.rxn(invalid), records.enzymes(invalid), ... + records.kcat(invalid), 'VariableNames', {'line','rxn','enzymes','kcat'}); +report.unmatched = table(records.line(unmatched), records.rxn(unmatched), records.enzymes(unmatched), ... + records.kcat(unmatched), 'VariableNames', {'line','rxn','enzymes','kcat'}); +report.partial = table(records.line(partialRec), records.rxn(partialRec), records.enzymes(partialRec), ... + reshape(model.ec.rxns(partialEc), [], 1), enzymeSets(partialEc), ... + 'VariableNames', {'line','rxn','enzymes','ecRxn','ecEnzymes'}); +if any(invalid) + printOrange(sprintf('%d of %d rows of %s have no valid kcat; see report.invalid.\n', ... + sum(invalid), nRec, filename)); +end +if any(unmatched) + printOrange(sprintf('%d of %d rows of %s apply to no reaction of this model; see report.unmatched.\n', ... + sum(unmatched), nRec, filename)); +end +end + +function map = groupIndex(keys, rows) +% containers.Map from each key to the rows (default 1:n) that have it. +if nargin < 2 + rows = (1:numel(keys))'; +end +map = containers.Map('KeyType', 'char', 'ValueType', 'any'); +if isempty(keys) + return +end +[u, ~, g] = unique(keys(:)); +groups = accumarray(g, rows(:), [], @(x) {sort(x)}); +for i = 1:numel(u) + if ~isempty(u{i}) + map(u{i}) = groups{i}; + end +end +end + +function rows = lookup(map, key) +if isKey(map, key) + rows = map(key); +else + rows = zeros(0, 1); +end +end + +function enz = normalizeEnzymes(enz) +% Sorted and joined as getEcRxnKeys does, so hand-written rows match too. +if isempty(enz) + return +end +parts = strtrim(strsplit(enz, '+')); +parts = unique(parts(~cellfun(@isempty, parts))); +enz = strjoin(sort(parts(:)'), ' + '); +end + +function [wc, origin] = parseDetail(detail) +wc = nan(numel(detail), 1); +origin = nan(numel(detail), 1); +for i = 1:numel(detail) + t = regexp(detail{i}, 'wc=(\d+)', 'tokens', 'once'); + if ~isempty(t) + wc(i) = str2double(t{1}); + end + t = regexp(detail{i}, 'origin=(\d+)', 'tokens', 'once'); + if ~isempty(t) + origin(i) = str2double(t{1}); + end +end +end diff --git a/src/gather_kcats/private/checkNameValue.m b/src/gather_kcats/private/checkNameValue.m new file mode 100644 index 00000000..ed2e7bb0 --- /dev/null +++ b/src/gather_kcats/private/checkNameValue.m @@ -0,0 +1,21 @@ +function checkNameValue(args, names, functionName) +% checkNameValue Accept only name-value pairs with known names. +% +% The kcat source functions have no positional form, so a misspelled +% option name must be an error rather than being read as a positional +% argument by parseGECKOargs. +id = [functionName ':arguments']; +if mod(numel(args), 2) ~= 0 + error(id, 'Options of %s are given as name-value pairs.', functionName) +end +for k = 1:2:numel(args) + name = args{k}; + if ~(ischar(name) || (isstring(name) && isscalar(name))) + error(id, 'Option names of %s are text, not %s. Options: %s.', functionName, ... + class(name), strjoin(names, ', ')) + elseif ~any(strcmp(char(name), names)) + error(id, 'Unknown option of %s: %s. Options: %s.', functionName, char(name), ... + strjoin(names, ', ')) + end +end +end diff --git a/src/gather_kcats/private/countRxns.m b/src/gather_kcats/private/countRxns.m new file mode 100644 index 00000000..fd812fab --- /dev/null +++ b/src/gather_kcats/private/countRxns.m @@ -0,0 +1,8 @@ +function n = countRxns(ecRxns) +% countRxns Number of reactions an ecRxns argument selects (logical or indices). +if islogical(ecRxns) + n = sum(ecRxns); +else + n = numel(unique(ecRxns)); +end +end diff --git a/src/gather_kcats/private/formatKcatNumber.m b/src/gather_kcats/private/formatKcatNumber.m new file mode 100644 index 00000000..f7251a3e --- /dev/null +++ b/src/gather_kcats/private/formatKcatNumber.m @@ -0,0 +1,24 @@ +function t = formatKcatNumber(x) +% formatKcatNumber Number text of kcat source files. +% +% Whole numbers as integers, others in the shortest decimal form that reads +% back as the same double; 'Inf' or '-Inf' for infinity, empty for NaN. +% geckopy's format_number writes the same text. +if isnan(x) + t = ''; +elseif isinf(x) + t = 'Inf'; + if x < 0 + t = '-Inf'; + end +elseif x == round(x) && abs(x) < 1e15 + t = sprintf('%d', x); +else + for digits = 1:17 + t = sprintf('%.*g', digits, x); + if str2double(t) == x + break + end + end +end +end diff --git a/src/gather_kcats/private/kcatSourceFile.m b/src/gather_kcats/private/kcatSourceFile.m new file mode 100644 index 00000000..3d112cf3 --- /dev/null +++ b/src/gather_kcats/private/kcatSourceFile.m @@ -0,0 +1,58 @@ +function filename = kcatSourceFile(name, folder, modelAdapter) +% kcatSourceFile Path of a kcat source file that may be written. +% +% Used by saveKcatSource and the kcatsFrom functions before any +% work is done: checks the name, resolves the folder (data/kcats in the +% project folder by default), and refuses a file that is curated by hand: +% named custom.tsv, or whose header says 'kind: curated' or 'source: +% custom' (in any case, also when a spreadsheet quoted the line). + +name = char(name); +if isempty(regexp(name, '^[A-Za-z0-9_.-]+$', 'once')) + error('saveKcatSource:invalidName', ... + 'The source name "%s" may only contain letters, digits, ''_'', ''-'' and ''.''.', name) +end +if isempty(folder) + 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 + folder = fullfile(modelAdapter.params.path, 'data', 'kcats'); +end +filename = fullfile(char(folder), [name '.tsv']); +if strcmpi(name, 'custom') + error('saveKcatSource:curated', ['custom.tsv holds curated kcats and is not ' ... + 'written by GECKO; edit it by hand.']) +end +if isfile(filename) + meta = readHeaderLines(filename); + if (isfield(meta, 'kind') && strcmpi(meta.kind, 'curated')) || ... + (isfield(meta, 'source') && strcmpi(meta.source, 'custom')) + error('saveKcatSource:curated', '%s is curated by hand and is not overwritten.', filename) + end +end +end + +function meta = readHeaderLines(filename) +% Only the header lines at the top, so that a damaged table does not stop +% a gathered file from being written again. +meta = struct(); +fid = fopen(filename, 'r', 'n', 'UTF-8'); +if fid < 0 + return +end +cleanup = onCleanup(@() fclose(fid)); +line = fgetl(fid); +if ischar(line) + line = regexprep(line, ['^' char(65279)], ''); +end +while ischar(line) && (isempty(strtrim(line)) || startsWith(unquoteField(line), '#')) + [key, value] = parseHeaderLine(line); + if ~isempty(key) + meta.(key) = value; + end + line = fgetl(fid); +end +end diff --git a/src/gather_kcats/private/kcatSourceJson.m b/src/gather_kcats/private/kcatSourceJson.m new file mode 100644 index 00000000..86808bd1 --- /dev/null +++ b/src/gather_kcats/private/kcatSourceJson.m @@ -0,0 +1,61 @@ +function j = kcatSourceJson(v) +% kcatSourceJson The parameters line of kcat source files as JSON. +% +% Compact JSON with keys sorted at every level. Numbers are written with +% formatKcatNumber, and NaN and infinity as null; logical values as true +% or false; char as a string (non-ASCII characters as they are); a cell +% array as a list; a numeric or logical vector with more than one element +% as a list. geckopy's save_kcat_source writes the same text. +if isstruct(v) && isscalar(v) + f = sort(fieldnames(v)); + parts = cell(1, numel(f)); + for k = 1:numel(f) + parts{k} = [jsonString(f{k}) ':' kcatSourceJson(v.(f{k}))]; + end + j = ['{' strjoin(parts, ',') '}']; +elseif ischar(v) || (isstring(v) && isscalar(v)) + j = jsonString(char(v)); +elseif iscell(v) + j = ['[' strjoin(cellfun(@kcatSourceJson, v(:)', 'UniformOutput', false), ',') ']']; +elseif (isnumeric(v) || islogical(v)) && numel(v) ~= 1 + j = ['[' strjoin(arrayfun(@kcatSourceJson, v(:)', 'UniformOutput', false), ',') ']']; +elseif islogical(v) + j = 'false'; + if v + j = 'true'; + end +elseif isnumeric(v) + if ~isfinite(v) + j = 'null'; + else + j = formatKcatNumber(double(v)); + end +elseif isempty(v) + j = 'null'; +else + error('kcatSourceJson:type', 'Parameters of class %s cannot be written.', class(v)) +end +end + +function s = jsonString(s) +% JSON escapes, as Python's json.dumps(ensure_ascii=False). +out = ''; +for c = s + switch c + case '"', out = [out '\"']; %#ok + case '\', out = [out '\\']; %#ok + case char(8), out = [out '\b']; %#ok + case char(9), out = [out '\t']; %#ok + case char(10), out = [out '\n']; %#ok + case char(12), out = [out '\f']; %#ok + case char(13), out = [out '\r']; %#ok + otherwise + if double(c) < 32 + out = [out sprintf('\\u%04x', double(c))]; %#ok + else + out = [out c]; %#ok + end + end +end +s = ['"' out '"']; +end diff --git a/src/gather_kcats/private/parseHeaderLine.m b/src/gather_kcats/private/parseHeaderLine.m new file mode 100644 index 00000000..ee17bbc4 --- /dev/null +++ b/src/gather_kcats/private/parseHeaderLine.m @@ -0,0 +1,14 @@ +function [key, value] = parseHeaderLine(line) +% parseHeaderLine Key and value of a '# key: value' header line. +% +% The line may be wrapped in quotes by a spreadsheet. The key is returned +% in lower case; it must start with a letter and be at most 63 letters, +% digits or underscores long. Other lines give an empty key. geckopy reads +% header lines the same way. +key = ''; value = ''; +tok = regexp(unquoteField(line), '^\s*#\s*([A-Za-z][A-Za-z0-9_]{0,62})\s*:\s?(.*)$', 'tokens', 'once'); +if ~isempty(tok) + key = lower(tok{1}); + value = strtrim(tok{2}); +end +end diff --git a/src/gather_kcats/private/parseKcatNumber.m b/src/gather_kcats/private/parseKcatNumber.m new file mode 100644 index 00000000..f17b49a0 --- /dev/null +++ b/src/gather_kcats/private/parseKcatNumber.m @@ -0,0 +1,11 @@ +function x = parseKcatNumber(texts) +% parseKcatNumber Read numbers written in kcat source files. +% +% Accepts plain decimal numbers with an optional exponent (12, 0.5, 1e-05, +% 2.5E+3) and Inf; anything else, such as '5,5' from a spreadsheet with a +% decimal comma, becomes NaN instead of being misread. geckopy's +% read_kcat_source accepts the same text. +x = nan(numel(texts), 1); +ok = ~cellfun(@isempty, regexp(texts(:), '^[-+]?(\d+\.?\d*|\.\d+)([eE][-+]?\d+)?$|^[-+]?Inf$', 'once')); +x(ok) = str2double(texts(ok)); +end diff --git a/src/gather_kcats/private/unquoteField.m b/src/gather_kcats/private/unquoteField.m new file mode 100644 index 00000000..36c95aa1 --- /dev/null +++ b/src/gather_kcats/private/unquoteField.m @@ -0,0 +1,11 @@ +function t = unquoteField(t) +% unquoteField 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; spaces inside the quotes are +% trimmed too. +t = strtrim(t); +if numel(t) >= 2 && t(1) == '"' && t(end) == '"' + t = strtrim(strrep(t(2:end-1), '""', '"')); +end +end diff --git a/src/gather_kcats/private/warnSubset.m b/src/gather_kcats/private/warnSubset.m new file mode 100644 index 00000000..00064a5b --- /dev/null +++ b/src/gather_kcats/private/warnSubset.m @@ -0,0 +1,7 @@ +function warnSubset(ecRxns, nRxns, file) +% warnSubset Say that gathering for some reactions replaces the whole 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 diff --git a/src/gather_kcats/readKcatSource.m b/src/gather_kcats/readKcatSource.m new file mode 100644 index 00000000..063ab8b0 --- /dev/null +++ b/src/gather_kcats/readKcatSource.m @@ -0,0 +1,137 @@ +function [records, meta] = readKcatSource(filename) +% readKcatSource Read a kcat source file as it is, without a model. +% +% Reads a file written by saveKcatSource or a kcatsFrom function, +% or curated by hand in the same format. Use loadKcatSource to map the +% rows to an ecModel. +% +% The file is read as UTF-8 (a leading byte-order mark is ignored). The +% '#' lines before the column header are the header (keys in lower case); +% '#' 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 loadKcatSource reports. A row with more filled +% fields than there are columns is an error. +% +% Parameters +% ---------- +% filename : char +% path to the kcat source file. +% +% Returns +% ------- +% records : struct +% one field per column, each with one entry per data row: kcat and +% log10_sd as double (NaN when empty or not a number), every other +% column as a cell of char. Columns other than the standard ones are +% kept. records.line is the line number of each row in the file. +% meta : struct +% the '#' header lines, as field: value (for instance meta.source, +% meta.kind). +% +% Examples +% -------- +% [records, meta] = readKcatSource(fullfile(projectPath,'data','kcats','brenda.tsv')); +% +% See also +% -------- +% loadKcatSource, saveKcatSource + +text = readUtf8(filename); +lines = regexp(text, '\r?\n', 'split'); + +meta = struct(); +k = 1; +while k <= numel(lines) && (isempty(strtrim(lines{k})) || startsWith(unquoteField(lines{k}), '#')) + [key, value] = parseHeaderLine(lines{k}); + if ~isempty(key) + meta.(key) = value; + end + k = k + 1; +end +if k > numel(lines) + error('readKcatSource:noHeader', '%s has no column header.', filename) +end +columns = cellfun(@unquoteField, strsplit(lines{k}, sprintf('\t'), ... + 'CollapseDelimiters', false), 'UniformOutput', false); +if ~all(ismember({'kcat','source'}, columns)) || ~any(ismember({'rxn','enzymes'}, columns)) + error('readKcatSource:columns', ['%s needs the columns kcat and source, and rxn ' ... + 'and/or enzymes.'], filename) +end +named = columns(~cellfun(@isempty, columns)); +if numel(unique(named)) < numel(named) || any(strcmp(columns, 'line')) + error('readKcatSource:columns', '%s has a column name twice, or a column named line.', filename) +end +nCol = numel(columns); + +trimmed = strtrim(lines); +skip = cellfun(@isempty, trimmed) | startsWith(trimmed, '#') | startsWith(trimmed, '"#'); +skip(1:k) = true; +kept = find(~skip); +L = string(lines(kept)); +L = L(:); +F = strings(numel(kept), nCol); +nTabs = count(L, sprintf('\t')); +regular = nTabs == nCol - 1; +if any(regular) + % Rows with exactly one field per column, split in one go. + F(regular, :) = reshape(split(L(regular), sprintf('\t'), 2), [], nCol); +end +for r = find(~regular)' + fields = strsplit(char(L(r)), sprintf('\t'), 'CollapseDelimiters', false); + if numel(fields) > nCol + if any(~cellfun(@isempty, strtrim(fields(nCol+1:end)))) + error('readKcatSource:fields', '%s, line %d has more fields than the %d columns.', ... + filename, kept(r), nCol) + end + fields = fields(1:nCol); + end + fields(end+1:nCol) = {''}; + F(r, :) = string(fields); +end +F = strtrim(F); +quoted = find(startsWith(F, '"')); +for q = quoted(:)' + F(q) = unquoteField(char(F(q))); +end +filled = any(F ~= "", 2); % a row of tabs only is skipped +values = cellstr(F(filled, :)); +if isempty(values) + values = cell(0, nCol); +end +lineNo = kept(filled); +lineNo = lineNo(:); +n = numel(lineNo); + +records = struct(); +for c = 1:nCol + if isempty(columns{c}) + continue + end + name = matlab.lang.makeValidName(columns{c}); + v = values(:, c); + if any(strcmp(columns{c}, {'kcat','log10_sd'})) + v = parseKcatNumber(v); + end + records.(name) = v; +end +for c = {'rxn','enzymes','substrate','source','detail','ec_code','note'} + if ~isfield(records, c{1}) + records.(c{1}) = repmat({''}, n, 1); + end +end +if ~isfield(records, 'log10_sd') + records.log10_sd = nan(n, 1); +end +records.line = lineNo; +end + +function text = readUtf8(filename) +fid = fopen(filename, 'r', 'n', 'UTF-8'); +if fid < 0 + error('readKcatSource:read', 'Cannot read %s', filename) +end +cleanup = onCleanup(@() fclose(fid)); +text = fread(fid, [1 Inf], '*char'); +text = regexprep(text, ['^' char(65279)], ''); +end diff --git a/src/gather_kcats/saveKcatSource.m b/src/gather_kcats/saveKcatSource.m new file mode 100644 index 00000000..37101553 --- /dev/null +++ b/src/gather_kcats/saveKcatSource.m @@ -0,0 +1,226 @@ +function filename = saveKcatSource(model, kcatList, name, varargin) +% saveKcatSource Write a kcatList to a kcat source file of the project. +% +% Writes the kcat values of kcatList to data/kcats/.tsv in the +% project folder, the file that loadKcatSource reads back. Rows are keyed +% by their conventional reaction and enzyme set (see getEcRxnKeys), not +% by model.ec.rxns, so the file stays valid when the ecModel is rebuilt. +% kcatsFromBrenda, kcatsFromDLKcat and kcatsFromOKP write these files for +% their sources; use this function for a kcatList made in another way. +% +% The file is UTF-8 text. 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 all columns from rxn onwards. +% Left out are: +% - rows whose kcat is not a positive finite number (0 means no value +% was found); +% - rows of model.ec reactions without enzymes, which have no enzyme +% constraint and could not be told apart from their isozymes. +% +% Parameters +% ---------- +% model : struct +% the ecModel the kcatList refers to. +% kcatList : struct +% kcat values as returned by kcatsFromBrenda, fuzzyKcatMatching, +% readDLKcatOutput, readOpenKineticsPredictorOutput or mergeKcats: +% fields rxns and kcats, optionally source, substrates, eccodes, +% wildcardLvl, origin, kcatSource and log10sd. +% name : char +% name of the source; the file is .tsv. +% +% Name-Value Arguments +% -------------------- +% modelAdapter : ModelAdapter +% a loaded model adapter (default: the current default model adapter). +% folder : char +% folder to write to (default: data/kcats in the project folder). +% generator : char +% what produced the values (default: 'saveKcatSource'). +% parameters : struct +% settings used to produce the values, written as compact JSON with +% sorted keys (default: none). +% data : char +% version of the data the values come from, e.g. a BRENDA release +% (default: empty). +% kind : char +% 'gathered' (default) or 'curated'. +% +% A file whose header says 'kind: curated' or 'source: custom', or named +% custom.tsv, is never written. +% +% Returns +% ------- +% filename : char +% the file written. +% +% Examples +% -------- +% kcatList = fuzzyKcatMatching(ecModel); +% saveKcatSource(ecModel, kcatList, 'mybrenda'); +% +% See also +% -------- +% kcatsFromBrenda, loadKcatSource, readKcatSource, getEcRxnKeys + +names = {'modelAdapter','folder','generator','parameters','data','kind'}; +checkNameValue(varargin, names, 'saveKcatSource'); +p = parseGECKOargs(varargin, { ... + 'modelAdapter', []; ... + 'folder', []; ... + 'generator', 'saveKcatSource'; ... + 'parameters', struct(); ... + 'data', ''; ... + 'kind', 'gathered'}); +kind = validatestring(p.kind, {'gathered','curated'}); +name = char(name); +filename = kcatSourceFile(name, p.folder, p.modelAdapter); + +[baseRxns, enzymeSets] = getEcRxnKeys(model); +[found, idx] = ismember(kcatList.rxns(:), model.ec.rxns); +if ~all(found) + error('saveKcatSource:unknownRxns', '%s', ravenList( ... + 'These kcatList reactions are not in model.ec.rxns:', kcatList.rxns(~found), false)); +end + +n = numel(kcatList.rxns); +kcats = double(kcatList.kcats(:)); +if numel(kcats) ~= n + error('saveKcatSource:fields', 'kcatList.kcats and kcatList.rxns differ in length.') +end +substrate = textColumn(kcatList, 'substrates', n); +ecCode = textColumn(kcatList, 'eccodes', n); +if isfield(kcatList, 'kcatSource') + source = textColumn(kcatList, 'kcatSource', n); +elseif isfield(kcatList, 'source') + source = repmat({char(kcatList.source)}, n, 1); +else + source = repmat({name}, n, 1); +end +wc = numColumn(kcatList, 'wildcardLvl', n); +origin = numColumn(kcatList, 'origin', n); +sd = numColumn(kcatList, 'log10sd', n); + +keep = kcats > 0 & isfinite(kcats); +noEnzymes = keep & cellfun(@isempty, enzymeSets(idx)); +if any(noEnzymes) + fprintf('%d kcats of reactions without enzymes are not written to %s.\n', sum(noEnzymes), filename); +end +keep = find(keep & ~noEnzymes); + +clean = @(c) regexprep(string(c(:)), '[\t\r\n]', ' '); +rows = strings(numel(keep), 9); +if ~isempty(keep) + wcText = strings(numel(keep), 1); + hasWc = ~isnan(wc(keep)); + wcText(hasWc) = "wc=" + string(fix(wc(keep(hasWc)))); + originText = strings(numel(keep), 1); + hasOrigin = ~isnan(origin(keep)); + originText(hasOrigin) = "origin=" + string(fix(origin(keep(hasOrigin)))); + detail = wcText + "; " + originText; + detail(~hasWc) = originText(~hasWc); + detail(~hasOrigin) = wcText(~hasOrigin); + rows(:, 1) = clean(baseRxns(idx(keep))); + rows(:, 2) = clean(enzymeSets(idx(keep))); + rows(:, 3) = clean(substrate(keep)); + rows(:, 4) = formatNumbers(kcats(keep)); + rows(:, 5) = formatNumbers(sd(keep)); + rows(:, 6) = clean(source(keep)); + rows(:, 7) = detail; + rows(:, 8) = clean(ecCode(keep)); + % Sorted by the columns in this order, as geckopy sorts tuples: one key + % per row, the fields joined by a character below any in the text. + [~, order] = sort(join(rows(:, [1 2 3 6 4 5 7 8 9]), char(1), 2)); + rows = rows(order, :); +end + +folder = fileparts(filename); +if ~isfolder(folder) + mkdir(folder); +end +fid = fopen(filename, 'w', 'n', 'UTF-8'); +if fid < 0 + error('saveKcatSource:write', 'Cannot write %s', filename) +end +cleanup = onCleanup(@() fclose(fid)); +fprintf(fid, '# source: %s\n', name); +fprintf(fid, '# generator: %s (GECKO %s)\n', cleanText(p.generator), geckoVersion()); +fprintf(fid, '# parameters: %s\n', kcatSourceJson(p.parameters)); +fprintf(fid, '# data: %s\n', cleanText(p.data)); +fprintf(fid, '# kind: %s\n', kind); +fprintf(fid, '# created: %s\n', char(datetime('now', 'TimeZone', 'UTC', 'Format', 'yyyy-MM-dd''T''HH:mm:ss''Z'''))); +fprintf(fid, 'rxn\tenzymes\tsubstrate\tkcat\tlog10_sd\tsource\tdetail\tec_code\tnote\n'); +if ~isempty(rows) + fprintf(fid, '%s\n', join(rows, sprintf('\t'), 2)); +end +end + +function col = textColumn(s, field, n) +% A kcatList field as an n x 1 cell of char; a cell of several entries +% (e.g. all substrates of a reaction) is joined by '; '. +col = repmat({''}, n, 1); +if ~isfield(s, field) || isempty(s.(field)) + return +end +v = s.(field); +if ischar(v) || isstring(v) + v = cellstr(v); +end +if numel(v) ~= n + error('saveKcatSource:fields', 'kcatList.%s and kcatList.rxns differ in length.', field) +end +for i = 1:n + e = v{i}; + if iscell(e) + e = strjoin(cellfun(@char, e(:)', 'UniformOutput', false), '; '); + elseif isstring(e) + e = char(e); + elseif ~ischar(e) + e = ''; + end + col{i} = e; +end +end + +function col = numColumn(s, field, n) +col = nan(n, 1); +if isfield(s, field) && ~isempty(s.(field)) + if numel(s.(field)) ~= n + error('saveKcatSource:fields', 'kcatList.%s and kcatList.rxns differ in length.', field) + end + col = double(s.(field)(:)); +end +end + +function t = formatNumbers(x) +% formatKcatNumber for a vector, in a few vectorized passes. +x = x(:); +t = strings(numel(x), 1); +whole = isfinite(x) & x == round(x) & abs(x) < 1e15; +t(whole) = compose('%d', x(whole)); +t(isinf(x) & x > 0) = "Inf"; +t(isinf(x) & x < 0) = "-Inf"; +open = find(isfinite(x) & ~whole); +for digits = 1:17 + if isempty(open) + break + end + candidate = compose(sprintf('%%.%dg', digits), x(open)); + hit = str2double(candidate) == x(open); + t(open(hit)) = candidate(hit); + open = open(~hit); +end +end + +function t = cleanText(t) +t = regexprep(char(t), '[\t\r\n]', ' '); +end + +function v = geckoVersion() +v = ''; +file = fullfile(findGECKOroot, 'version.txt'); +if isfile(file) + v = strtrim(fileread(file)); +end +end diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index fad92a5a..56eca3ff 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -15,6 +15,13 @@ gather_kcats/readDLKcatOutput.m tcm0001, tcm0002 - The test doesn't really check gather_kcats/runDLKcat.m tcm0001, tcm0002 - The test doesn't really check the function, only that we get something gather_kcats/assignKcatValues.m tc0011 - we did not test all possible parameterizations here gather_kcats/writeDLKcatInput.m partly by tc0011 +gather_kcats/getEcRxnKeys.m tc0070, tc0076 +gather_kcats/saveKcatSource.m tc0071, tc0073, tc0076, tc0078 +gather_kcats/readKcatSource.m tc0071, tc0077 +gather_kcats/loadKcatSource.m tc0071, tc0072, tc0073, tc0076, tc0077, tc0078 +gather_kcats/kcatsFromBrenda.m tc0074, tc0078 +gather_kcats/kcatsFromDLKcat.m tc0074, tc0078 - reads an existing DLKcat output; running DLKcat needs Docker +gather_kcats/kcatsFromOKP.m tc0075 - reads a stored result; submitting needs the OKP service gather_kcats/reviewKcatAssignment.m Has only been tested by running it manually, not tested by test cases gather_kcats/ecStatsFromBrenda.m Not verified against real BRENDA data in this environment -- see its own TODO comment get_enzyme_data/calculateMW.m Not explicitly tested - The code has existed for a long time diff --git a/test/unit_tests/ecTestGEM/data/kcatSourceFixture.tsv b/test/unit_tests/ecTestGEM/data/kcatSourceFixture.tsv new file mode 100644 index 00000000..4cfdc9f1 --- /dev/null +++ b/test/unit_tests/ecTestGEM/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/test/unit_tests/kcatSourceTests.m b/test/unit_tests/kcatSourceTests.m new file mode 100644 index 00000000..4228cd74 --- /dev/null +++ b/test/unit_tests/kcatSourceTests.m @@ -0,0 +1,279 @@ +%run this test case with the command +%results = runtests('kcatSourceTests.m') + +function tests = kcatSourceTests + tests = functiontests(localfunctions); +end + +function testGetEcRxnKeys_tc0070(testCase) + % Full and light models give the same conventional reactions and + % enzyme sets. + [full, light] = buildKcatModels(); + [baseF, enzF] = getEcRxnKeys(full); + [baseL, enzL] = getEcRxnKeys(light); + verifyEqual(testCase, baseF(strcmp(full.ec.rxns, 'R2_EXP_1')), {'R2'}) + verifyEqual(testCase, baseF(strcmp(full.ec.rxns, 'R2_REV_EXP_2')), {'R2_REV'}) + verifyEqual(testCase, sort(unique(strcat(baseF, '|', enzF))), sort(unique(strcat(baseL, '|', enzL)))) + for i = 1:numel(enzF) + parts = strsplit(enzF{i}, ' + '); + verifyEqual(testCase, parts, sort(parts)) + end +end + +function testSaveLoadRoundTrip_tc0071(testCase) + % What saveKcatSource writes, loadKcatSource maps back to the same + % model.ec rows and values. + [full, ~, adapter] = buildKcatModels(); + list = fuzzyKcatMatching(full, [], adapter); + folder = tempname; + file = saveKcatSource(full, list, 'brenda', 'folder', folder, ... + 'parameters', struct('aggregation','max'), 'data', 'BRENDA release test'); + [records, meta] = readKcatSource(file); + verifyEqual(testCase, meta.source, 'brenda') + verifyEqual(testCase, meta.kind, 'gathered') + verifyEqual(testCase, meta.data, 'BRENDA release test') + verifyEqual(testCase, meta.parameters, '{"aggregation":"max"}') + keep = list.kcats > 0; + verifyEqual(testCase, numel(records.kcat), sum(keep)) + loaded = loadKcatSource(full, file); + verifyEqual(testCase, sort(loaded.rxns), sort(list.rxns(keep))) + [~, a] = sort(list.rxns(keep)); [~, b] = sort(loaded.rxns); + kept = find(keep); + verifyEqual(testCase, loaded.kcats(b), list.kcats(kept(a))) + verifyEqual(testCase, loaded.wildcardLvl(b), list.wildcardLvl(kept(a))) + verifyEqual(testCase, loaded.origin(b), list.origin(kept(a))) + verifyEqual(testCase, loaded.source, 'brenda') + subs = list.substrates(kept(a)); + verifyEqual(testCase, loaded.substrates(b), cellfun(@(c) strjoin(c, '; '), subs, 'UniformOutput', false)) + verifyEqual(testCase, loaded.eccodes(b), list.eccodes(kept(a))) + % Rows are sorted, so a rebuilt model writes the same file. + sorted = sortrows([records.rxn records.enzymes records.substrate records.source]); + verifyEqual(testCase, [records.rxn records.enzymes records.substrate records.source], sorted) +end + +function testFullFileOnLightModel_tc0072(testCase) + % A file written from the full model applies to the light model. + [full, light, adapter] = buildKcatModels(); + list = fuzzyKcatMatching(full, [], adapter); + file = saveKcatSource(full, list, 'brenda', 'folder', tempname); + [loaded, report] = loadKcatSource(light, file); + verifyEmpty(testCase, report.unmatched) + [baseL, enzL] = getEcRxnKeys(light); + [records] = readKcatSource(file); + for r = 1:numel(records.kcat) + hit = strcmp(baseL, records.rxn{r}) & strcmp(enzL, records.enzymes{r}); + verifyEqual(testCase, sort(loaded.rxns(loaded.line == records.line(r))), sort(light.ec.rxns(hit))) + end +end + +function testMatchingRules_tc0073(testCase) + % Rows with only enzymes or only a reaction, unmatched and partial rows. + [full, ~, adapter] = buildKcatModels(); + [base, enz] = getEcRxnKeys(full); + multi = find(contains(enz, ' + '), 1); + parts = strsplit(enz{multi}, ' + '); + folder = tempname; mkdir(folder); + file = fullfile(folder, 'custom.tsv'); + fid = fopen(file, 'w'); + fprintf(fid, '# source: custom\n# kind: curated\n'); + fprintf(fid, 'rxn\tenzymes\tsubstrate\tkcat\tlog10_sd\tsource\tdetail\tec_code\tnote\n'); + fprintf(fid, '\t%s\t\t10\t\tcustom\t\t\tenzymes only, reversed order\n', strjoin(fliplr(parts), ' + ')); + fprintf(fid, '%s\t\t\t20\t\tcustom\t\t\treaction only\n', base{multi}); + fprintf(fid, 'R999\tX1\t\t30\t\tcustom\t\t\tunmatched\n'); + fprintf(fid, '%s\t%s + XNEW\t\t40\t\tcustom\t\t\tpartial\n', base{multi}, enz{multi}); + fclose(fid); + [loaded, report] = loadKcatSource(full, file); + verifyEqual(testCase, sort(loaded.rxns(loaded.line == 4)), sort(full.ec.rxns(strcmp(enz, enz{multi})))) + verifyEqual(testCase, sort(loaded.rxns(loaded.line == 5)), sort(full.ec.rxns(strcmp(base, base{multi})))) + verifyEqual(testCase, report.unmatched.line, [6; 7]) + verifyEqual(testCase, unique(report.partial.line), 7) + % Curated files are never written: by name, by 'kind: curated' and by + % 'source: custom', in any case, and the file is left as it was. + list = fuzzyKcatMatching(full, [], adapter); + before = fileread(file); + verifyError(testCase, @() saveKcatSource(full, list, 'custom', 'folder', folder), 'saveKcatSource:curated') + verifyEqual(testCase, fileread(file), before) + writeLines(fullfile(folder, 'mine.tsv'), {'# source: Custom', 'rxn kcat source'}); + verifyError(testCase, @() saveKcatSource(full, list, 'mine', 'folder', folder), 'saveKcatSource:curated') + writeLines(fullfile(folder, 'mine2.tsv'), {'# kind: Curated', 'rxn kcat source'}); + verifyError(testCase, @() saveKcatSource(full, list, 'mine2', 'folder', folder), 'saveKcatSource:curated') +end + +function testKcatsFromBrendaAndDLKcat_tc0074(testCase) + % The source functions write their files into data/kcats. + [full, ~, adapter] = buildKcatModels(); + project = tempname; + copyfile(fullfile(adapter.params.path, 'data'), fullfile(project, 'data')); + adapter.params.path = project; + list = kcatsFromBrenda(full, 'aggregation', 'median', 'modelAdapter', adapter); + [records, meta] = readKcatSource(fullfile(project, 'data', 'kcats', 'brenda.tsv')); + verifyEqual(testCase, meta.generator(1:15), 'kcatsFromBrenda') + verifyEqual(testCase, meta.parameters, '{"aggregation":"median","force_wildcard_level":0}') + verifyEqual(testCase, meta.data, 'BRENDA release test-fixture') + verifyEqual(testCase, numel(records.kcat), sum(list.kcats > 0)) + + % DLKcat output for this model, read without running DLKcat. + rows = {'R2_EXP_1','G1','m1','C','MRAL','12.5'; 'R3','G4','m1','C','MDFM','3.25'}; + fid = fopen(fullfile(project, 'data', 'DLKcat.tsv'), 'w'); + fprintf(fid, 'rxns\tgenes\tsubstrates\tsmiles\tsequence\tkcat\n'); + for i = 1:size(rows, 1) + fprintf(fid, '%s\t%s\t%s\t%s\t%s\t%s\n', rows{i,:}); + end + fclose(fid); + kcatsFromDLKcat(full, 'modelAdapter', adapter); + [records, meta] = readKcatSource(fullfile(project, 'data', 'kcats', 'dlkcat.tsv')); + verifyEqual(testCase, meta.source, 'dlkcat') + verifyEqual(testCase, sort(records.kcat), [3.25; 12.5]) + verifyEqual(testCase, unique(records.source), {'DLKcat'}) +end + +function testKcatsFromOKPStored_tc0075(testCase) + % Stored OKP results: a method's own result of its current job, or a + % result stored by fetchOpenKineticsPredictor alone when its predictions + % are of the requested method. + [full, ~, adapter] = buildKcatModels(); + full.metSmiles = repmat({''}, numel(full.mets), 1); + full.metSmiles(strcmp(full.mets, 'm1c')) = {'C(C1C)O'}; + project = tempname; + copyfile(fullfile(adapter.params.path, 'data'), fullfile(project, 'data')); + adapter.params.path = project; + data = fullfile(project, 'data'); + okpCsv = {'kcat (1/s),Source kcat,Extra Info kcat,Protein Sequence,Substrate', ... + '12.5,Prediction from CataPro,,MRAL,C(C1C)O'}; + writeLines(fullfile(data, 'OKP_output.csv'), okpCsv); + writeLines(fullfile(data, 'OKP_job.txt'), {'jobId: abc123', 'method: CataPro'}); + [list, done] = kcatsFromOKP(full, 'method', "CataPro", 'useStored', true, 'modelAdapter', adapter); + verifyTrue(testCase, done) + [records, meta] = readKcatSource(fullfile(data, 'kcats', 'okp_catapro.tsv')); + verifyEqual(testCase, meta.data, 'OKP job abc123') + verifyEqual(testCase, numel(records.kcat), sum(list.kcats > 0)) + verifyEqual(testCase, unique(records.source), {'OKP-CataPro'}) + verifyTrue(testCase, isfile(fullfile(data, 'OKP_output_catapro.csv'))) + + % The shared job file names UniKP, but the shared result holds CataPro + % predictions (fetched after UniKP was submitted): refused. + writeLines(fullfile(data, 'OKP_job.txt'), {'jobId: uni1', 'method: UniKP'}); + verifyError(testCase, @() kcatsFromOKP(full, 'method', 'UniKP', 'useStored', true, ... + 'modelAdapter', adapter), 'kcatsFromOKP:noStored') + verifyFalse(testCase, isfile(fullfile(data, 'kcats', 'okp_unikp.tsv'))) + % A method's stored result of an earlier job is refused too. + writeLines(fullfile(data, 'OKP_job_kinform_h.txt'), {'jobId: new2', 'method: KinForm-H', 'resultOf: old1'}); + writeLines(fullfile(data, 'OKP_output_kinform_h.csv'), okpCsv); + verifyError(testCase, @() kcatsFromOKP(full, 'method', 'KinForm-H', 'useStored', true, ... + 'modelAdapter', adapter), 'kcatsFromOKP:noStored') + % ecRxns when no job is submitted, and a misspelled option, are errors. + verifyError(testCase, @() kcatsFromOKP(full, 'method', 'CataPro', 'useStored', true, ... + 'ecRxns', true(numel(full.ec.rxns), 1), 'modelAdapter', adapter), 'kcatsFromOKP:ecRxns') + verifyError(testCase, @() kcatsFromOKP(full, 'method', 'UniKP', 'wiat', true, ... + 'modelAdapter', adapter), 'kcatsFromOKP:arguments') +end + +function testGuardsNumbersAndJson_tc0078(testCase) + [full, ~, adapter] = buildKcatModels(); + project = tempname; + copyfile(fullfile(adapter.params.path, 'data'), fullfile(project, 'data')); + adapter.params.path = project; + folder = fullfile(project, 'data', 'kcats'); + mkdir(folder); + % The curated guard reads quoted header lines and ignores the case of keys. + writeLines(fullfile(folder, 'a.tsv'), {'"# kind: curated"', 'rxn kcat source'}); + writeLines(fullfile(folder, 'b.tsv'), {'# Kind: Curated', 'rxn kcat source'}); + writeLines(fullfile(folder, 'c.tsv'), {'# SOURCE: custom', 'rxn kcat source'}); + list = fuzzyKcatMatching(full, [], adapter); + for f = {'a','b','c'} + verifyError(testCase, @() saveKcatSource(full, list, f{1}, 'modelAdapter', adapter), 'saveKcatSource:curated') + end + % The source functions check the file before they gather. + verifyError(testCase, @() kcatsFromBrenda(full, 'saveAs', "custom", 'modelAdapter', adapter), 'saveKcatSource:curated') + verifyError(testCase, @() kcatsFromDLKcat(full, 'ecRxns', true(numel(full.ec.rxns), 1), ... + 'modelAdapter', adapter), 'kcatsFromDLKcat:ecRxns') + + % Infinite, NaN and non-positive kcats are not written, and not applied. + list.kcats(1) = Inf; list.kcats(2) = NaN; list.kcats(3) = -2; + file = saveKcatSource(full, list, "numbers", 'modelAdapter', adapter, 'parameters', ... + struct('z', 1, 'a', struct('y', ['caf' char(233)], 'b', true), 'n', NaN, ... + 'f', 1/3, 'g', 1e-7, 'h', 1e15, 'k', 1e20, 'l', {{5}}, 's', 'tab quote"')); + [records, meta] = readKcatSource(file); + verifyTrue(testCase, all(isfinite(records.kcat) & records.kcat > 0)) + % The same text as geckopy's test_parameters_json_matches_gecko. + verifyEqual(testCase, meta.parameters, ['{"a":{"b":true,"y":"caf' char(233) '"},"f":0.3333333333333333,' ... + '"g":1e-07,"h":1e+15,"k":1e+20,"l":[5],"n":null,"s":"tab\tquote\"","z":1}']) + writeLines(fullfile(folder, 'inf.tsv'), {'rxn enzymes kcat source', 'R3 P4 Inf custom', 'R3 P4 1e999 custom'}); + [loaded, report] = loadKcatSource(full, fullfile(folder, 'inf.tsv')); + verifyEmpty(testCase, loaded.rxns) + verifyEqual(testCase, report.invalid.line, [2; 3]) + + % Rows are sorted by every column, in the same order as geckopy. + writeLines(fullfile(folder, 'unsorted.tsv'), {'rxn enzymes substrate kcat source detail', ... + 'R3 P4 m1 2 b ', 'R3 P4 m1 2 a wc=1', 'R3 P4 5 a ', 'R2_REV P3 m2 1 a '}); + l2 = loadKcatSource(full, fullfile(folder, 'unsorted.tsv')); + file = saveKcatSource(full, l2, 'sorted', 'modelAdapter', adapter); + records = readKcatSource(file); + verifyEqual(testCase, records.rxn', {'R2_REV','R3','R3','R3'}) + verifyEqual(testCase, records.substrate', {'m2','','m1','m1'}) + verifyEqual(testCase, records.source', {'a','a','a','b'}) +end + +function testEnzymelessRowsAndLightToFull_tc0076(testCase) + % Rows of reactions without enzymes are not written, so they cannot + % spread to their isozymes; a light-model file applies to the full model. + [full, light, adapter] = buildKcatModels(); + list = fuzzyKcatMatching(full, [], adapter); + i = find(strcmp(full.ec.rxns, 'R2_EXP_2')); + full.ec.rxnEnzMat(i, :) = 0; + file = saveKcatSource(full, list, 'brenda', 'folder', tempname); + loaded = loadKcatSource(full, file); + verifyFalse(testCase, any(strcmp(loaded.rxns, 'R2_EXP_2'))) + verifyEqual(testCase, sum(strcmp(loaded.rxns, 'R2_EXP_1')), 1) + + listL = fuzzyKcatMatching(light, [], adapter); + fileL = saveKcatSource(light, listL, 'brenda', 'folder', tempname); + fullAgain = buildKcatModels(); + [~, report] = loadKcatSource(fullAgain, fileL); + verifyEmpty(testCase, report.unmatched) +end + +function testSharedFixture_tc0077(testCase) + % The fixture file, read the same way by geckopy: 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. + file = fullfile(findGECKOroot, 'test', 'unit_tests', 'ecTestGEM', 'data', 'kcatSourceFixture.tsv'); + [records, meta] = readKcatSource(file); + verifyEqual(testCase, meta.source, 'fixture') + verifyEqual(testCase, meta.parameters, '{"aggregation":"median"}') + verifyEqual(testCase, meta.kind, 'curated') + verifyEqual(testCase, records.line', [7 9 11 12 13 14 15]) + verifyEqual(testCase, records.rxn', {'R2','R3','R5','R2_REV','','R9','R9'}) + verifyEqual(testCase, records.substrate{2}, '2,3-bisphospho-D-glycerate') + verifyEqual(testCase, records.kcat', [12.5 1e-05 NaN 7 3 -1 4]) + verifyEqual(testCase, records.log10_sd(1), 0.3) + + [full] = buildKcatModels(); + [loaded, report] = loadKcatSource(full, file); + verifyEqual(testCase, loaded.rxns(loaded.line == 7), {'R2_EXP_1'}) + verifyEqual(testCase, sort(loaded.rxns(loaded.line == 12)), {'R2_REV_EXP_1';'R2_REV_EXP_2'}) + verifyEqual(testCase, sort(loaded.rxns(loaded.line == 13)), {'R2_EXP_2';'R2_REV_EXP_2'}) + verifyEqual(testCase, report.invalid.line, [11; 14]) + verifyEqual(testCase, report.unmatched.line, 15) + + % A row with more filled fields than columns is an error. + bad = [tempname '.tsv']; + writeLines(bad, {'rxn kcat source', 'R2 1 brenda extra'}); + verifyError(testCase, @() readKcatSource(bad), 'readKcatSource:fields') +end + +function writeLines(file, lines) +fid = fopen(file, 'w'); +fprintf(fid, '%s\n', lines{:}); +fclose(fid); +end + +function [full, light, adapter] = buildKcatModels() + geckoPath = findGECKOroot; + adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); + model = getGeckoTestModel(); + full = makeEcModel(model, false, adapter); + full = getECfromGEM(full); + light = makeEcModel(model, true, adapter); + light = getECfromGEM(light); +end From 7c316352bb3ba75681decbec71a3b5a85e7b4220 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 07:40:04 +0200 Subject: [PATCH 2/3] Assign kcats from the kcat source files, as model_adapter.toml sets collectKcats decides one kcat per reaction from the files in data/kcats, in the stages of [kcat_assignment] (a preset, or [[kcat_assignment.stage]] tables), without changing the model, and writes output/_kcatDecisions.tsv. assignKcats applies the decisions and the kcat constraints. updateKcatSources gathers the source files the configuration needs that are missing, or outdated against [kcat_sources.]. The default preset, gecko_median, orders custom kcats, the BRENDA median of matches without EC wildcards at origins 1 to 4, DLKcat, then the other BRENDA matches. gecko_classic keeps the GECKO 3 order with the BRENDA maximum; the full_ecModel tutorial uses it and gives the same model as before. --- src/gather_kcats/assignKcats.m | 42 ++++ src/gather_kcats/collectKcats.m | 236 ++++++++++++++++++ .../private/kcatAssignmentConfig.m | 92 +++++++ src/gather_kcats/updateKcatSources.m | 206 +++++++++++++++ src/model_adapter/adapterParameters.tsv | 5 + src/model_adapter/adapterTemplate.toml | 15 ++ test/VerificationMatrix.tsv | 3 + test/unit_tests/kcatAssignmentTests.m | 150 +++++++++++ tutorials/full_ecModel/model_adapter.toml | 6 + tutorials/full_ecModel/protocol.m | 94 +++---- 10 files changed, 803 insertions(+), 46 deletions(-) create mode 100644 src/gather_kcats/assignKcats.m create mode 100644 src/gather_kcats/collectKcats.m create mode 100644 src/gather_kcats/private/kcatAssignmentConfig.m create mode 100644 src/gather_kcats/updateKcatSources.m create mode 100644 test/unit_tests/kcatAssignmentTests.m diff --git a/src/gather_kcats/assignKcats.m b/src/gather_kcats/assignKcats.m new file mode 100644 index 00000000..ca945081 --- /dev/null +++ b/src/gather_kcats/assignKcats.m @@ -0,0 +1,42 @@ +function [model, decisions, report] = assignKcats(model, varargin) +% assignKcats Assign kcats from the kcat source files and apply them. +% +% Decides one kcat per reaction from the kcat source files in data/kcats, +% as [kcat_assignment] in model_adapter.toml sets (see collectKcats), writes +% model.ec.kcat and model.ec.source, and applies the kcat constraints +% (applyKcatConstraints). Any kcats the model had are replaced: the source +% files are the input, and nothing in the model is written back to them. +% +% Parameters +% ---------- +% model : struct +% an ecModel in GECKO 3 format (with ecModel.ec structure). +% +% Name-Value Arguments +% -------------------- +% modelAdapter : ModelAdapter +% a loaded model adapter (default: the current default model adapter). +% writeReport : logical +% write the decisions to output/_kcatDecisions.tsv (default +% true). +% +% Returns +% ------- +% model : struct +% the ecModel with the assigned kcats applied. +% decisions, report +% as returned by collectKcats. +% +% Examples +% -------- +% updateKcatSources(ecModel); % gather what the configuration needs +% ecModel = assignKcats(ecModel); +% +% See also +% -------- +% collectKcats, updateKcatSources, applyKcatConstraints + +checkNameValue(varargin, {'modelAdapter','writeReport'}, 'assignKcats'); +[decisions, report, model] = collectKcats(model, varargin{:}); +model = applyKcatConstraints(model); +end diff --git a/src/gather_kcats/collectKcats.m b/src/gather_kcats/collectKcats.m new file mode 100644 index 00000000..5d860243 --- /dev/null +++ b/src/gather_kcats/collectKcats.m @@ -0,0 +1,236 @@ +function [decisions, report, ecModel] = collectKcats(model, varargin) +% collectKcats Decide one kcat per reaction from the kcat source files. +% +% Combines the kcat source files in data/kcats (see kcatsFromBrenda, +% kcatsFromDLKcat, kcatsFromOKP and saveKcatSource) as [kcat_assignment] in +% model_adapter.toml sets, without changing the model given. assignKcats +% applies the decisions. +% +% The order is a preset or a list of stages: +% +% [kcat_assignment] +% preset = "gecko_median" # the default, or gecko_classic +% isozymes = true # fill from isozymes after the stages +% standard = true # standard kcat for the rest +% +% # instead of a preset: +% [[kcat_assignment.stage]] +% action = "select" +% sources = ["custom", "database_top", "dlkcat", "database_bottom"] +% criteria = "max" # among several values of one source +% top_origin_limit = 4 +% bottom_origin_limit = 6 +% wildcard_limit = 3 +% +% A select stage decides each reaction that is still open by the first +% 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). +% +% After the stages: +% - custom kcats of data/customKcats.tsv are applied (applyCustomKcats), +% when a stage lists custom but there is no data/kcats/custom.tsv; +% - with isozymes, reactions without a kcat take one from their isozymes +% (getKcatAcrossIsozymes); +% - with standard, the remaining reactions get the standard kcat +% (getStandardKcat). +% +% A source file a stage names that does not exist is reported and left out. +% +% Parameters +% ---------- +% model : struct +% an ecModel in GECKO 3 format (with ecModel.ec structure). +% +% Name-Value Arguments +% -------------------- +% modelAdapter : ModelAdapter +% a loaded model adapter (default: the current default model adapter). +% writeReport : logical +% write the decisions to output/_kcatDecisions.tsv (default +% true). +% +% Returns +% ------- +% decisions : table +% one row per model.ec reaction: ecRxn, stage (the stage number, or +% custom, isozymes, standard, or empty when no kcat was found), source +% (as in model.ec.source), kcat, candidates (the number of values the +% deciding source had for the reaction) and spreadLog10 (their range, +% in log10 units). +% report : struct +% - preset: the preset used, or '' with stages; +% - missing: source files the stages name that do not exist; +% - sources: per source file loaded, its unmatched, partial and invalid +% rows (see loadKcatSource). +% ecModel : struct +% the model with the decided kcats, before applyKcatConstraints; this +% is what assignKcats returns, after applying the constraints. +% +% Examples +% -------- +% decisions = collectKcats(ecModel); +% decisions(strcmp(decisions.stage, ''), :) % reactions without a kcat +% +% See also +% -------- +% assignKcats, updateKcatSources, loadKcatSource, mergeKcats, assignKcatValues + +checkNameValue(varargin, {'modelAdapter','writeReport'}, 'collectKcats'); +p = parseGECKOargs(varargin, {'modelAdapter', []; 'writeReport', true}); +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 +params = modelAdapter.params; +config = kcatAssignmentConfig(params); +kcatDir = fullfile(params.path, 'data', 'kcats'); + +ecModel = model; +nEc = numel(ecModel.ec.rxns); +ecModel.ec.kcat(:) = 0; +ecModel.ec.source(:) = {''}; +stage = repmat({''}, nEc, 1); +candidates = zeros(nEc, 1); +spread = nan(nEc, 1); + +report.preset = config.preset; +report.missing = {}; +report.sources = struct(); +loaded = struct(); +legacyCustom = false; + + function list = loadSource(name) + % The kcat list of a source file, or [] when it does not exist. + list = []; + file = fullfile(kcatDir, [name '.tsv']); + if isfile(file) + [list, rep] = loadKcatSource(ecModel, file); + report.sources.(name) = rep; + elseif strcmp(name, 'custom') && isfile(fullfile(params.path, 'data', 'customKcats.tsv')) + legacyCustom = true; + else + report.missing{end+1} = name; + printOrange(sprintf(['No kcat source file %s; it is left out. ' ... + 'updateKcatSources gathers it.\n'], file)); + end + end + +for s = 1:numel(config.stages) + st = config.stages(s); + lists = {}; + priority = {}; + for t = 1:numel(st.sources) + token = st.sources{t}; + name = sourceFile(token); + if ~isfield(loaded, name) + loaded.(name) = loadSource(name); + end + list = loaded.(name); + if isempty(list) + continue + end + if ~any(cellfun(@(l) isequal(l, list), lists)) + lists{end+1} = list; %#ok + end + if startsWith(token, 'database_') + priority{end+1} = token; %#ok + else + % The labels of the file's rows, at the file's place in the order. + labels = unique(cellfun(@normalizeSource, list.kcatSource, 'UniformOutput', false), 'stable'); + priority = [priority, labels(~ismember(labels, {'brenda','sabio_rk'}))']; %#ok + end + end + if isempty(lists) || isempty(priority) + continue + end + merged = mergeKcats(lists, priority, 'topOriginLimit', st.topOriginLimit, ... + 'bottomOriginLimit', st.bottomOriginLimit, 'wildcardLimit', st.wildcardLimit); + if isempty(merged.rxns) + continue + end + open = ecModel.ec.kcat == 0; + [ecModel, idx] = assignKcatValues(ecModel, merged, 'criteria', st.criteria, 'overwrite', 'false'); + idx = idx(open(idx)); + stage(idx) = {sprintf('%d', s)}; + [~, where] = ismember(merged.rxns, ecModel.ec.rxns); + valid = merged.kcats > 0; + for i = idx(:)' + values = merged.kcats(where == i & valid); + candidates(i) = numel(values); + spread(i) = log10(max(values)) - log10(min(values)); + end +end + +if legacyCustom + before = ecModel.ec.kcat; + ecModel = applyCustomKcats(ecModel, 'modelAdapter', modelAdapter); + changed = ecModel.ec.kcat ~= before; + stage(changed) = {'custom'}; + candidates(changed) = 1; + spread(changed) = 0; +end +if config.isozymes && ~ecModel.ec.geckoLight + before = ecModel.ec.kcat == 0; + ecModel = getKcatAcrossIsozymes(ecModel); + stage(before & ecModel.ec.kcat > 0) = {'isozymes'}; +end +if config.standard + ecModel = getStandardKcat(ecModel, 'modelAdapter', modelAdapter); + n = numel(ecModel.ec.rxns); + stage(end+1:n) = {''}; + candidates(end+1:n) = 0; + spread(end+1:n) = NaN; + stage(strcmp(ecModel.ec.source, 'standard') & ecModel.ec.kcat > 0) = {'standard'}; +end + +decisions = table(ecModel.ec.rxns(:), stage(:), ecModel.ec.source(:), ecModel.ec.kcat(:), ... + candidates(:), spread(:), 'VariableNames', {'ecRxn','stage','source','kcat','candidates','spreadLog10'}); +if p.writeReport + outDir = fullfile(params.path, 'output'); + if ~isfolder(outDir) + mkdir(outDir); + end + id = 'ecModel'; + if isfield(model, 'id') && ~isempty(model.id) + id = char(model.id); + end + writetable(decisions, fullfile(outDir, [id '_kcatDecisions.tsv']), 'FileType', 'text', ... + 'Delimiter', '\t'); +end +nOpen = sum(strcmp(stage, '')); +fprintf('Kcats decided for %d of %d reactions.\n', numel(stage) - nOpen, numel(stage)); +end + +function s = normalizeSource(label) +% As mergeKcats folds source labels: lower case, other characters as '_'. +s = lower(strtrim(char(label))); +s = regexprep(s, '[^0-9a-z]+', '_'); +s = regexprep(s, '^_+|_+$', ''); +end + +function name = sourceFile(token) +% The source file a stage token reads: the BRENDA tiers read brenda.tsv. +if startsWith(token, 'database_') + name = 'brenda'; +else + name = token; +end +end diff --git a/src/gather_kcats/private/kcatAssignmentConfig.m b/src/gather_kcats/private/kcatAssignmentConfig.m new file mode 100644 index 00000000..c0592e3b --- /dev/null +++ b/src/gather_kcats/private/kcatAssignmentConfig.m @@ -0,0 +1,92 @@ +function config = kcatAssignmentConfig(params) +% kcatAssignmentConfig The kcat assignment settings of a model adapter. +% +% 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. Without a preset or stages the preset +% is gecko_median. +% +% Returns a struct with: +% - preset: the preset name, or '' when stages are given; +% - isozymes, standard: logical; +% - stages: struct array with action, sources (cell), criteria, +% topOriginLimit, bottomOriginLimit, wildcardLimit; +% - sources: struct, per source name its options (from [kcat_sources]). +% geckopy's kcat_assignment_config gives the same settings. + +a = struct(); +if isfield(params, 'kcat_assignment') && isstruct(params.kcat_assignment) + a = params.kcat_assignment; +end +config.preset = ''; +config.isozymes = getOr(a, 'isozymes', true); +config.standard = getOr(a, 'standard', true); +config.sources = struct(); +hasPreset = isfield(a, 'preset') && ~isempty(a.preset); +hasStages = isfield(a, 'stage') && ~isempty(a.stage); +if hasPreset && hasStages + error('kcatAssignment:config', ['[kcat_assignment] gives both a preset and stages; ' ... + 'give one of the two.']) +end +if hasStages + config.stages = arrayfun(@readStage, a.stage(:)'); +else + preset = 'gecko_median'; + if hasPreset + preset = char(a.preset); + end + [config.stages, config.sources] = presetStages(preset); + config.preset = preset; +end +if isfield(params, 'kcat_sources') && isstruct(params.kcat_sources) + names = fieldnames(params.kcat_sources); + for k = 1:numel(names) + config.sources.(names{k}) = params.kcat_sources.(names{k}); + end +end +end + +function [stages, sources] = presetStages(preset) +% The presets of kcat assignment; geckopy defines the same. +select = @(srcs, top) struct('action', 'select', 'sources', {srcs}, 'criteria', 'max', ... + 'topOriginLimit', top, 'bottomOriginLimit', 6, 'wildcardLimit', 3); +order = {'custom', 'database_top', 'dlkcat', 'database_bottom'}; +switch preset + case 'gecko_classic' + stages = select(order, 6); + sources.brenda = struct('aggregation', 'max'); + case 'gecko_median' + stages = select(order, 4); + sources.brenda = struct('aggregation', 'median'); + otherwise + error('kcatAssignment:config', 'Unknown kcat assignment preset "%s". Presets: %s.', ... + preset, 'gecko_classic, gecko_median') +end +end + +function stage = readStage(s) +stage.action = getOr(s, 'action', 'select'); +if ~strcmp(stage.action, 'select') + error('kcatAssignment:config', ['Unknown stage action "%s"; this GECKO version knows ' ... + '"select".'], stage.action) +end +sources = getOr(s, 'sources', {}); +if ischar(sources) + sources = {sources}; +end +if isempty(sources) + error('kcatAssignment:config', 'A [[kcat_assignment.stage]] needs sources.') +end +stage.sources = sources(:)'; +stage.criteria = validatestring(getOr(s, 'criteria', 'max'), {'max','min','median','mean'}); +stage.topOriginLimit = getOr(s, 'top_origin_limit', 4); +stage.bottomOriginLimit = getOr(s, 'bottom_origin_limit', 6); +stage.wildcardLimit = getOr(s, 'wildcard_limit', 3); +end + +function v = getOr(s, field, default) +v = default; +if isstruct(s) && isfield(s, field) && ~isempty(s.(field)) + v = s.(field); +end +end diff --git a/src/gather_kcats/updateKcatSources.m b/src/gather_kcats/updateKcatSources.m new file mode 100644 index 00000000..e189a422 --- /dev/null +++ b/src/gather_kcats/updateKcatSources.m @@ -0,0 +1,206 @@ +function report = updateKcatSources(model, varargin) +% updateKcatSources Gather the kcat sources that the configuration needs. +% +% Lists the kcat source files that [kcat_assignment] in model_adapter.toml +% uses (see collectKcats), and gathers those that are missing or outdated +% with kcatsFromBrenda, kcatsFromDLKcat and kcatsFromOKP: +% +% - 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 GECKO's BRENDA files; +% - brenda is matched again; 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 the job runs; +% - custom.tsv and other files it cannot gather are left as they are. +% +% Requirements are checked first: Docker for a DLKcat run and an +% OpenKineticsPredictor API key for OKP. Calling it again is safe: files +% that are current are not gathered again. +% +% Per-source settings are read from model_adapter.toml, for instance: +% +% [kcat_sources.brenda] +% aggregation = "median" +% +% [kcat_sources.okp_catapro] +% method = "CataPro" +% +% Parameters +% ---------- +% model : struct +% an ecModel in GECKO 3 format (with ecModel.ec structure). +% +% Name-Value Arguments +% -------------------- +% dryRun : logical +% only report what would be gathered (default false). +% modelAdapter : ModelAdapter +% a loaded model adapter (default: the current default model adapter). +% +% Returns +% ------- +% report : table +% one row per source: source, status (current, gathered, pending, +% outdated, missing, failed or skipped) and note. +% +% Examples +% -------- +% updateKcatSources(ecModel, 'dryRun', true) +% updateKcatSources(ecModel); +% +% See also +% -------- +% collectKcats, assignKcats, kcatsFromBrenda, kcatsFromDLKcat, kcatsFromOKP + +checkNameValue(varargin, {'dryRun','modelAdapter'}, 'updateKcatSources'); +p = parseGECKOargs(varargin, {'dryRun', false; '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 +params = modelAdapter.params; +config = kcatAssignmentConfig(params); +kcatDir = fullfile(params.path, 'data', 'kcats'); + +names = {}; +for s = 1:numel(config.stages) + for t = 1:numel(config.stages(s).sources) + token = config.stages(s).sources{t}; + if startsWith(token, 'database_') + token = 'brenda'; + end + if ~any(strcmp(names, token)) + names{end+1} = token; %#ok + end + end +end + +status = cell(numel(names), 1); +note = cell(numel(names), 1); +for k = 1:numel(names) + name = names{k}; + file = fullfile(kcatDir, [name '.tsv']); + wanted = struct(); + if isfield(config.sources, name) + wanted = config.sources.(name); + end + [status{k}, note{k}] = check(name, file, wanted); + if any(strcmp(status{k}, {'missing','outdated'})) && ~p.dryRun + [status{k}, note{k}] = gather(model, name, wanted, modelAdapter); + end +end +report = table(names(:), status, note, 'VariableNames', {'source','status','note'}); +disp(report) + + function [st, nt] = check(name, file, wanted) + st = 'current'; nt = ''; + if strcmp(name, 'custom') + st = 'skipped'; nt = 'curated by hand'; + return + elseif ~any(strcmp(name, {'brenda','dlkcat'})) && ~startsWith(name, 'okp_') + st = 'skipped'; nt = 'not a source GECKO gathers'; + return + end + if ~isfile(file) + st = 'missing'; + return + end + [~, meta] = readKcatSource(file); + recorded = struct(); + if isfield(meta, 'parameters') && ~isempty(meta.parameters) + recorded = jsondecode(meta.parameters); + end + f = fieldnames(wanted); + for i = 1:numel(f) + if ~isfield(recorded, f{i}) || ~isequal(normalizeValue(recorded.(f{i})), normalizeValue(wanted.(f{i}))) + st = 'outdated'; + nt = sprintf('%s differs from model_adapter.toml', f{i}); + return + end + end + if strcmp(name, 'brenda') + release = brendaRelease(modelAdapter); + if isfield(meta, 'data') && ~isempty(release) && ~strcmp(meta.data, release) + st = 'outdated'; + nt = sprintf('made from %s; GECKO has %s', meta.data, release); + end + end + end +end + +function [st, nt] = gather(model, name, wanted, modelAdapter) +st = 'gathered'; nt = ''; +try + if strcmp(name, 'brenda') + args = {}; + if isfield(wanted, 'aggregation'), args = [args, {'aggregation', wanted.aggregation}]; end + if isfield(wanted, 'force_wildcard_level'), args = [args, {'forceWClvl', wanted.force_wildcard_level}]; end + kcatsFromBrenda(model, args{:}, 'modelAdapter', modelAdapter); + elseif strcmp(name, 'dlkcat') + run = ~isfile(fullfile(modelAdapter.params.path, 'data', 'DLKcat.tsv')); + if run + [status, ~] = system('docker --version'); + if status ~= 0 + st = 'failed'; nt = 'DLKcat needs Docker to run'; + return + end + end + kcatsFromDLKcat(model, 'run', run, 'modelAdapter', modelAdapter); + else % okp_ + method = okpMethod(name, wanted); + resolveOkpApiKey('', modelAdapter.params.path); % errors without a key + [~, done] = kcatsFromOKP(model, 'method', method, 'saveAs', name, 'modelAdapter', modelAdapter); + if ~done + st = 'pending'; nt = sprintf('OKP job for %s still running; call again later', method); + end + end +catch err + st = 'failed'; nt = err.message; +end +end + +function method = okpMethod(name, wanted) +% The OKP method of okp_: from [kcat_sources.okp_] method, +% or the known method whose name gives this file name. +if isfield(wanted, 'method') + method = char(wanted.method); + return +end +known = {'CataPro','CatPred','CatRange','DLKcat','EITLEM','KinForm-H','KinForm-L', ... + 'OmniESI','OmniESI-O2DENet','TurNuP','UniKP'}; +tags = cellfun(@(m) ['okp_' lower(regexprep(m, '[^A-Za-z0-9]+', '_'))], known, 'UniformOutput', false); +hit = strcmp(tags, name); +if ~any(hit) + error('updateKcatSources:okpMethod', ['Set the OKP method of %s in model_adapter.toml: ' ... + '[kcat_sources.%s] method = "..."'], name, name) +end +method = known{hit}; +end + +function v = normalizeValue(v) +if isstring(v) + v = char(v); +end +if ischar(v) + v = lower(v); +elseif islogical(v) + v = double(v); +end +end + +function release = brendaRelease(modelAdapter) +release = ''; +fid = fopen(fullfile(modelAdapter.getBrendaDBFolder(), 'kcat.tsv'), '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 diff --git a/src/model_adapter/adapterParameters.tsv b/src/model_adapter/adapterParameters.tsv index c4c45f80..7ad87b1c 100644 --- a/src/model_adapter/adapterParameters.tsv +++ b/src/model_adapter/adapterParameters.tsv @@ -32,5 +32,10 @@ okp.method okp.method "CataPro" OpenKineticsPredictor method okp.handle_long_sequences okp.handleLongSequences "truncate" How OpenKineticsPredictor handles sequences that are too long okp.include_similarity_columns okp.includeSimilarityColumns true Return similarity to the method's training data okp.canonicalize_substrates okp.canonicalizeSubstrates true Canonicalize substrate SMILES before prediction +kcat_assignment.preset Kcat assignment preset: gecko_median (used when neither preset nor stages are given) or gecko_classic +kcat_assignment.isozymes true Fill reactions without a kcat from their isozymes after the stages +kcat_assignment.standard true Give the remaining reactions the standard kcat +kcat_assignment.stage Stages instead of a preset: action, sources, criteria and BRENDA limits ([[kcat_assignment.stage]]) +kcat_sources Settings per kcat source file, e.g. [kcat_sources.brenda] aggregation matlab.make_model_anaerobic matlab.makeModelAnaerobic MATLAB function that constrains the ecModel to anaerobic growth matlab.change_protein_biomass matlab.changeProteinBiomass MATLAB function that changes the protein content of the biomass diff --git a/src/model_adapter/adapterTemplate.toml b/src/model_adapter/adapterTemplate.toml index f75c2143..d05d8d06 100644 --- a/src/model_adapter/adapterTemplate.toml +++ b/src/model_adapter/adapterTemplate.toml @@ -104,6 +104,21 @@ dlkcat = 0.4 brenda = 0.2 custom = 0.1 +# How kcats are assigned from the kcat source files in data/kcats (see +# collectKcats and assignKcats): a preset, or [[kcat_assignment.stage]] +# tables instead. gecko_median orders custom kcats, the BRENDA median of +# matches without EC wildcards or specific activities, DLKcat, then the +# other BRENDA matches; gecko_classic uses the BRENDA maximum and ranks the +# specific-activity matches above DLKcat, as the GECKO 3 tutorial. +[kcat_assignment] +preset = "gecko_median" +isozymes = true +standard = true + +# Settings per kcat source, used by updateKcatSources, for example: +# [kcat_sources.okp_catapro] +# method = "CataPro" + # Read by GECKO only. Names of MATLAB functions, for instance in the code # folder, used during kcat tuning: one that constrains the ecModel to # anaerobic growth, and one that changes the protein content of the diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 36118cbd..43edd03b 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -22,6 +22,9 @@ gather_kcats/loadKcatSource.m tc0071, tc0072, tc0073, tc0076, tc0077, tc0078 gather_kcats/kcatsFromBrenda.m tc0074, tc0078 gather_kcats/kcatsFromDLKcat.m tc0074, tc0078 - reads an existing DLKcat output; running DLKcat needs Docker gather_kcats/kcatsFromOKP.m tc0075 - reads a stored result; submitting needs the OKP service +gather_kcats/collectKcats.m tc0079, tc0080, tc0081, tc0083 +gather_kcats/assignKcats.m tc0079, tc0083 +gather_kcats/updateKcatSources.m tc0082 gather_kcats/reviewKcatAssignment.m Has only been tested by running it manually, not tested by test cases gather_kcats/ecStatsFromBrenda.m Not verified against real BRENDA data in this environment -- see its own TODO comment get_enzyme_data/calculateMW.m Not explicitly tested - The code has existed for a long time diff --git a/test/unit_tests/kcatAssignmentTests.m b/test/unit_tests/kcatAssignmentTests.m new file mode 100644 index 00000000..d119b085 --- /dev/null +++ b/test/unit_tests/kcatAssignmentTests.m @@ -0,0 +1,150 @@ +%run this test case with the command +%results = runtests('kcatAssignmentTests.m') + +function tests = kcatAssignmentTests + tests = functiontests(localfunctions); +end + +function testAssignKcatsReproducesChain_tc0079(testCase) + % The gecko_classic preset gives the same model as the steps it + % replaces: fuzzyKcatMatching, readDLKcatOutput, mergeKcats, + % assignKcatValues, applyCustomKcats, getKcatAcrossIsozymes, + % getStandardKcat and applyKcatConstraints. + [base, adapter, project] = buildProject(); + adapter.params.kcat_assignment = struct('preset', 'gecko_classic'); + fuzzy = fuzzyKcatMatching(base, [], adapter); + dlkcat = readDLKcatOutput(base, 'modelAdapter', adapter); + merged = mergeKcats({fuzzy, dlkcat}, {'database_top','dlkcat','database_bottom'}); + A = assignKcatValues(base, merged); + A = applyCustomKcats(A, 'modelAdapter', adapter); + A = getKcatAcrossIsozymes(A); + A = getStandardKcat(A, 'modelAdapter', adapter); + A = applyKcatConstraints(A); + + kcatsFromBrenda(base, 'modelAdapter', adapter); + kcatsFromDLKcat(base, 'modelAdapter', adapter); + [B, decisions, report] = assignKcats(base, 'modelAdapter', adapter); + verifyEqual(testCase, B.ec.rxns, A.ec.rxns) + verifyEqual(testCase, B.ec.kcat, A.ec.kcat) + verifyEqual(testCase, B.ec.source, A.ec.source) + verifyEqual(testCase, B.S, A.S) + verifyEqual(testCase, report.preset, 'gecko_classic') + verifyEqual(testCase, decisions.ecRxn, B.ec.rxns) + verifyTrue(testCase, isfile(fullfile(project, 'output', [char(base.id) '_kcatDecisions.tsv']))) +end + +function testCollectKcatsLeavesModel_tc0080(testCase) + % collectKcats decides without changing the model; a missing source is + % reported; custom.tsv decides first. + [base, adapter, project] = buildProject(); + kcatsFromBrenda(base, 'modelAdapter', adapter); + [decisions, report, decided] = collectKcats(base, 'modelAdapter', adapter, 'writeReport', false); + verifyEqual(testCase, base.ec.kcat, zeros(numel(base.ec.rxns), 1)) + verifyEqual(testCase, report.missing, {'dlkcat'}) + verifyTrue(testCase, any(decided.ec.kcat > 0)) + + [baseRxns, enz] = getEcRxnKeys(base); + i = find(strcmp(base.ec.rxns, 'R3')); + writeLines(fullfile(project, 'data', 'kcats', 'custom.tsv'), {'# source: custom', '# kind: curated', ... + 'rxn enzymes kcat source', sprintf('%s %s 777 custom', baseRxns{i}, enz{i})}); + decisions = collectKcats(base, 'modelAdapter', adapter, 'writeReport', false); + verifyEqual(testCase, decisions.kcat(i), 777) + verifyEqual(testCase, decisions.source{i}, 'custom') + verifyEqual(testCase, decisions.stage{i}, '1') +end + +function testStagesFromToml_tc0081(testCase) + % Stages from model_adapter.toml, and their errors. + [base, adapter, project] = buildProject(); + kcatsFromBrenda(base, 'modelAdapter', adapter); + kcatsFromDLKcat(base, 'modelAdapter', adapter); + adapter.params.kcat_assignment = struct('stage', struct('action', 'select', ... + 'sources', {{'dlkcat'}}, 'criteria', 'min'), 'isozymes', false, 'standard', false); + decisions = collectKcats(base, 'modelAdapter', adapter, 'writeReport', false); + verifyEqual(testCase, unique(decisions.source(decisions.kcat > 0)), {'DLKcat'}) + verifyEqual(testCase, decisions.kcat(strcmp(decisions.ecRxn, 'R2_EXP_1')), 6.5) + adapter.params.kcat_assignment.preset = 'gecko_median'; + verifyError(testCase, @() collectKcats(base, 'modelAdapter', adapter), 'kcatAssignment:config') + adapter.params.kcat_assignment = struct('preset', 'gecko_fast'); + verifyError(testCase, @() collectKcats(base, 'modelAdapter', adapter), 'kcatAssignment:config') + + % The same stages read from a TOML file. + toml = fullfile(project, 'model_adapter.toml'); + writeLines(toml, {'conv_gem = "models/testModel.xml"', 'org_name = "testus testus"', ... + '[kcat_assignment]', 'isozymes = false', '[[kcat_assignment.stage]]', 'sources = ["dlkcat"]', ... + 'criteria = "min"', '[kcat_sources.brenda]', 'aggregation = "median"'}); + params = loadAdapterParams(toml); + verifyEqual(testCase, params.kcat_assignment.stage.sources, {'dlkcat'}) + verifyEqual(testCase, params.kcat_assignment.standard, true) + verifyEqual(testCase, params.kcat_sources.brenda.aggregation, 'median') +end + +function testUpdateKcatSources_tc0082(testCase) + % Missing and outdated sources are gathered; current ones are not. + [base, adapter, project] = buildProject(); + adapter.params.kcat_assignment = struct('preset', 'gecko_median'); + report = updateKcatSources(base, 'dryRun', true, 'modelAdapter', adapter); + verifyEqual(testCase, report.source', {'custom','brenda','dlkcat'}) + verifyEqual(testCase, report.status', {'skipped','missing','missing'}) + verifyFalse(testCase, isfolder(fullfile(project, 'data', 'kcats'))) + report = updateKcatSources(base, 'modelAdapter', adapter); + verifyEqual(testCase, report.status', {'skipped','gathered','gathered'}) + [~, meta] = readKcatSource(fullfile(project, 'data', 'kcats', 'brenda.tsv')); + verifyEqual(testCase, meta.parameters, '{"aggregation":"median","force_wildcard_level":0}') + report = updateKcatSources(base, 'modelAdapter', adapter); + verifyEqual(testCase, report.status', {'skipped','current','current'}) + % gecko_classic asks for the BRENDA maximum: brenda.tsv is outdated. + adapter.params.kcat_assignment = struct('preset', 'gecko_classic'); + report = updateKcatSources(base, 'dryRun', true, 'modelAdapter', adapter); + verifyEqual(testCase, report.status{2}, 'outdated') +end + +function testDefaultPresetIsGeckoMedian_tc0083(testCase) + % Without [kcat_assignment] the preset is gecko_median: the BRENDA median, + % and BRENDA matches beyond origin 4 rank below DLKcat. + [base, adapter] = buildProject(); + verifyFalse(testCase, isfield(adapter.params, 'kcat_assignment')) + fuzzy = kcatsFromBrenda(base, 'aggregation', 'median', 'modelAdapter', adapter); + dlkcat = kcatsFromDLKcat(base, 'modelAdapter', adapter); + merged = mergeKcats({fuzzy, dlkcat}, {'database_top','dlkcat','database_bottom'}, 'topOriginLimit', 4); + A = assignKcatValues(base, merged); + A = applyCustomKcats(A, 'modelAdapter', adapter); + A = getKcatAcrossIsozymes(A); + A = getStandardKcat(A, 'modelAdapter', adapter); + A = applyKcatConstraints(A); + [B, ~, report] = assignKcats(base, 'modelAdapter', adapter, 'writeReport', false); + verifyEqual(testCase, report.preset, 'gecko_median') + verifyEqual(testCase, B.ec.kcat, A.ec.kcat) + verifyEqual(testCase, B.ec.source, A.ec.source) + verifyEqual(testCase, B.S, A.S) +end + +function [base, adapter, project] = buildProject() + geckoPath = findGECKOroot; + adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); + project = tempname; + copyfile(fullfile(adapter.params.path, 'data'), fullfile(project, 'data')); + adapter.params.path = project; + model = getGeckoTestModel(); + base = makeEcModel(model, false, adapter); + base = getECfromGEM(base); + % DLKcat output for this model. + rows = {'R2_EXP_1','G1','m1','C','MRAL','12.5'; 'R2_EXP_1','G2','m1','C','MNTD','6.5'; ... + 'R3','G4','m1','C','MDFM','3.25'}; + fid = fopen(fullfile(project, 'data', 'DLKcat.tsv'), 'w'); + fprintf(fid, 'rxns\tgenes\tsubstrates\tsmiles\tsequence\tkcat\n'); + for i = 1:size(rows, 1) + fprintf(fid, '%s\t%s\t%s\t%s\t%s\t%s\n', rows{i,:}); + end + fclose(fid); +end + +function writeLines(file, lines) +folder = fileparts(file); +if ~isfolder(folder) + mkdir(folder); +end +fid = fopen(file, 'w'); +fprintf(fid, '%s\n', lines{:}); +fclose(fid); +end diff --git a/tutorials/full_ecModel/model_adapter.toml b/tutorials/full_ecModel/model_adapter.toml index ea7d466a..b22add3a 100644 --- a/tutorials/full_ecModel/model_adapter.toml +++ b/tutorials/full_ecModel/model_adapter.toml @@ -79,6 +79,12 @@ dlkcat = 0.4 brenda = 0.2 custom = 0.1 +# Kcats are assigned as in the GECKO 3 tutorial: custom kcats, BRENDA +# without wildcards, DLKcat, then the other BRENDA matches, using the BRENDA +# maximum (see collectKcats). +[kcat_assignment] +preset = "gecko_classic" + # Read by GECKO only. During kcat tuning, anaerobic conditions are set by # anaerobicModel_GECKO in the code folder, which GECKO adds to the path. [matlab] diff --git a/tutorials/full_ecModel/protocol.m b/tutorials/full_ecModel/protocol.m index 84395d3a..0165004e 100644 --- a/tutorials/full_ecModel/protocol.m +++ b/tutorials/full_ecModel/protocol.m @@ -128,7 +128,11 @@ % STEP 15 Different kcat sources % Decide which kcat source to use. In the steps below, all options are % shown. Not all are required, it is up to the user to decide which ones -% they want to use. +% they want to use. Each source is kept as a file in data/kcats (brenda.tsv, +% dlkcat.tsv, okp_.tsv, custom.tsv), which assignKcats combines in +% STEP 26-31. [kcat_assignment] in model_adapter.toml sets how; this +% tutorial uses the gecko_classic preset. updateKcatSources(ecModel) +% gathers all sources that the configuration needs in one call. % STEP 16-17 Gather EC numbers, required for BRENDA as kcat source % Requires EC numbers, which are here first taken from the starting model, @@ -142,9 +146,9 @@ ecModel = getECfromDatabase(ecModel); % STEP 18-19 Gather kcat values from BRENDA -kcatList_fuzzy = fuzzyKcatMatching(ecModel); -% Now we have a kcatList, which will be used to update ecModel in a later -% step. +% Fuzzy matching against the BRENDA data (see fuzzyKcatMatching), kept as +% data/kcats/brenda.tsv. The kcatList is returned as well. +kcatList_fuzzy = kcatsFromBrenda(ecModel); % STEP 20-22 Gather metabolite SMILES, required for DLKcat as kcat source % Metabolite SMILES are gathered from PubChem. @@ -165,36 +169,37 @@ %runDLKcat(); % STEP 25 Load DLKcat output -kcatList_DLKcat = readDLKcatOutput(ecModel); - -% 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 doc mergeKcats for the tier definitions. -kcatList_merged = mergeKcats({kcatList_fuzzy, kcatList_DLKcat}, ... - {'database_top', 'dlkcat', 'database_bottom'}); - -% STEP 27 Take kcatList and populate edModel.ec.kcat -ecModel = assignKcatValues(ecModel, kcatList_merged); - -% STEP 28 Apply custom kcat values -% During the development of yeast-GEM ecModels (through GECKO 1 & 2), -% various kcat values have been manually curated, to result in a model that -% is able to validate a wide range of phenotypes. These curations are -% summarized under /data/customKcats.tsv, and applied here. -[ecModel, rxnUpdated, notMatch] = applyCustomKcats(ecModel); - -% To modify the S-matrix, to actually implement the kcat/MW constraints, -% applyKcatConstraints should be run. This is done in STEP 31. - -% STEP 29 Get kcat values across isozymes -ecModel = getKcatAcrossIsozymes(ecModel); - -% STEP 30 Get standard kcat -% Assign a protein cost to reactions without gene assocation. These -% reactions are identified as those without a corresponding entry in -% ecModel.grRules. The following reactions are exempted: +% The DLKcat output in data/DLKcat.tsv, kept as data/kcats/dlkcat.tsv. +kcatList_DLKcat = kcatsFromDLKcat(ecModel); + +% STEP 26-31 Assign kcat values +% assignKcats combines the kcat source files and applies the result to +% ecModel. With the gecko_classic preset of this tutorial: +% 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 6 ('database_top'), then DLKcat, then +% any remaining BRENDA match ('database_bottom'). See doc mergeKcats for +% the tier definitions. +% STEP 27 The kcats are written to ecModel.ec.kcat. +% STEP 28 Custom kcat values. During the development of yeast-GEM ecModels +% (through GECKO 1 & 2), various kcat values have been manually curated, +% to result in a model that is able to validate a wide range of +% phenotypes. These curations are summarized under data/customKcats.tsv. +% A data/kcats/custom.tsv file (see saveKcatSource) is used instead when +% it exists. +% STEP 29 Reactions without a kcat take one from their isozymes. +% STEP 30 Get standard kcat. Assign a protein cost to reactions without +% gene assocation, see below. +% STEP 31 Apply kcat constraints from ecModel.ec.kcat to ecModel.S, see +% below. +% The decision per reaction is in output/_kcatDecisions.tsv: +% which stage and source gave its kcat, and how many values there were. +% collectKcats gives the same decisions without changing ecModel. +[ecModel, kcatDecisions] = assignKcats(ecModel); + +% About STEP 30: the standard kcat (getStandardKcat) is assigned to +% reactions without gene assocation. These reactions are identified as +% those without a corresponding entry in ecModel.grRules. The following reactions are exempted: % A Exchange reactions: exchanging a metabolite across the model boundary, % not representing a real enzymatic reaction. % B Transport reactions: transporting a metabolite with the same name from @@ -210,20 +215,17 @@ % correctly identify all non-enzyme reactions that should be ignored by % getStandardKcat, /data/pseudoRxns.tsv can be specified in the adapter % folder, containing the relevant reaction identifiers. -[ecModel, rxnsMissingGPR, standardMW, standardKcat] = getStandardKcat(ecModel); -% STEP 31 Apply kcat constraints from ecModel.ec.kcat to ecModel.S -% While STEP 27-30 add or modify values in ecModel.ec.kcat, these -% contraints are not yet applied to the S-matrix: the enzyme is not yet set -% as pseudo- substrate with the -MW/kcat stoichiometry. Now, -% applyKcatConstraints will take the values from ecModel.ec.kcat, considers +% About STEP 31: while STEP 27-30 add or modify values in ecModel.ec.kcat, +% these contraints are not yet applied to the S-matrix: the enzyme is not +% yet set as pseudo- substrate with the -MW/kcat stoichiometry. Then, +% applyKcatConstraints takes the values from ecModel.ec.kcat, considers % any complex/subunit data that is tracked in ecModel.ec.rxnEnzMat, % together with the MW in ecModel.ec.mw, and uses this to modify the enzyme % usage coefficients directly in ecModel.S. Any time a change is made to % the .kcat, .rxnEnzMat or .mw fields, the applyKcatConstraints function % should be run again to reapply the new constraints onto the metabolic -% model. -ecModel = applyKcatConstraints(ecModel); +% model. assignKcats ran it already. % STEP 32 Set upper bound of protein pool % The protein pool exchange is constrained by the total protein content @@ -337,11 +339,11 @@ % As example, the kcat of 5'-phosphoribosylformyl glycinamidine synthetase % (reaction r_0079) was increased from 0.05 to 0.5 by sensitivityTuning. % Inspecting the kcat source might help to determine if this is reasonable. -rxnIdx = find(strcmp(kcatList_merged.rxns,'r_0079')); +rxnIdx = find(strcmp(kcatList_fuzzy.rxns,'r_0079')); doc fuzzyKcatMatching % To check the meaning of wildcardLvl and origin. -kcatList_merged.wildcardLvl(rxnIdx) % 0: no EC number wildcard. -kcatList_merged.origin(rxnIdx) % 4: any organism, any substrate, kcat. -kcatList_merged.eccodes(rxnIdx) % EC number 6.3.5.3. +kcatList_fuzzy.wildcardLvl(rxnIdx) % 0: no EC number wildcard. +kcatList_fuzzy.origin(rxnIdx) % 4: any organism, any substrate, kcat. +kcatList_fuzzy.eccodes(rxnIdx) % EC number 6.3.5.3. % On BRENDA https://www.brenda-enzymes.org/enzyme.php?ecno=6.3.5.3#TURNOVER%20NUMBER%20[1/s] % the kcat value is from E. coli with NH4+ as substrate. The reaction From 4f06b59bb4fbdd47a0a38e5d15c25ac5b73c6566 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 22:10:47 +0200 Subject: [PATCH 3/3] Read kcat source files in loadKcatSource; no private helper folder loadKcatSource reads the file itself and returns the header lines and the rows as read in report.header and report.records; readKcatSource is gone. The helpers of src/gather_kcats/private are local functions of the functions that use them. kcatsFrom no longer check the file name before gathering; saveKcatSource does when it writes. --- src/gather_kcats/kcatsFromBrenda.m | 35 ++++- src/gather_kcats/kcatsFromDLKcat.m | 35 ++++- src/gather_kcats/kcatsFromOKP.m | 35 ++++- src/gather_kcats/loadKcatSource.m | 162 +++++++++++++++++++- src/gather_kcats/private/checkNameValue.m | 21 --- src/gather_kcats/private/countRxns.m | 8 - src/gather_kcats/private/formatKcatNumber.m | 24 --- src/gather_kcats/private/kcatSourceFile.m | 58 ------- src/gather_kcats/private/kcatSourceJson.m | 61 -------- src/gather_kcats/private/parseHeaderLine.m | 14 -- src/gather_kcats/private/parseKcatNumber.m | 11 -- src/gather_kcats/private/unquoteField.m | 11 -- src/gather_kcats/private/warnSubset.m | 7 - src/gather_kcats/readKcatSource.m | 137 ----------------- src/gather_kcats/saveKcatSource.m | 161 ++++++++++++++++++- test/VerificationMatrix.tsv | 1 - test/unit_tests/kcatSourceTests.m | 27 ++-- 17 files changed, 433 insertions(+), 375 deletions(-) delete mode 100644 src/gather_kcats/private/checkNameValue.m delete mode 100644 src/gather_kcats/private/countRxns.m delete mode 100644 src/gather_kcats/private/formatKcatNumber.m delete mode 100644 src/gather_kcats/private/kcatSourceFile.m delete mode 100644 src/gather_kcats/private/kcatSourceJson.m delete mode 100644 src/gather_kcats/private/parseHeaderLine.m delete mode 100644 src/gather_kcats/private/parseKcatNumber.m delete mode 100644 src/gather_kcats/private/unquoteField.m delete mode 100644 src/gather_kcats/private/warnSubset.m delete mode 100644 src/gather_kcats/readKcatSource.m diff --git a/src/gather_kcats/kcatsFromBrenda.m b/src/gather_kcats/kcatsFromBrenda.m index 39721861..64d29428 100644 --- a/src/gather_kcats/kcatsFromBrenda.m +++ b/src/gather_kcats/kcatsFromBrenda.m @@ -44,7 +44,7 @@ % -------- % fuzzyKcatMatching, kcatsFromDLKcat, kcatsFromOKP, saveKcatSource, loadKcatSource -checkNameValue(varargin, {'aggregation','forceWClvl','ecRxns','save','saveAs','modelAdapter'}, 'kcatsFromBrenda'); +checkOptions(varargin, {'aggregation','forceWClvl','ecRxns','save','saveAs','modelAdapter'}, 'kcatsFromBrenda'); p = parseGECKOargs(varargin, { ... 'aggregation', []; ... 'forceWClvl', 0; ... @@ -61,7 +61,7 @@ end aggregation = resolveKcatAggregation(p.aggregation, modelAdapter); if p.save - file = kcatSourceFile(p.saveAs, [], modelAdapter); % fails before the matching, not after + file = fullfile(modelAdapter.params.path, 'data', 'kcats', [char(p.saveAs) '.tsv']); warnSubset(p.ecRxns, numel(model.ec.rxns), file); end @@ -95,3 +95,34 @@ 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 diff --git a/src/gather_kcats/kcatsFromDLKcat.m b/src/gather_kcats/kcatsFromDLKcat.m index fb50a39e..4608e3d2 100644 --- a/src/gather_kcats/kcatsFromDLKcat.m +++ b/src/gather_kcats/kcatsFromDLKcat.m @@ -44,7 +44,7 @@ % -------- % writeDLKcatInput, runDLKcat, readDLKcatOutput, kcatsFromBrenda, saveKcatSource -checkNameValue(varargin, {'run','ecRxns','save','saveAs','modelAdapter'}, 'kcatsFromDLKcat'); +checkOptions(varargin, {'run','ecRxns','save','saveAs','modelAdapter'}, 'kcatsFromDLKcat'); p = parseGECKOargs(varargin, { ... 'run', false; ... 'ecRxns', []; ... @@ -63,7 +63,7 @@ 'it needs ''run'', true.']) end if p.save - file = kcatSourceFile(p.saveAs, [], modelAdapter); + file = fullfile(modelAdapter.params.path, 'data', 'kcats', [char(p.saveAs) '.tsv']); if p.run warnSubset(p.ecRxns, numel(model.ec.rxns), file); end @@ -112,3 +112,34 @@ 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 diff --git a/src/gather_kcats/kcatsFromOKP.m b/src/gather_kcats/kcatsFromOKP.m index 4e0ed5d1..5fe891d4 100644 --- a/src/gather_kcats/kcatsFromOKP.m +++ b/src/gather_kcats/kcatsFromOKP.m @@ -67,7 +67,7 @@ % -------- % submitOpenKineticsPredictor, fetchOpenKineticsPredictor, kcatsFromBrenda, saveKcatSource -checkNameValue(varargin, {'method','wait','resubmit','useStored','ecRxns','save','saveAs','modelAdapter'}, 'kcatsFromOKP'); +checkOptions(varargin, {'method','wait','resubmit','useStored','ecRxns','save','saveAs','modelAdapter'}, 'kcatsFromOKP'); p = parseGECKOargs(varargin, { ... 'method', []; ... 'wait', false; ... @@ -108,7 +108,7 @@ 'submits none (a job for %s exists; use ''resubmit'', true for a new one).'], method) end if p.save - file = kcatSourceFile(name, [], modelAdapter); % fails before submitting, not after + file = fullfile(modelAdapter.params.path, 'data', 'kcats', [char(name) '.tsv']); if submitting warnSubset(p.ecRxns, numel(model.ec.rxns), file); end @@ -212,3 +212,34 @@ function writeJobFile(file, job) end fclose(fid); 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 diff --git a/src/gather_kcats/loadKcatSource.m b/src/gather_kcats/loadKcatSource.m index 87d99867..56e2d950 100644 --- a/src/gather_kcats/loadKcatSource.m +++ b/src/gather_kcats/loadKcatSource.m @@ -19,6 +19,14 @@ % % The enzymes of a row may be listed in any order. % +% 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 reported as invalid. A row with more filled fields than there are +% columns is an error. +% % Parameters % ---------- % model : struct @@ -44,7 +52,12 @@ % - invalid: table of the file rows with a missing or invalid kcat; % - unmatched: table of the file rows that apply to no model.ec row; % - partial: table of rows that only partly match an enzyme set, with -% the model.ec reaction and its enzymes. +% the model.ec reaction and its enzymes; +% - header: the '#' header lines of the file, as field: value (for +% instance header.source, header.parameters); +% - records: the rows of the file as read, one field per column (kcat +% and log10_sd as double, other columns as cell of char) and line, +% the line number of each row. % % Examples % -------- @@ -53,9 +66,9 @@ % % See also % -------- -% saveKcatSource, readKcatSource, kcatsFromBrenda, getEcRxnKeys +% saveKcatSource, kcatsFromBrenda, getEcRxnKeys -checkNameValue(varargin, {'modelAdapter'}, 'loadKcatSource'); +checkOptions(varargin, {'modelAdapter'}, 'loadKcatSource'); p = parseGECKOargs(varargin, {'modelAdapter', []}); if isfile(name) @@ -75,7 +88,7 @@ 'kcatsFromOKP), or with saveKcatSource.'], filename) end end -[records, meta] = readKcatSource(filename); +[records, meta] = readSourceFile(filename); sourceName = ''; if isfield(meta, 'source') sourceName = meta.source; @@ -164,6 +177,8 @@ printOrange(sprintf('%d of %d rows of %s apply to no reaction of this model; see report.unmatched.\n', ... sum(unmatched), nRec, filename)); end +report.header = meta; +report.records = records; end function map = groupIndex(keys, rows) @@ -216,3 +231,142 @@ end end end + +function [records, meta] = readSourceFile(filename) +% The rows and '#' header lines of a kcat source file, as they are. +text = readUtf8(filename); +lines = regexp(text, '\r?\n', 'split'); + +meta = struct(); +k = 1; +while k <= numel(lines) && (isempty(strtrim(lines{k})) || startsWith(unquoteField(lines{k}), '#')) + [key, value] = parseHeaderLine(lines{k}); + if ~isempty(key) + meta.(key) = value; + end + k = k + 1; +end +if k > numel(lines) + error('loadKcatSource:noHeader', '%s has no column header.', filename) +end +columns = cellfun(@unquoteField, strsplit(lines{k}, sprintf('\t'), ... + 'CollapseDelimiters', false), 'UniformOutput', false); +if ~all(ismember({'kcat','source'}, columns)) || ~any(ismember({'rxn','enzymes'}, columns)) + error('loadKcatSource:columns', ['%s needs the columns kcat and source, and rxn ' ... + 'and/or enzymes.'], filename) +end +named = columns(~cellfun(@isempty, columns)); +if numel(unique(named)) < numel(named) || any(strcmp(columns, 'line')) + error('loadKcatSource:columns', '%s has a column name twice, or a column named line.', filename) +end +nCol = numel(columns); + +trimmed = strtrim(lines); +skip = cellfun(@isempty, trimmed) | startsWith(trimmed, '#') | startsWith(trimmed, '"#'); +skip(1:k) = true; +kept = find(~skip); +L = string(lines(kept)); +L = L(:); +F = strings(numel(kept), nCol); +nTabs = count(L, sprintf('\t')); +regular = nTabs == nCol - 1; +if any(regular) + % Rows with exactly one field per column, split in one go. + F(regular, :) = reshape(split(L(regular), sprintf('\t'), 2), [], nCol); +end +for r = find(~regular)' + fields = strsplit(char(L(r)), sprintf('\t'), 'CollapseDelimiters', false); + if numel(fields) > nCol + if any(~cellfun(@isempty, strtrim(fields(nCol+1:end)))) + error('loadKcatSource:fields', '%s, line %d has more fields than the %d columns.', ... + filename, kept(r), nCol) + end + fields = fields(1:nCol); + end + fields(end+1:nCol) = {''}; + F(r, :) = string(fields); +end +F = strtrim(F); +quoted = find(startsWith(F, '"')); +for q = quoted(:)' + F(q) = unquoteField(char(F(q))); +end +filled = any(F ~= "", 2); % a row of tabs only is skipped +values = cellstr(F(filled, :)); +if isempty(values) + values = cell(0, nCol); +end +lineNo = kept(filled); +lineNo = lineNo(:); +n = numel(lineNo); + +records = struct(); +for c = 1:nCol + if isempty(columns{c}) + continue + end + name = matlab.lang.makeValidName(columns{c}); + v = values(:, c); + if any(strcmp(columns{c}, {'kcat','log10_sd'})) + v = parseKcatNumber(v); + end + records.(name) = v; +end +for c = {'rxn','enzymes','substrate','source','detail','ec_code','note'} + if ~isfield(records, c{1}) + records.(c{1}) = repmat({''}, n, 1); + end +end +if ~isfield(records, 'log10_sd') + records.log10_sd = nan(n, 1); +end +records.line = lineNo; +end + +function text = readUtf8(filename) +fid = fopen(filename, 'r', 'n', 'UTF-8'); +if fid < 0 + error('loadKcatSource:read', 'Cannot read %s', filename) +end +cleanup = onCleanup(@() fclose(fid)); +text = fread(fid, [1 Inf], '*char'); +text = regexprep(text, ['^' char(65279)], ''); +end + +function x = parseKcatNumber(texts) +% Plain numbers only: '5,5' or '12 /s' read as NaN. +x = nan(numel(texts), 1); +ok = ~cellfun(@isempty, regexp(texts(:), '^[-+]?(\d+\.?\d*|\.\d+)([eE][-+]?\d+)?$|^[-+]?Inf$', 'once')); +x(ok) = str2double(texts(ok)); +end + +function [key, value] = parseHeaderLine(line) +key = ''; value = ''; +tok = regexp(unquoteField(line), '^\s*#\s*([A-Za-z][A-Za-z0-9_]{0,62})\s*:\s?(.*)$', 'tokens', 'once'); +if ~isempty(tok) + key = lower(tok{1}); + value = strtrim(tok{2}); +end +end + +function t = unquoteField(t) +t = strtrim(t); +if numel(t) >= 2 && t(1) == '"' && t(end) == '"' + t = strtrim(strrep(t(2:end-1), '""', '"')); +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 diff --git a/src/gather_kcats/private/checkNameValue.m b/src/gather_kcats/private/checkNameValue.m deleted file mode 100644 index ed2e7bb0..00000000 --- a/src/gather_kcats/private/checkNameValue.m +++ /dev/null @@ -1,21 +0,0 @@ -function checkNameValue(args, names, functionName) -% checkNameValue Accept only name-value pairs with known names. -% -% The kcat source functions have no positional form, so a misspelled -% option name must be an error rather than being read as a positional -% argument by parseGECKOargs. -id = [functionName ':arguments']; -if mod(numel(args), 2) ~= 0 - error(id, 'Options of %s are given as name-value pairs.', functionName) -end -for k = 1:2:numel(args) - name = args{k}; - if ~(ischar(name) || (isstring(name) && isscalar(name))) - error(id, 'Option names of %s are text, not %s. Options: %s.', functionName, ... - class(name), strjoin(names, ', ')) - elseif ~any(strcmp(char(name), names)) - error(id, 'Unknown option of %s: %s. Options: %s.', functionName, char(name), ... - strjoin(names, ', ')) - end -end -end diff --git a/src/gather_kcats/private/countRxns.m b/src/gather_kcats/private/countRxns.m deleted file mode 100644 index fd812fab..00000000 --- a/src/gather_kcats/private/countRxns.m +++ /dev/null @@ -1,8 +0,0 @@ -function n = countRxns(ecRxns) -% countRxns Number of reactions an ecRxns argument selects (logical or indices). -if islogical(ecRxns) - n = sum(ecRxns); -else - n = numel(unique(ecRxns)); -end -end diff --git a/src/gather_kcats/private/formatKcatNumber.m b/src/gather_kcats/private/formatKcatNumber.m deleted file mode 100644 index f7251a3e..00000000 --- a/src/gather_kcats/private/formatKcatNumber.m +++ /dev/null @@ -1,24 +0,0 @@ -function t = formatKcatNumber(x) -% formatKcatNumber Number text of kcat source files. -% -% Whole numbers as integers, others in the shortest decimal form that reads -% back as the same double; 'Inf' or '-Inf' for infinity, empty for NaN. -% geckopy's format_number writes the same text. -if isnan(x) - t = ''; -elseif isinf(x) - t = 'Inf'; - if x < 0 - t = '-Inf'; - end -elseif x == round(x) && abs(x) < 1e15 - t = sprintf('%d', x); -else - for digits = 1:17 - t = sprintf('%.*g', digits, x); - if str2double(t) == x - break - end - end -end -end diff --git a/src/gather_kcats/private/kcatSourceFile.m b/src/gather_kcats/private/kcatSourceFile.m deleted file mode 100644 index 3d112cf3..00000000 --- a/src/gather_kcats/private/kcatSourceFile.m +++ /dev/null @@ -1,58 +0,0 @@ -function filename = kcatSourceFile(name, folder, modelAdapter) -% kcatSourceFile Path of a kcat source file that may be written. -% -% Used by saveKcatSource and the kcatsFrom functions before any -% work is done: checks the name, resolves the folder (data/kcats in the -% project folder by default), and refuses a file that is curated by hand: -% named custom.tsv, or whose header says 'kind: curated' or 'source: -% custom' (in any case, also when a spreadsheet quoted the line). - -name = char(name); -if isempty(regexp(name, '^[A-Za-z0-9_.-]+$', 'once')) - error('saveKcatSource:invalidName', ... - 'The source name "%s" may only contain letters, digits, ''_'', ''-'' and ''.''.', name) -end -if isempty(folder) - 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 - folder = fullfile(modelAdapter.params.path, 'data', 'kcats'); -end -filename = fullfile(char(folder), [name '.tsv']); -if strcmpi(name, 'custom') - error('saveKcatSource:curated', ['custom.tsv holds curated kcats and is not ' ... - 'written by GECKO; edit it by hand.']) -end -if isfile(filename) - meta = readHeaderLines(filename); - if (isfield(meta, 'kind') && strcmpi(meta.kind, 'curated')) || ... - (isfield(meta, 'source') && strcmpi(meta.source, 'custom')) - error('saveKcatSource:curated', '%s is curated by hand and is not overwritten.', filename) - end -end -end - -function meta = readHeaderLines(filename) -% Only the header lines at the top, so that a damaged table does not stop -% a gathered file from being written again. -meta = struct(); -fid = fopen(filename, 'r', 'n', 'UTF-8'); -if fid < 0 - return -end -cleanup = onCleanup(@() fclose(fid)); -line = fgetl(fid); -if ischar(line) - line = regexprep(line, ['^' char(65279)], ''); -end -while ischar(line) && (isempty(strtrim(line)) || startsWith(unquoteField(line), '#')) - [key, value] = parseHeaderLine(line); - if ~isempty(key) - meta.(key) = value; - end - line = fgetl(fid); -end -end diff --git a/src/gather_kcats/private/kcatSourceJson.m b/src/gather_kcats/private/kcatSourceJson.m deleted file mode 100644 index 86808bd1..00000000 --- a/src/gather_kcats/private/kcatSourceJson.m +++ /dev/null @@ -1,61 +0,0 @@ -function j = kcatSourceJson(v) -% kcatSourceJson The parameters line of kcat source files as JSON. -% -% Compact JSON with keys sorted at every level. Numbers are written with -% formatKcatNumber, and NaN and infinity as null; logical values as true -% or false; char as a string (non-ASCII characters as they are); a cell -% array as a list; a numeric or logical vector with more than one element -% as a list. geckopy's save_kcat_source writes the same text. -if isstruct(v) && isscalar(v) - f = sort(fieldnames(v)); - parts = cell(1, numel(f)); - for k = 1:numel(f) - parts{k} = [jsonString(f{k}) ':' kcatSourceJson(v.(f{k}))]; - end - j = ['{' strjoin(parts, ',') '}']; -elseif ischar(v) || (isstring(v) && isscalar(v)) - j = jsonString(char(v)); -elseif iscell(v) - j = ['[' strjoin(cellfun(@kcatSourceJson, v(:)', 'UniformOutput', false), ',') ']']; -elseif (isnumeric(v) || islogical(v)) && numel(v) ~= 1 - j = ['[' strjoin(arrayfun(@kcatSourceJson, v(:)', 'UniformOutput', false), ',') ']']; -elseif islogical(v) - j = 'false'; - if v - j = 'true'; - end -elseif isnumeric(v) - if ~isfinite(v) - j = 'null'; - else - j = formatKcatNumber(double(v)); - end -elseif isempty(v) - j = 'null'; -else - error('kcatSourceJson:type', 'Parameters of class %s cannot be written.', class(v)) -end -end - -function s = jsonString(s) -% JSON escapes, as Python's json.dumps(ensure_ascii=False). -out = ''; -for c = s - switch c - case '"', out = [out '\"']; %#ok - case '\', out = [out '\\']; %#ok - case char(8), out = [out '\b']; %#ok - case char(9), out = [out '\t']; %#ok - case char(10), out = [out '\n']; %#ok - case char(12), out = [out '\f']; %#ok - case char(13), out = [out '\r']; %#ok - otherwise - if double(c) < 32 - out = [out sprintf('\\u%04x', double(c))]; %#ok - else - out = [out c]; %#ok - end - end -end -s = ['"' out '"']; -end diff --git a/src/gather_kcats/private/parseHeaderLine.m b/src/gather_kcats/private/parseHeaderLine.m deleted file mode 100644 index ee17bbc4..00000000 --- a/src/gather_kcats/private/parseHeaderLine.m +++ /dev/null @@ -1,14 +0,0 @@ -function [key, value] = parseHeaderLine(line) -% parseHeaderLine Key and value of a '# key: value' header line. -% -% The line may be wrapped in quotes by a spreadsheet. The key is returned -% in lower case; it must start with a letter and be at most 63 letters, -% digits or underscores long. Other lines give an empty key. geckopy reads -% header lines the same way. -key = ''; value = ''; -tok = regexp(unquoteField(line), '^\s*#\s*([A-Za-z][A-Za-z0-9_]{0,62})\s*:\s?(.*)$', 'tokens', 'once'); -if ~isempty(tok) - key = lower(tok{1}); - value = strtrim(tok{2}); -end -end diff --git a/src/gather_kcats/private/parseKcatNumber.m b/src/gather_kcats/private/parseKcatNumber.m deleted file mode 100644 index f17b49a0..00000000 --- a/src/gather_kcats/private/parseKcatNumber.m +++ /dev/null @@ -1,11 +0,0 @@ -function x = parseKcatNumber(texts) -% parseKcatNumber Read numbers written in kcat source files. -% -% Accepts plain decimal numbers with an optional exponent (12, 0.5, 1e-05, -% 2.5E+3) and Inf; anything else, such as '5,5' from a spreadsheet with a -% decimal comma, becomes NaN instead of being misread. geckopy's -% read_kcat_source accepts the same text. -x = nan(numel(texts), 1); -ok = ~cellfun(@isempty, regexp(texts(:), '^[-+]?(\d+\.?\d*|\.\d+)([eE][-+]?\d+)?$|^[-+]?Inf$', 'once')); -x(ok) = str2double(texts(ok)); -end diff --git a/src/gather_kcats/private/unquoteField.m b/src/gather_kcats/private/unquoteField.m deleted file mode 100644 index 36c95aa1..00000000 --- a/src/gather_kcats/private/unquoteField.m +++ /dev/null @@ -1,11 +0,0 @@ -function t = unquoteField(t) -% unquoteField 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; spaces inside the quotes are -% trimmed too. -t = strtrim(t); -if numel(t) >= 2 && t(1) == '"' && t(end) == '"' - t = strtrim(strrep(t(2:end-1), '""', '"')); -end -end diff --git a/src/gather_kcats/private/warnSubset.m b/src/gather_kcats/private/warnSubset.m deleted file mode 100644 index 00064a5b..00000000 --- a/src/gather_kcats/private/warnSubset.m +++ /dev/null @@ -1,7 +0,0 @@ -function warnSubset(ecRxns, nRxns, file) -% warnSubset Say that gathering for some reactions replaces the whole 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 diff --git a/src/gather_kcats/readKcatSource.m b/src/gather_kcats/readKcatSource.m deleted file mode 100644 index 063ab8b0..00000000 --- a/src/gather_kcats/readKcatSource.m +++ /dev/null @@ -1,137 +0,0 @@ -function [records, meta] = readKcatSource(filename) -% readKcatSource Read a kcat source file as it is, without a model. -% -% Reads a file written by saveKcatSource or a kcatsFrom function, -% or curated by hand in the same format. Use loadKcatSource to map the -% rows to an ecModel. -% -% The file is read as UTF-8 (a leading byte-order mark is ignored). The -% '#' lines before the column header are the header (keys in lower case); -% '#' 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 loadKcatSource reports. A row with more filled -% fields than there are columns is an error. -% -% Parameters -% ---------- -% filename : char -% path to the kcat source file. -% -% Returns -% ------- -% records : struct -% one field per column, each with one entry per data row: kcat and -% log10_sd as double (NaN when empty or not a number), every other -% column as a cell of char. Columns other than the standard ones are -% kept. records.line is the line number of each row in the file. -% meta : struct -% the '#' header lines, as field: value (for instance meta.source, -% meta.kind). -% -% Examples -% -------- -% [records, meta] = readKcatSource(fullfile(projectPath,'data','kcats','brenda.tsv')); -% -% See also -% -------- -% loadKcatSource, saveKcatSource - -text = readUtf8(filename); -lines = regexp(text, '\r?\n', 'split'); - -meta = struct(); -k = 1; -while k <= numel(lines) && (isempty(strtrim(lines{k})) || startsWith(unquoteField(lines{k}), '#')) - [key, value] = parseHeaderLine(lines{k}); - if ~isempty(key) - meta.(key) = value; - end - k = k + 1; -end -if k > numel(lines) - error('readKcatSource:noHeader', '%s has no column header.', filename) -end -columns = cellfun(@unquoteField, strsplit(lines{k}, sprintf('\t'), ... - 'CollapseDelimiters', false), 'UniformOutput', false); -if ~all(ismember({'kcat','source'}, columns)) || ~any(ismember({'rxn','enzymes'}, columns)) - error('readKcatSource:columns', ['%s needs the columns kcat and source, and rxn ' ... - 'and/or enzymes.'], filename) -end -named = columns(~cellfun(@isempty, columns)); -if numel(unique(named)) < numel(named) || any(strcmp(columns, 'line')) - error('readKcatSource:columns', '%s has a column name twice, or a column named line.', filename) -end -nCol = numel(columns); - -trimmed = strtrim(lines); -skip = cellfun(@isempty, trimmed) | startsWith(trimmed, '#') | startsWith(trimmed, '"#'); -skip(1:k) = true; -kept = find(~skip); -L = string(lines(kept)); -L = L(:); -F = strings(numel(kept), nCol); -nTabs = count(L, sprintf('\t')); -regular = nTabs == nCol - 1; -if any(regular) - % Rows with exactly one field per column, split in one go. - F(regular, :) = reshape(split(L(regular), sprintf('\t'), 2), [], nCol); -end -for r = find(~regular)' - fields = strsplit(char(L(r)), sprintf('\t'), 'CollapseDelimiters', false); - if numel(fields) > nCol - if any(~cellfun(@isempty, strtrim(fields(nCol+1:end)))) - error('readKcatSource:fields', '%s, line %d has more fields than the %d columns.', ... - filename, kept(r), nCol) - end - fields = fields(1:nCol); - end - fields(end+1:nCol) = {''}; - F(r, :) = string(fields); -end -F = strtrim(F); -quoted = find(startsWith(F, '"')); -for q = quoted(:)' - F(q) = unquoteField(char(F(q))); -end -filled = any(F ~= "", 2); % a row of tabs only is skipped -values = cellstr(F(filled, :)); -if isempty(values) - values = cell(0, nCol); -end -lineNo = kept(filled); -lineNo = lineNo(:); -n = numel(lineNo); - -records = struct(); -for c = 1:nCol - if isempty(columns{c}) - continue - end - name = matlab.lang.makeValidName(columns{c}); - v = values(:, c); - if any(strcmp(columns{c}, {'kcat','log10_sd'})) - v = parseKcatNumber(v); - end - records.(name) = v; -end -for c = {'rxn','enzymes','substrate','source','detail','ec_code','note'} - if ~isfield(records, c{1}) - records.(c{1}) = repmat({''}, n, 1); - end -end -if ~isfield(records, 'log10_sd') - records.log10_sd = nan(n, 1); -end -records.line = lineNo; -end - -function text = readUtf8(filename) -fid = fopen(filename, 'r', 'n', 'UTF-8'); -if fid < 0 - error('readKcatSource:read', 'Cannot read %s', filename) -end -cleanup = onCleanup(@() fclose(fid)); -text = fread(fid, [1 Inf], '*char'); -text = regexprep(text, ['^' char(65279)], ''); -end diff --git a/src/gather_kcats/saveKcatSource.m b/src/gather_kcats/saveKcatSource.m index 37101553..330025c0 100644 --- a/src/gather_kcats/saveKcatSource.m +++ b/src/gather_kcats/saveKcatSource.m @@ -62,10 +62,10 @@ % % See also % -------- -% kcatsFromBrenda, loadKcatSource, readKcatSource, getEcRxnKeys +% kcatsFromBrenda, loadKcatSource, getEcRxnKeys names = {'modelAdapter','folder','generator','parameters','data','kind'}; -checkNameValue(varargin, names, 'saveKcatSource'); +checkOptions(varargin, names, 'saveKcatSource'); p = parseGECKOargs(varargin, { ... 'modelAdapter', []; ... 'folder', []; ... @@ -224,3 +224,160 @@ v = strtrim(fileread(file)); end end + +function filename = kcatSourceFile(name, folder, modelAdapter) +name = char(name); +if isempty(regexp(name, '^[A-Za-z0-9_.-]+$', 'once')) + error('saveKcatSource:invalidName', ... + 'The source name "%s" may only contain letters, digits, ''_'', ''-'' and ''.''.', name) +end +if isempty(folder) + 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 + folder = fullfile(modelAdapter.params.path, 'data', 'kcats'); +end +filename = fullfile(char(folder), [name '.tsv']); +if strcmpi(name, 'custom') + error('saveKcatSource:curated', ['custom.tsv holds curated kcats and is not ' ... + 'written by GECKO; edit it by hand.']) +end +if isfile(filename) + meta = readHeaderLines(filename); + if (isfield(meta, 'kind') && strcmpi(meta.kind, 'curated')) || ... + (isfield(meta, 'source') && strcmpi(meta.source, 'custom')) + error('saveKcatSource:curated', '%s is curated by hand and is not overwritten.', filename) + end +end +end + +function meta = readHeaderLines(filename) +% Only the header lines at the top, so that a damaged table does not stop +% a gathered file from being written again. +meta = struct(); +fid = fopen(filename, 'r', 'n', 'UTF-8'); +if fid < 0 + return +end +cleanup = onCleanup(@() fclose(fid)); +line = fgetl(fid); +if ischar(line) + line = regexprep(line, ['^' char(65279)], ''); +end +while ischar(line) && (isempty(strtrim(line)) || startsWith(unquoteField(line), '#')) + [key, value] = parseHeaderLine(line); + if ~isempty(key) + meta.(key) = value; + end + line = fgetl(fid); +end +end + +function j = kcatSourceJson(v) +if isstruct(v) && isscalar(v) + f = sort(fieldnames(v)); + parts = cell(1, numel(f)); + for k = 1:numel(f) + parts{k} = [jsonString(f{k}) ':' kcatSourceJson(v.(f{k}))]; + end + j = ['{' strjoin(parts, ',') '}']; +elseif ischar(v) || (isstring(v) && isscalar(v)) + j = jsonString(char(v)); +elseif iscell(v) + j = ['[' strjoin(cellfun(@kcatSourceJson, v(:)', 'UniformOutput', false), ',') ']']; +elseif (isnumeric(v) || islogical(v)) && numel(v) ~= 1 + j = ['[' strjoin(arrayfun(@kcatSourceJson, v(:)', 'UniformOutput', false), ',') ']']; +elseif islogical(v) + j = 'false'; + if v + j = 'true'; + end +elseif isnumeric(v) + if ~isfinite(v) + j = 'null'; + else + j = formatKcatNumber(double(v)); + end +elseif isempty(v) + j = 'null'; +else + error('kcatSourceJson:type', 'Parameters of class %s cannot be written.', class(v)) +end +end + +function s = jsonString(s) +% JSON escapes, as Python's json.dumps(ensure_ascii=False). +out = ''; +for c = s + switch c + case '"', out = [out '\"']; %#ok + case '\', out = [out '\\']; %#ok + case char(8), out = [out '\b']; %#ok + case char(9), out = [out '\t']; %#ok + case char(10), out = [out '\n']; %#ok + case char(12), out = [out '\f']; %#ok + case char(13), out = [out '\r']; %#ok + otherwise + if double(c) < 32 + out = [out sprintf('\\u%04x', double(c))]; %#ok + else + out = [out c]; %#ok + end + end +end +s = ['"' out '"']; +end + +function t = formatKcatNumber(x) +if isnan(x) + t = ''; +elseif isinf(x) + t = 'Inf'; + if x < 0 + t = '-Inf'; + end +elseif x == round(x) && abs(x) < 1e15 + t = sprintf('%d', x); +else + for digits = 1:17 + t = sprintf('%.*g', digits, x); + if str2double(t) == x + break + end + end +end +end + +function [key, value] = parseHeaderLine(line) +key = ''; value = ''; +tok = regexp(unquoteField(line), '^\s*#\s*([A-Za-z][A-Za-z0-9_]{0,62})\s*:\s?(.*)$', 'tokens', 'once'); +if ~isempty(tok) + key = lower(tok{1}); + value = strtrim(tok{2}); +end +end + +function t = unquoteField(t) +t = strtrim(t); +if numel(t) >= 2 && t(1) == '"' && t(end) == '"' + t = strtrim(strrep(t(2:end-1), '""', '"')); +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 diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 36c28481..f6df8d6b 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -17,7 +17,6 @@ gather_kcats/assignKcatValues.m tc0011 - we did not test all possible parameteri gather_kcats/writeDLKcatInput.m partly by tc0011 gather_kcats/getEcRxnKeys.m tc0070, tc0076 gather_kcats/saveKcatSource.m tc0071, tc0073, tc0076, tc0078 -gather_kcats/readKcatSource.m tc0071, tc0077 gather_kcats/loadKcatSource.m tc0071, tc0072, tc0073, tc0076, tc0077, tc0078 gather_kcats/kcatsFromBrenda.m tc0074, tc0078 gather_kcats/kcatsFromDLKcat.m tc0074, tc0078 - reads an existing DLKcat output; running DLKcat needs Docker diff --git a/test/unit_tests/kcatSourceTests.m b/test/unit_tests/kcatSourceTests.m index 4228cd74..d40829f1 100644 --- a/test/unit_tests/kcatSourceTests.m +++ b/test/unit_tests/kcatSourceTests.m @@ -28,7 +28,7 @@ function testSaveLoadRoundTrip_tc0071(testCase) folder = tempname; file = saveKcatSource(full, list, 'brenda', 'folder', folder, ... 'parameters', struct('aggregation','max'), 'data', 'BRENDA release test'); - [records, meta] = readKcatSource(file); + [records, meta] = rowsOf(full, file); verifyEqual(testCase, meta.source, 'brenda') verifyEqual(testCase, meta.kind, 'gathered') verifyEqual(testCase, meta.data, 'BRENDA release test') @@ -59,7 +59,7 @@ function testFullFileOnLightModel_tc0072(testCase) [loaded, report] = loadKcatSource(light, file); verifyEmpty(testCase, report.unmatched) [baseL, enzL] = getEcRxnKeys(light); - [records] = readKcatSource(file); + records = rowsOf(light, file); for r = 1:numel(records.kcat) hit = strcmp(baseL, records.rxn{r}) & strcmp(enzL, records.enzymes{r}); verifyEqual(testCase, sort(loaded.rxns(loaded.line == records.line(r))), sort(light.ec.rxns(hit))) @@ -106,7 +106,7 @@ function testKcatsFromBrendaAndDLKcat_tc0074(testCase) copyfile(fullfile(adapter.params.path, 'data'), fullfile(project, 'data')); adapter.params.path = project; list = kcatsFromBrenda(full, 'aggregation', 'median', 'modelAdapter', adapter); - [records, meta] = readKcatSource(fullfile(project, 'data', 'kcats', 'brenda.tsv')); + [records, meta] = rowsOf(full, fullfile(project, 'data', 'kcats', 'brenda.tsv')); verifyEqual(testCase, meta.generator(1:15), 'kcatsFromBrenda') verifyEqual(testCase, meta.parameters, '{"aggregation":"median","force_wildcard_level":0}') verifyEqual(testCase, meta.data, 'BRENDA release test-fixture') @@ -121,7 +121,7 @@ function testKcatsFromBrendaAndDLKcat_tc0074(testCase) end fclose(fid); kcatsFromDLKcat(full, 'modelAdapter', adapter); - [records, meta] = readKcatSource(fullfile(project, 'data', 'kcats', 'dlkcat.tsv')); + [records, meta] = rowsOf(full, fullfile(project, 'data', 'kcats', 'dlkcat.tsv')); verifyEqual(testCase, meta.source, 'dlkcat') verifyEqual(testCase, sort(records.kcat), [3.25; 12.5]) verifyEqual(testCase, unique(records.source), {'DLKcat'}) @@ -144,7 +144,7 @@ function testKcatsFromOKPStored_tc0075(testCase) writeLines(fullfile(data, 'OKP_job.txt'), {'jobId: abc123', 'method: CataPro'}); [list, done] = kcatsFromOKP(full, 'method', "CataPro", 'useStored', true, 'modelAdapter', adapter); verifyTrue(testCase, done) - [records, meta] = readKcatSource(fullfile(data, 'kcats', 'okp_catapro.tsv')); + [records, meta] = rowsOf(full, fullfile(data, 'kcats', 'okp_catapro.tsv')); verifyEqual(testCase, meta.data, 'OKP job abc123') verifyEqual(testCase, numel(records.kcat), sum(list.kcats > 0)) verifyEqual(testCase, unique(records.source), {'OKP-CataPro'}) @@ -193,7 +193,7 @@ function testGuardsNumbersAndJson_tc0078(testCase) file = saveKcatSource(full, list, "numbers", 'modelAdapter', adapter, 'parameters', ... struct('z', 1, 'a', struct('y', ['caf' char(233)], 'b', true), 'n', NaN, ... 'f', 1/3, 'g', 1e-7, 'h', 1e15, 'k', 1e20, 'l', {{5}}, 's', 'tab quote"')); - [records, meta] = readKcatSource(file); + [records, meta] = rowsOf(full, file); verifyTrue(testCase, all(isfinite(records.kcat) & records.kcat > 0)) % The same text as geckopy's test_parameters_json_matches_gecko. verifyEqual(testCase, meta.parameters, ['{"a":{"b":true,"y":"caf' char(233) '"},"f":0.3333333333333333,' ... @@ -208,7 +208,7 @@ function testGuardsNumbersAndJson_tc0078(testCase) 'R3 P4 m1 2 b ', 'R3 P4 m1 2 a wc=1', 'R3 P4 5 a ', 'R2_REV P3 m2 1 a '}); l2 = loadKcatSource(full, fullfile(folder, 'unsorted.tsv')); file = saveKcatSource(full, l2, 'sorted', 'modelAdapter', adapter); - records = readKcatSource(file); + records = rowsOf(full, file); verifyEqual(testCase, records.rxn', {'R2_REV','R3','R3','R3'}) verifyEqual(testCase, records.substrate', {'m2','','m1','m1'}) verifyEqual(testCase, records.source', {'a','a','a','b'}) @@ -238,7 +238,8 @@ function testSharedFixture_tc0077(testCase) % a quoted header line and field, comments and blank lines between % rows, a row of tabs, a decimal comma, a negative kcat. file = fullfile(findGECKOroot, 'test', 'unit_tests', 'ecTestGEM', 'data', 'kcatSourceFixture.tsv'); - [records, meta] = readKcatSource(file); + [full] = buildKcatModels(); + [records, meta] = rowsOf(full, file); verifyEqual(testCase, meta.source, 'fixture') verifyEqual(testCase, meta.parameters, '{"aggregation":"median"}') verifyEqual(testCase, meta.kind, 'curated') @@ -248,7 +249,6 @@ function testSharedFixture_tc0077(testCase) verifyEqual(testCase, records.kcat', [12.5 1e-05 NaN 7 3 -1 4]) verifyEqual(testCase, records.log10_sd(1), 0.3) - [full] = buildKcatModels(); [loaded, report] = loadKcatSource(full, file); verifyEqual(testCase, loaded.rxns(loaded.line == 7), {'R2_EXP_1'}) verifyEqual(testCase, sort(loaded.rxns(loaded.line == 12)), {'R2_REV_EXP_1';'R2_REV_EXP_2'}) @@ -259,7 +259,14 @@ function testSharedFixture_tc0077(testCase) % A row with more filled fields than columns is an error. bad = [tempname '.tsv']; writeLines(bad, {'rxn kcat source', 'R2 1 brenda extra'}); - verifyError(testCase, @() readKcatSource(bad), 'readKcatSource:fields') + verifyError(testCase, @() loadKcatSource(full, bad), 'loadKcatSource:fields') +end + +function [records, header] = rowsOf(model, file) + % The rows and header of a kcat source file, as loadKcatSource reads them. + [~, report] = loadKcatSource(model, file); + records = report.records; + header = report.header; end function writeLines(file, lines)