diff --git a/src/gather_kcats/assignKcats.m b/src/gather_kcats/assignKcats.m new file mode 100644 index 00000000..25554ef3 --- /dev/null +++ b/src/gather_kcats/assignKcats.m @@ -0,0 +1,57 @@ +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 + +checkOptions(varargin, {'modelAdapter','writeReport'}, 'assignKcats'); +[decisions, report, model] = collectKcats(model, varargin{:}); +model = applyKcatConstraints(model); +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/collectKcats.m b/src/gather_kcats/collectKcats.m new file mode 100644 index 00000000..5a46b89f --- /dev/null +++ b/src/gather_kcats/collectKcats.m @@ -0,0 +1,251 @@ +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 + +checkOptions(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 + +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/kcatAssignmentConfig.m b/src/gather_kcats/kcatAssignmentConfig.m new file mode 100644 index 00000000..36d22ec7 --- /dev/null +++ b/src/gather_kcats/kcatAssignmentConfig.m @@ -0,0 +1,109 @@ +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, as collectKcats and updateKcatSources +% use them. Without a preset or stages the preset is gecko_median. +% geckopy's kcat_assignment_config gives the same settings. +% +% Parameters +% ---------- +% params : struct +% the parameters of a model adapter (modelAdapter.params). +% +% Returns +% ------- +% config : struct +% - preset: the preset name, or '' when stages are given; +% - isozymes, standard: logical; +% - stages: struct array with action, sources (cell), criteria, +% topOriginLimit, bottomOriginLimit and wildcardLimit; +% - sources: per source name its options (from [kcat_sources] and the +% preset). +% +% Examples +% -------- +% config = kcatAssignmentConfig(modelAdapter.params); +% config.stages(1).sources % the order of the sources +% +% See also +% -------- +% collectKcats, assignKcats, updateKcatSources + +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..ed8523f6 --- /dev/null +++ b/src/gather_kcats/updateKcatSources.m @@ -0,0 +1,259 @@ +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 + +checkOptions(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 = readHeaderLines(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 + +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 [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/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 f6df8d6b..55de64b3 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -21,6 +21,10 @@ 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/kcatAssignmentConfig.m tc0079, tc0081, tc0083 gather_kcats/readOpenKineticsPredictorOutput.m tc0030, tc0069 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 diff --git a/test/unit_tests/kcatAssignmentTests.m b/test/unit_tests/kcatAssignmentTests.m new file mode 100644 index 00000000..3bdeb46b --- /dev/null +++ b/test/unit_tests/kcatAssignmentTests.m @@ -0,0 +1,151 @@ +%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'}) + [~, loaded] = loadKcatSource(base, fullfile(project, 'data', 'kcats', 'brenda.tsv')); + meta = loaded.header; + 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 8530c9dc..aa434321 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