From c1c7eee92514ea33976cccf2ca8a2b711377bc1a Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Fri, 23 Oct 2026 22:22:01 +0200 Subject: [PATCH 01/13] saveEcModel: drop the description assignment, which neither writer reads --- src/io/saveEcModel.m | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/io/saveEcModel.m b/src/io/saveEcModel.m index ae51dd54..9b1d8f65 100644 --- a/src/io/saveEcModel.m +++ b/src/io/saveEcModel.m @@ -59,8 +59,6 @@ function saveEcModel(ecModel,varargin) end filename = fullfile(params.path,'models',filename); -ecModel.description = ['Enzyme-constrained model of ' ecModel.id]; - switch filename(end-3:end) case {'.xml','sbml'} exportModel(ecModel, filename); From 682fe752654031f0d8a0008821392de355b69186 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Fri, 23 Oct 2026 22:27:46 +0200 Subject: [PATCH 02/13] ecFSEOF: select targets by the correlation with the enforced production, as RAVEN's FSEOF does A reaction is a target when its absolute flux correlates with the enforced production flux (absolute Pearson correlation of at least corrThreshold, default 0.9), instead of requiring a strictly monotonic flux and a slope above the 75th percentile. The targets are ranked by slope times correlation. --- src/flux_analysis/ecFSEOF.m | 107 ++++++++++++----------- test/unit_tests/geckoCoreFunctionTests.m | 37 ++++++++ 2 files changed, 95 insertions(+), 49 deletions(-) diff --git a/src/flux_analysis/ecFSEOF.m b/src/flux_analysis/ecFSEOF.m index 0ec57b93..d4b10c26 100644 --- a/src/flux_analysis/ecFSEOF.m +++ b/src/flux_analysis/ecFSEOF.m @@ -3,7 +3,9 @@ % % Runs Flux-Scanning with Enforced Objective Function (FSEOF) for a % specified production target, to identify reaction and gene targets for -% metabolic engineering. +% metabolic engineering. Reactions are classified as targets with the same +% rule as RAVEN's FSEOF: the correlation between the absolute flux of a +% reaction and the enforced production flux. % % Parameters % ---------- @@ -28,6 +30,10 @@ % modelAdapter, e.g. output/ecFSEOF_rxns.tsv). % modelAdapter : ModelAdapter % a loaded model adapter (default: the current default model adapter). +% corrThreshold : double +% minimum absolute Pearson correlation between the absolute flux of a +% reaction and the enforced production flux for the reaction to be a +% target (default 0.9). % % Returns % ------- @@ -40,25 +46,36 @@ % % - alpha : target production used for enforced objective limits (from minimum to maximum production). % - v_matrix : fluxes for each target reaction predicted and each alpha. -% - rxnTargets : a list with all reactions with fluxes that change consistently as target production increases. Contains ID, name, slope, gene rule, and equation. -% - transportTargets : a list with all transport reactions with fluxes that change consistently as target production increases. Contains ID, name, slope, gene rule, and equation. -% - geneTargets : a list with all selected targets that increase production. Contains gene, shortName, slope, action, and essentiality. +% - rxnTargets : a list with all reactions whose absolute flux is correlated with the enforced production, with an absolute correlation of at least corrThreshold. Contains ID, name, slope, gene rule, and equation. The slope is the absolute slope of the regression of the absolute flux on the enforced production flux. The list is ordered by the product of slope and correlation, highest first. +% - transportTargets : the same for transport reactions. +% - geneTargets : a list with all genes of the targets. Contains gene, shortName, slope, action, and essentiality. +% +% Only reactions with a gene rule are considered, and not those that are +% assigned to the standard enzyme. The enforced production reaction is +% included when it has a gene rule, as its flux follows the enforced +% production by definition. The reactions that draw enzyme usage +% (usage_prot_*) take part in the classification and are not reported. +% A reaction is over expressed (OE) when its flux increases with production, +% knocked out (KO) when it decreases and is at most 1e-8 at the last step, +% and knocked down (KD) when it decreases otherwise. % % Examples % -------- % % optional arguments may be given positionally or as name-value pairs: -% fseof = ecFSEOF(model,prodTargetRxn,csRxn,nSteps,outputFile,filePath,modelAdapter); +% fseof = ecFSEOF(model,prodTargetRxn,csRxn,nSteps,outputFile,filePath,modelAdapter,corrThreshold); % fseof = ecFSEOF(model,prodTargetRxn,csRxn,'nSteps',32); p = parseGECKOargs(varargin, { ... 'nSteps', 16; ... 'outputFile', false; ... 'filePath', []; ... - 'modelAdapter', []}); + 'modelAdapter', []; ... + 'corrThreshold',0.9}); nSteps = p.nSteps; outputFile = p.outputFile; filePath = p.filePath; -modelAdapter = p.modelAdapter; +modelAdapter = p.modelAdapter; +corrThreshold = p.corrThreshold; if isempty(modelAdapter) modelAdapter = ModelAdapterManager.getDefault(); @@ -133,66 +150,58 @@ v_matrix = v_matrix(withGR,:); rxnGeneM = model.rxnGeneMat(withGR,:); -% Filter out rxns that are always zero -zero_flux = ~all(abs(v_matrix(:,1:nSteps))==0,2); -target_rxns = target_rxns(zero_flux,:); -v_matrix = v_matrix(zero_flux,:); -rxnGeneM = rxnGeneM(zero_flux,:); - -% Identify those rxns that always increase or decrease, and calculate the -% slope as the difference in the flux when enforce objetive target -% production is set to 90%% of the maximum teorethical yield -% << v_matrix(i,nSteps-1) >> and the flux when the enforce objetive target -% production is set to the minimum << v_matrix(i,1) >> for the reaction i, -% divided by maxTarget-maxTarget/nSteps-1. +% Classify the reactions as in RAVEN's FSEOF. Slope and correlation are taken +% on the absolute flux against the enforced production flux, so a reversible +% reaction running backwards counts by its magnitude. A reaction is a target +% when some absolute flux exceeds 1e-10, the standard deviation of its +% absolute flux is above 0 and the absolute Pearson correlation reaches +% corrThreshold. It is suggested for over expression (OE) when the slope is +% positive, for knock out (KO) when the slope is negative and the absolute +% flux at the last step is at most 1e-8, and for knock down (KD) for any +% other negative slope. +absV = abs(v_matrix); slope_rxns = zeros(size(target_rxns)); -targets = logical(size(target_rxns)); +corr_rxns = zeros(size(target_rxns)); target_type = cell(size(target_rxns)); -for i = 1:length(target_rxns) - if issorted(abs(v_matrix(i,1:nSteps)),'strictascend') - % Those reactions that always increase while enforcing target - % production are suggested for Over Expression - targets(i) = true; - slope_rxns(i) = abs(v_matrix(i,nSteps-1)-v_matrix(i,1))/abs(maxTarget-maxTarget/nSteps-1); - target_type(i) = {'OE'}; - elseif issorted(abs(v_matrix(i,1:nSteps)),'strictdescend') - % Those reactions that always decrease while enforcing target - % production are suggested for KnockDown or KnockOut. KO are those - % reactions which have zero flux when enforcing target production - % to 90%% of the maximum theoretical yield. - targets(i) = true; - slope_rxns(i) = abs(v_matrix(i,nSteps-1)-v_matrix(i,1))/abs(maxTarget-maxTarget/nSteps-1); - if v_matrix(i,nSteps) == 0 - target_type(i) = {'KO'}; - else - target_type(i) = {'KD'}; +targets = false(size(target_rxns)); +if numel(alpha) > 1 && std(alpha) > 0 + for i = 1:length(target_rxns) + track = absV(i,:); + if any(track > 1e-10) && std(track) > 0 + fit = polyfit(alpha,track,1); + r = corrcoef(alpha,track); + if abs(r(1,2)) >= corrThreshold && fit(1) ~= 0 + targets(i) = true; + slope_rxns(i) = abs(fit(1)); + corr_rxns(i) = abs(r(1,2)); + if fit(1) > 0 + target_type(i) = {'OE'}; + elseif track(end) <= 1e-8 + target_type(i) = {'KO'}; + else + target_type(i) = {'KD'}; + end + end end end end -% Only keep those reaction that shows an increase or decrease pattern. +% Only keep those reactions that are classified as a target target_rxns = target_rxns(targets); v_matrix = v_matrix(targets,:); rxnGeneM = rxnGeneM(targets,:); slope_rxns = slope_rxns(targets); +corr_rxns = corr_rxns(targets); target_type = target_type(targets); -% Order from highest to lowest slope -[~,order] = sort(slope_rxns,'descend'); +% Order from highest to lowest score, the product of slope and correlation +[~,order] = sort(slope_rxns .* corr_rxns,'descend'); target_rxns = target_rxns(order); v_matrix = v_matrix(order,:); rxnGeneM = rxnGeneM(order,:); slope_rxns = slope_rxns(order); target_type = target_type(order); -% Filter out reactions with slope < quantile -non_zero_slope = slope_rxns > quantile(slope_rxns,0.75); -target_rxns = target_rxns(non_zero_slope); -v_matrix = v_matrix(non_zero_slope,:); -rxnGeneM = rxnGeneM(non_zero_slope,:); -slope_rxns = slope_rxns(non_zero_slope); -target_type = target_type(non_zero_slope); - % Create gene list of those connected to the remaining rxns genes = model.genes(sum(rxnGeneM,1) > 0); slope_genes = zeros(size(genes)); diff --git a/test/unit_tests/geckoCoreFunctionTests.m b/test/unit_tests/geckoCoreFunctionTests.m index d0a0e7fd..c1b1e0f1 100644 --- a/test/unit_tests/geckoCoreFunctionTests.m +++ b/test/unit_tests/geckoCoreFunctionTests.m @@ -2192,3 +2192,40 @@ function testGECKOInstallerIsOlderVersion_tc0059(testCase) verifyFalse(testCase, older('4.0.0', '4.0.0b1')) verifyFalse(testCase, older('4.0.0b1', '4.0.0b2')) end + +function testEcFSEOFSelectsTargetsByCorrelation_tc0099(testCase) + % Verifies ecFSEOF classifies a reaction as a target when the absolute + % flux correlates with the enforced production flux (default + % corrThreshold 0.9), as RAVEN's FSEOF does, and does not cut the list to + % the top quartile by slope. With R3 enforced, the flux of R5 falls as + % the shared enzyme pool goes to R3, and R3 itself follows the enforced + % flux; both are targets, the ones without a gene rule are not. + geckoPath = findGECKOroot; + adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM', 'TestGEMAdapter.m')); + model = getGeckoTestModel(); + model.geneShortNames = model.genes; + ecModel = makeEcModel(model, false, adapter); + ecModel = getECfromGEM(ecModel); + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R2_EXP_1')) = 50; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R2_EXP_2')) = 10; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R2_REV_EXP_1')) = 50; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R2_REV_EXP_2')) = 10; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R3')) = 20; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R5')) = 100; + ecModel.ec.source(:) = {'manual'}; + ecModel = applyKcatConstraints(ecModel); + ecModel = setProtPoolSize(ecModel, [], [], [], adapter); + adapter.params.bioRxn = 'R5'; + + fseof = ecFSEOF(ecModel, 'R3', 'S1', 'nSteps', 6, 'modelAdapter', adapter); + targets = sort([fseof.rxnTargets(:,1); fseof.transportTargets(:,1)]); + verifyEqual(testCase, targets, {'R3';'R5'}) + verifyEqual(testCase, numel(fseof.alpha), 6) + + % The slope column is the absolute regression slope, highest score first + verifyTrue(testCase, all(cell2mat(fseof.rxnTargets(:,3)) > 0)) + + % The default corrThreshold is 0.9 + sameAsDefault = ecFSEOF(ecModel, 'R3', 'S1', 'nSteps', 6, 'modelAdapter', adapter, 'corrThreshold', 0.9); + verifyEqual(testCase, sameAsDefault.rxnTargets, fseof.rxnTargets) +end From 587d43ba9571fd453bebce67305b74cc83d02a5c Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 24 Oct 2026 12:50:44 +0200 Subject: [PATCH 03/13] Proteomics filters, complex stoichiometry and kcat lists: loadProtData blanks the whole row above maxRSD and uses cutLowest, applyComplexData assigns the stoichiometry per protein ID, assignKcatValues drops the source of a zero kcat with its row, loadBRENDAdata finds lowercase EC codes --- src/change_model/applyComplexData.m | 8 +++- src/download_databases/loadBRENDAdata.m | 4 +- src/gather_kcats/assignKcatValues.m | 5 ++ src/io/loadProtData.m | 6 +-- test/unit_tests/kcatListTests.m | 57 ++++++++++++++++++++++ test/unit_tests/protDataTests.m | 63 +++++++++++++++++++++++++ 6 files changed, 137 insertions(+), 6 deletions(-) create mode 100644 test/unit_tests/kcatListTests.m create mode 100644 test/unit_tests/protDataTests.m diff --git a/src/change_model/applyComplexData.m b/src/change_model/applyComplexData.m index b5a60944..68bcc976 100644 --- a/src/change_model/applyComplexData.m +++ b/src/change_model/applyComplexData.m @@ -37,6 +37,9 @@ % % Notes % ----- +% The subunit stoichiometry is assigned per protein ID: the proteins of a +% complex need not be listed in the order of model.ec.enzymes. +% % foundComplex and proposedComplex share the following columns: rxn % (reaction identifier), complexID and name (Complex Portal identifiers), % genes (gene names of the complex subunits), protID_model (matched protein @@ -154,8 +157,11 @@ foundComplex{foundComplexCount,7} = complexData(complexMatch).stochiometry; % Apply in rxnEnzMat + % The stoichiometry is taken per protein ID, as the proteins of the + % complex need not be listed in the order of model.ec.enzymes [~, modelProtsIdx] = ismember(modelProts, model.ec.enzymes); - model.ec.rxnEnzMat(i,modelProtsIdx) = complexData(complexMatch).stochiometry; + [~, complexPos] = ismember(modelProts, complexData(complexMatch).protID); + model.ec.rxnEnzMat(i,modelProtsIdx) = complexData(complexMatch).stochiometry(complexPos); elseif modComplexUnits > 1 % Only suggest potential complex if grRule > 1 protein moreUnitsInData = find(percMatch == 1 & totalMatch > 1); diff --git a/src/download_databases/loadBRENDAdata.m b/src/download_databases/loadBRENDAdata.m index 6844ae25..ec69d64c 100644 --- a/src/download_databases/loadBRENDAdata.m +++ b/src/download_databases/loadBRENDAdata.m @@ -85,11 +85,11 @@ %build an index on MW{1} to speed things up a bit %first just extract the genus (i.e. the first part of the name) -MWECNum = upper(unique(MW{1})); +MWECNum = unique(upper(MW{1})); MWECNumIndices = cell(length(MWECNum),1); MWECNumHashMap = containers.Map(MWECNum,1:length(MWECNum)); for i = 1:length(MW{1}) - matchInd = cell2mat(values(MWECNumHashMap, MW{1}(i))); + matchInd = cell2mat(values(MWECNumHashMap, upper(MW{1}(i)))); MWECNumIndices{matchInd} = [MWECNumIndices{matchInd};i]; end diff --git a/src/gather_kcats/assignKcatValues.m b/src/gather_kcats/assignKcatValues.m index c7ec1c14..8188a3b0 100644 --- a/src/gather_kcats/assignKcatValues.m +++ b/src/gather_kcats/assignKcatValues.m @@ -21,6 +21,8 @@ % - genes : gene identifiers, matching model.ec.genes. % - substrates : substrate names, matching model.metNames. % - kcats : predicted kcat value in /sec. +% - kcatSource : optional, the source of each entry; entries with a zero +% kcat are discarded together with their source. % % Name-Value Arguments % -------------------- @@ -77,6 +79,9 @@ removeZero = kcatList.kcats == 0; kcatList.kcats(removeZero) = []; kcatList.rxns(removeZero) = []; +if isfield(kcatList,'kcatSource') + kcatList.kcatSource(removeZero) = []; +end % Map to model.ec.rxns [sanityCheck,idxInModel] = ismember(kcatList.rxns,model.ec.rxns); diff --git a/src/io/loadProtData.m b/src/io/loadProtData.m index bc3ea4b2..e4e3be27 100644 --- a/src/io/loadProtData.m +++ b/src/io/loadProtData.m @@ -143,16 +143,16 @@ end % Filter by RSD remData = (std(condAbund,0,2,'omitnan')./mean(condAbund,2,'omitnan'))>maxRSD; - condAbund(remData) = nan; + condAbund(remData,:) = nan; % Add stdevs condAbund = mean(condAbund,2,'omitnan')+(addStdevs*std(condAbund,0,2,'omitnan')); % Filter by minVal remData = mean(condAbund,2,'omitnan') Date: Sat, 24 Oct 2026 12:56:55 +0200 Subject: [PATCH 04/13] Gene rules, partial complexes, reaction id suffixes and isozymes without enzymes makeEcModel (light) and addNewRxnsToEC expand a gene rule to its isozymes in disjunctive normal form, so 'G1 and (G2 or G3)' is the isozymes {G1,G2} and {G1,G3}, and read gene ids as whole words, so an id such as orfX is not split at 'or'. makeEcModel couples a reaction to the subunits of its complex that the databases know, also when another subunit is missing from them. mapRxnsToConv and ecFVA strip the _REV and _EXP_n suffixes only at the end of a reaction id, so an id such as R_REVX is kept. applyKcatConstraints (light) ignores an isozyme row without enzymes when it selects the cheapest isozyme of a reaction. --- src/change_model/addNewRxnsToEC.m | 67 +++++++++++-- src/change_model/applyKcatConstraints.m | 13 +-- src/change_model/makeEcModel.m | 120 ++++++++++++------------ src/change_model/mapRxnsToConv.m | 9 +- src/flux_analysis/ecFVA.m | 2 +- test/unit_tests/ecModelBuildTests.m | 104 ++++++++++++++++++++ 6 files changed, 238 insertions(+), 77 deletions(-) create mode 100644 test/unit_tests/ecModelBuildTests.m diff --git a/src/change_model/addNewRxnsToEC.m b/src/change_model/addNewRxnsToEC.m index 7d7e99f5..7673c93b 100644 --- a/src/change_model/addNewRxnsToEC.m +++ b/src/change_model/addNewRxnsToEC.m @@ -35,7 +35,7 @@ % - rxns : cell array with unique strings that identifies each reaction. % - rxnNames : cell array with the names of each reaction. % - equations : cell array with equation strings. Decimal coefficients are expressed as "1.2". Reversibility is indicated by "<=>" or "=>". -% - grRules : cell array with the gene-reaction relationship for each reaction. For example "(A and B) or (C)" means that the reaction could be catalyzed by a complex between A & B or by C on its own. All the genes have to be present in model.genes. Add genes with addGenesRaven before calling this function if needed (opt, default ''). +% - grRules : cell array with the gene-reaction relationship for each reaction. For example "(A and B) or (C)" means that the reaction could be catalyzed by a complex between A & B or by C on its own. A rule with brackets is split into its isozymes in disjunctive normal form, so "A and (B or C)" is the two isozymes "A and B" and "A and C". All the genes have to be present in model.genes. Add genes with addGenesRaven before calling this function if needed (opt, default ''). % % The newEnzymes structure has the following fields: % @@ -139,8 +139,8 @@ % Before add the reactions validate that all the genes in grRules are % present in the model. genesInRules = strjoin(convertCharArray(newRxns.grRules)); - genesInRules = regexp(genesInRules,' |)|(|and|or','split'); % Remove all grRule punctuation - genesInRules = genesInRules(~cellfun(@isempty,genesInRules)); % Remove spaces and empty genes + genesInRules = regexp(genesInRules,'[^\s()]+','match'); % Gene ids and the words 'and' and 'or', without spaces and brackets + genesInRules = genesInRules(~ismember(genesInRules,{'and','or'})); % Only the gene ids genesInRules = setdiff(unique(genesInRules),model.genes); % Only keep new genes % If there is any gene which is not present, as to be pass as @@ -168,13 +168,14 @@ end % Get reactions to separate isozymes - toSplit = find(cellfun(@(x) contains(x,'or'),newRxns.grRules)); + isozymes = cellfun(@gprToDnf, newRxns.grRules, 'UniformOutput', false); + toSplit = find(cellfun(@numel,isozymes) > 1); if ~isempty(toSplit) splitRxns = {}; splitRxnNames = {}; splitEquations = {}; splitGrRules = {}; for i=1:numel(toSplit) idx = toSplit(i); - newRules = split(newRxns.grRules(idx),'or'); - newRules = strtrim(erase(newRules, {'(' ')'})); + % One reaction per AND-clause of the gene rule in disjunctive normal form + newRules = cellfun(@(x) strjoin(x,' and '), isozymes{idx}, 'UniformOutput', false); for j=1:numel(newRules) splitRxns{end+1} = [newRxns.rxns{idx} '_EXP_' num2str(j)]; @@ -260,3 +261,57 @@ rxnsAdded = newRxns.rxns(:); end end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%Function that expands a gene rule into disjunctive normal form: a cell array +%with one cell per isozyme, each holding the gene ids that are ANDed in it. +%Gene ids are whole words, so an id that contains 'or' or 'and' is kept +%intact. An empty rule gives an empty cell array. +function dnf = gprToDnf(rule) +tokens = regexp(rule,'\(|\)|[^\s()]+','match'); +if isempty(tokens) + dnf = {}; + return +end +[dnf,pos] = parseOr(tokens,1,rule); +if pos <= numel(tokens) + error('Cannot parse the gene rule "%s"',rule) +end +end + +function [dnf,pos] = parseOr(tokens,pos,rule) +[dnf,pos] = parseAnd(tokens,pos,rule); +while pos <= numel(tokens) && strcmp(tokens{pos},'or') + [rhs,pos] = parseAnd(tokens,pos+1,rule); + dnf = [dnf,rhs]; +end +end + +function [dnf,pos] = parseAnd(tokens,pos,rule) +[dnf,pos] = parseAtom(tokens,pos,rule); +while pos <= numel(tokens) && strcmp(tokens{pos},'and') + [rhs,pos] = parseAtom(tokens,pos+1,rule); + product = cell(1,0); + for a = 1:numel(dnf) + for b = 1:numel(rhs) + product{end+1} = [dnf{a},rhs{b}]; %#ok + end + end + dnf = product; +end +end + +function [dnf,pos] = parseAtom(tokens,pos,rule) +if pos > numel(tokens) || any(strcmp(tokens{pos},{')','and','or'})) + error('Cannot parse the gene rule "%s"',rule) +elseif strcmp(tokens{pos},'(') + [dnf,pos] = parseOr(tokens,pos+1,rule); + if pos > numel(tokens) || ~strcmp(tokens{pos},')') + error('Cannot parse the gene rule "%s"',rule) + end + pos = pos + 1; +else + dnf = {tokens(pos)}; + pos = pos + 1; +end +end diff --git a/src/change_model/applyKcatConstraints.m b/src/change_model/applyKcatConstraints.m index d8c25de2..d0536b7c 100644 --- a/src/change_model/applyKcatConstraints.m +++ b/src/change_model/applyKcatConstraints.m @@ -135,12 +135,13 @@ % summed per reaction), and divide by their kcat, to get a vector % of MW/kcat values. MWkcat = (model.ec.rxnEnzMat(ecIdx,:) * model.ec.mw) ./ model.ec.kcat(ecIdx); - % If kcat was zero, MWkcat is Inf. If no enzyme info was found, - % MWkcat is NaN. Drop those isozymes *before* selecting the - % cheapest one: correcting them to zero and including them in the - % min() below would let an isozyme with no kcat assigned always - % win over a real kcat on another isozyme of the same reaction. - validKcat = ~isinf(MWkcat) & ~isnan(MWkcat); + % If kcat was zero, MWkcat is Inf (or NaN without enzymes). An + % isozyme row without any enzyme has MWkcat zero. Drop those + % isozymes *before* selecting the cheapest one: including them in + % the min() below would let an isozyme with no kcat or no enzyme + % always win over a real kcat on another isozyme of the same reaction. + hasEnzyme = full(sum(model.ec.rxnEnzMat(ecIdx,:) ~= 0,2)) > 0; + validKcat = ~isinf(MWkcat) & ~isnan(MWkcat) & hasEnzyme; if any(validKcat) % Select the lowest MW/kcat (= most efficient), and convert to /hour model.S(prot_pool_idx, hasEc(i)) = -min(MWkcat(validKcat)/3600); diff --git a/src/change_model/makeEcModel.m b/src/change_model/makeEcModel.m index f8662b7e..b8737793 100644 --- a/src/change_model/makeEcModel.m +++ b/src/change_model/makeEcModel.m @@ -226,10 +226,11 @@ else %Different strategy for GECKO light: Each reaction can exist multiple times in %ec.rxns and similar fields - one time per isozyme. The number of copies is - %the number of ORs in the GPR + 1 - numOrs = count(model.grRules(rxnWithGene), ' or '); - cpys = numOrs + 1; - prevNumRxns = length(numOrs); + %the number of AND-clauses of the GPR in disjunctive normal form, so + %'G1 and (G2 or G3)' gives the isozymes {G1,G2} and {G1,G3}. + isozymes = cellfun(@gprToDnf, model.grRules(rxnWithGene), 'UniformOutput', false); + cpys = cellfun(@numel, isozymes); + prevNumRxns = numel(isozymes); cpyIndices = repelem(rxnWithGene, cpys).'; %.' only matter when number of rxns is 1 %loop through and add a prefix with an isozyme index to the rxns %we just give a fixed-length number as prefix, and assume that 999 is enough @@ -362,51 +363,18 @@ for r=1:numel(rxnWithGene) rxnGenes = model.genes(find(model.rxnGeneMat(rxnWithGene(r),:))); [~,locEnz] = ismember(rxnGenes,ec.genes); % Could also parse directly from rxnGeneMat, but some genes might be missing from Uniprot DB - if locEnz ~= 0 - ec.rxnEnzMat(r,locEnz) = 1; %Assume 1 copy per subunit or enzyme, can be modified later - end + locEnz = locEnz(locEnz ~= 0); % Subunits missing from the databases are left out; the others still couple the reaction + ec.rxnEnzMat(r,locEnz) = 1; %Assume 1 copy per subunit or enzyme, can be modified later end else %For light models, we need to split up all GPRs ec.rxnEnzMat = zeros(numel(ec.rxns),numel(ec.genes)); % Non-zeros will indicate the number of subunits nextIndex = 1; - %For full model generation, the GPRs are controlled in expandModel, but - %here we need to make an explicit format check - indexes2check = findPotentialErrors(model.grRules,model); - if ~isempty(indexes2check) - printOrange('Run standardizeGrRules(model) for a more detailed warning.\n'); - printOrange('For Human-GEM, these reactions can be corrected using simplifyGrRules.\n'); - end - for i=1:prevNumRxns - %ind is the index in the model, not to confuse with the index in the ec struct (i), - %which only contains reactions with GPRs. - ind = rxnWithGene(i); - %Get rid of all '(' and ')' since I'm not looking at complex stuff - %anyways - geneString=model.grRules{ind}; - geneString=strrep(geneString,'(',''); - geneString=strrep(geneString,')',''); - geneString=strrep(geneString,' or ',';'); - - if (numOrs(i) == 0) - geneNames = {geneString}; - else - %Split the string into gene names - geneNames=regexp(geneString,';','split'); - end - - %Now loop through the isozymes and set the rxnGeneMat - for j = 1:length(geneNames) - %Find the gene in the gene list If ' and ' relationship, first - %split the genes - fnd = strfind(geneNames{j},' and '); - if ~isempty(fnd) - andGenes=regexp(geneNames{j},' and ','split'); - ec.rxnEnzMat(nextIndex,ismember(ec.genes,andGenes)) = 1; %should be subunit stoichoimetry - else - ec.rxnEnzMat(nextIndex,ismember(ec.genes,geneNames(j)))=1;%should be subunit stoichoimetry - end + %Each isozyme is one AND-clause of the GPR and gets its own row of + %ec.rxnEnzMat, in the order of the copies of the reaction in ec.rxns. + for j = 1:numel(isozymes{i}) + ec.rxnEnzMat(nextIndex,ismember(ec.genes,isozymes{i}{j})) = 1;%should be subunit stoichoimetry nextIndex = nextIndex + 1; end end @@ -485,23 +453,55 @@ end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%Function that gets the model field grRules and returns the indexes of the -%rules in which the pattern ") and (" is present. -%Copied from standardizeGrRules -% TODO: Make this an accessible function in a separate file in RAVEN and remove this -%implementation. -function indexes2check = findPotentialErrors(grRules,model) -indxs_l = find(~cellfun(@isempty,strfind(grRules,') and ('))); -indxs_l_L = find(~cellfun(@isempty,strfind(grRules,') and'))); -indxs_l_R = find(~cellfun(@isempty,strfind(grRules,'and ('))); -indexes2check = vertcat(indxs_l,indxs_l_L,indxs_l_R); -indexes2check = unique(indexes2check); - -if ~isempty(indexes2check) - textToPrint = 'WARNING: Potentially problematic ") AND (" in the grRules for reaction(s):\n'; - for i=1:numel(indexes2check) - textToPrint=[textToPrint '\t' model.rxns{indexes2check(i)} '\n']; +%Function that expands a gene rule into disjunctive normal form: a cell array +%with one cell per isozyme, each holding the gene ids that are ANDed in it. +%Gene ids are whole words, so an id that contains 'or' or 'and' is kept +%intact. An empty rule gives an empty cell array. +function dnf = gprToDnf(rule) +tokens = regexp(rule,'\(|\)|[^\s()]+','match'); +if isempty(tokens) + dnf = {}; + return +end +[dnf,pos] = parseOr(tokens,1,rule); +if pos <= numel(tokens) + error('Cannot parse the gene rule "%s"',rule) +end +end + +function [dnf,pos] = parseOr(tokens,pos,rule) +[dnf,pos] = parseAnd(tokens,pos,rule); +while pos <= numel(tokens) && strcmp(tokens{pos},'or') + [rhs,pos] = parseAnd(tokens,pos+1,rule); + dnf = [dnf,rhs]; +end +end + +function [dnf,pos] = parseAnd(tokens,pos,rule) +[dnf,pos] = parseAtom(tokens,pos,rule); +while pos <= numel(tokens) && strcmp(tokens{pos},'and') + [rhs,pos] = parseAtom(tokens,pos+1,rule); + product = cell(1,0); + for a = 1:numel(dnf) + for b = 1:numel(rhs) + product{end+1} = [dnf{a},rhs{b}]; %#ok + end end - printOrange(textToPrint); + dnf = product; +end +end + +function [dnf,pos] = parseAtom(tokens,pos,rule) +if pos > numel(tokens) || any(strcmp(tokens{pos},{')','and','or'})) + error('Cannot parse the gene rule "%s"',rule) +elseif strcmp(tokens{pos},'(') + [dnf,pos] = parseOr(tokens,pos+1,rule); + if pos > numel(tokens) || ~strcmp(tokens{pos},')') + error('Cannot parse the gene rule "%s"',rule) + end + pos = pos + 1; +else + dnf = {tokens(pos)}; + pos = pos + 1; end end diff --git a/src/change_model/mapRxnsToConv.m b/src/change_model/mapRxnsToConv.m index e68f3dff..ee599626 100644 --- a/src/change_model/mapRxnsToConv.m +++ b/src/change_model/mapRxnsToConv.m @@ -43,11 +43,12 @@ rxnIDs = ecModel.rxns; % Invert flux of _REV reactions -revRxns = endsWith(rxnIDs,'_REV') | contains(rxnIDs,'_REV_EXP_'); +% The suffixes are added in the order _REV, _EXP_n, so strip _EXP_n first and +% then _REV, each only at the end of the id (an id like R_REVX is left intact). +rxnIDs = regexprep(rxnIDs,'_EXP_\d+$',''); +revRxns = endsWith(rxnIDs,'_REV'); fluxes(revRxns,:) = -fluxes(revRxns,:); -rxnIDs(revRxns) = replace(rxnIDs(revRxns),'_REV',''); -% Remove _EXP_. suffixes -rxnIDs = regexprep(rxnIDs,'_EXP_\d+',''); +rxnIDs(revRxns) = regexprep(rxnIDs(revRxns),'_REV$',''); % Map and sum fluxes to converted reaction IDs [rxnIDmap, convRxnID] = findgroups(rxnIDs); diff --git a/src/flux_analysis/ecFVA.m b/src/flux_analysis/ecFVA.m index dfee0ca4..e4aee30a 100644 --- a/src/flux_analysis/ecFVA.m +++ b/src/flux_analysis/ecFVA.m @@ -42,7 +42,7 @@ 'runParallel', true}); runParallel = p.runParallel; -rxnIDs = regexprep(ecModel.rxns,'(_REV)?(_EXP_\d+)?',''); +rxnIDs = regexprep(ecModel.rxns,'(_REV)?(_EXP_\d+)?$',''); [rxnIDmap, convRxnID] = findgroups(rxnIDs); N = numel(convRxnID); diff --git a/test/unit_tests/ecModelBuildTests.m b/test/unit_tests/ecModelBuildTests.m new file mode 100644 index 00000000..b1e39373 --- /dev/null +++ b/test/unit_tests/ecModelBuildTests.m @@ -0,0 +1,104 @@ +%run this test case with the command +%results = runtests('ecModelBuildTests.m') + +function tests = ecModelBuildTests + tests = functiontests(localfunctions); +end + +function testLightNestedGeneRuleExpandsToIsozymes_tc0210(testCase) + % 'G1 and (G2 or G3)' is the two isozymes {G1,G2} and {G1,G3}; the + % brackets must not be dropped, which would make the second one {G3}. + adapter = getAdapter(); + model = getGeckoTestModel(); + model = changeGrRules(model, {'R2'}, {'G1 and (G2 or G3)'}, true); + ecModel = makeEcModel(model, true, adapter); + verifyEqual(testCase, sort(coupledGenes(ecModel, '001_R2')), {'G1';'G2'}) + verifyEqual(testCase, sort(coupledGenes(ecModel, '002_R2')), {'G1';'G3'}) + verifyEqual(testCase, sum(strcmp(ecModel.ec.rxns, '003_R2')), 0) +end + +function testAddNewRxnsToECNestedRuleAndGeneIdWithOr_tc0211(testCase) + % The gene ids of the rule are whole words, so 'orfX' is a gene and not + % the operator 'or' followed by 'fX', and the nested rule is split into + % one reaction per isozyme: 'orfX and GNEW2' and 'orfX and GNEW3'. + adapter = getAdapter(); + model = getGeckoTestModel(); + ecModel = makeEcModel(model, false, adapter); + + newRxns.rxns = {'RNEWA'}; + newRxns.rxnNames = {'RNEWA'}; + newRxns.equations = {'m1[c] => e2[e]'}; + newRxns.grRules = {'orfX and (GNEW2 or GNEW3)'}; + newEnzymes.enzymes = {'ENEW1';'ENEW2';'ENEW3'}; + newEnzymes.genes = {'orfX';'GNEW2';'GNEW3'}; + newEnzymes.mw = [15000;16000;17000]; + + ecModel = addNewRxnsToEC(ecModel, newRxns, newEnzymes, adapter); + verifyEqual(testCase, ecModel.grRules(strcmp(ecModel.rxns, 'RNEWA_EXP_1')), {'orfX and GNEW2'}) + verifyEqual(testCase, ecModel.grRules(strcmp(ecModel.rxns, 'RNEWA_EXP_2')), {'orfX and GNEW3'}) +end + +function testMakeEcModelCouplesTheKnownSubunitsOfAComplex_tc0212(testCase) + % GX is in neither database: the reaction is still coupled to the enzyme + % of G4, as the databases know it. + adapter = getAdapter(); + model = getGeckoTestModel(); + model = changeGrRules(model, {'R3'}, {'G4 and GX'}, true); + ecModel = makeEcModel(model, false, adapter); + verifyEqual(testCase, coupledGenes(ecModel, 'R3'), {'G4'}) +end + +function testMapRxnsToConvStripsOnlyTheSuffixes_tc0213(testCase) + % R_REVX is a conventional reaction whose id contains _REV; its reverse + % variant R_REVX_REV carries flux that maps back to it with the opposite sign. + adapter = getAdapter(); + model = getGeckoTestModel(); + model.rxns{strcmp(model.rxns, 'R1')} = 'R_REVX'; + ecModel = makeEcModel(model, false, adapter); + fluxVect = zeros(numel(ecModel.rxns), 1); + fluxVect(strcmp(ecModel.rxns, 'R_REVX')) = 5; + fluxVect(strcmp(ecModel.rxns, 'R_REVX_REV')) = 2; + mappedFlux = mapRxnsToConv(ecModel, model, fluxVect); + verifyEqual(testCase, mappedFlux(strcmp(model.rxns, 'R_REVX')), 3) +end + +function testEcFVAStripsOnlyTheSuffixes_tc0214(testCase) + % The same reaction id in ecFVA, which groups the variants of each + % conventional reaction by id. + adapter = getAdapter(); + model = getGeckoTestModel(); + model.rxns{strcmp(model.rxns, 'R1')} = 'R_REVX'; + ecModel = makeEcModel(model, false, adapter); + [minFlux, maxFlux] = ecFVA(ecModel, model, 'runParallel', false); + verifySize(testCase, minFlux, [numel(model.rxns) 1]) + verifySize(testCase, maxFlux, [numel(model.rxns) 1]) +end + +function testApplyKcatConstraintsLightSkipsIsozymeWithoutEnzymes_tc0215(testCase) + % R2 has two isozymes; the second one has no enzyme but a kcat. It cannot + % be the cheapest isozyme: the cost of R2 is the one of the first. + adapter = getAdapter(); + model = getGeckoTestModel(); + lecModel = makeEcModel(model, true, adapter); + lecModel.ec.kcat(:) = 10; + first = strcmp(lecModel.ec.rxns, '001_R2'); + second = strcmp(lecModel.ec.rxns, '002_R2'); + lecModel.ec.rxnEnzMat(second,:) = 0; + lecModel = applyKcatConstraints(lecModel); + expected = -(lecModel.ec.rxnEnzMat(first,:) * lecModel.ec.mw) / 10 / 3600; + actual = full(lecModel.S(strcmp(lecModel.mets, 'prot_pool'), strcmp(lecModel.rxns, 'R2'))); + verifyGreaterThan(testCase, abs(expected), 0) + verifyEqual(testCase, actual, expected, 'AbsTol', 1e-10) +end + + +function adapter = getAdapter() +geckoPath = findGECKOroot; +adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); +end + +function genes = coupledGenes(ecModel, ecRxn) +% The genes of the enzymes that ec reaction ecRxn is coupled to. +row = strcmp(ecModel.ec.rxns, ecRxn); +genes = ecModel.ec.genes(logical(ecModel.ec.rxnEnzMat(row,:))); +end From fd2bf17e8528c5e134521056dcdc864f1364855f Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 24 Oct 2026 13:17:04 +0200 Subject: [PATCH 05/13] Verification matrix: list the tests of the changed functions; give the ecFSEOF test its own label --- test/VerificationMatrix.tsv | 15 ++++++++++----- test/unit_tests/geckoCoreFunctionTests.m | 2 +- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index d027bbbe..91c07dce 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -1,10 +1,10 @@ Function/File (in src) Verification Comments BrendaParser/* Not explicitly tested - The code has existed for a long time -change_model/applyComplexData.m tc0003, tc0004 -change_model/applyKcatConstraints.m tc0011 +change_model/applyComplexData.m tc0003, tc0004, tc0103 +change_model/applyKcatConstraints.m tc0011, tc0215 change_model/findMetSmiles.m tc0012 change_model/getComplexData.m tc0003, tc0004 - File download not covered by test cases but is tested manually and works. -change_model/makeEcModel.m tc0001, tc0002 +change_model/makeEcModel.m tc0001, tc0002, tc0210, tc0212 change_model/setProtPoolSize.m tc0005 change_model/getKcatAcrossIsoenzymes.m tc0011 change_model/applyCustomKcats.m @@ -13,7 +13,7 @@ gather_kcats/getStandardKcat.m gather_kcats/mergeDlkcatAndFuzzyKcats.m tc0011 gather_kcats/readDLKcatOutput.m tcm0001, tcm0002 - The test doesn't really check the function, only that some output is read 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/assignKcatValues.m tc0011 - we did not test all possible parameterizations here; also tc0104 gather_kcats/writeDLKcatInput.m partly by tc0011 gather_kcats/getEcRxnKeys.m tc0070, tc0076 gather_kcats/saveKcatSource.m tc0071, tc0073, tc0076, tc0078 @@ -32,6 +32,11 @@ gather_kcats/kcatsFromHomologs.m tc0095, tc0096 gather_kcats/extractEnzymeSubstratePairs.m tc0096, tc0015 gather_kcats/convertCustomKcats.m tc0091, tc0092, tc0094 change_model/upgradeEcModel.m tc0093, tc0094, tc0099 +io/loadProtData.m tc0101, tc0102 +change_model/addNewRxnsToEC.m tc0211 +change_model/mapRxnsToConv.m tc0213 +flux_analysis/ecFVA.m tc0214 +flux_analysis/ecFSEOF.m tc0100 utilities/upgradeGECKOproject.m tc0094, tc0099 gather_kcats/readOpenKineticsPredictorOutput.m tc0030, tc0069 gather_kcats/reviewKcatAssignment.m Has only been tested by running it manually, not tested by test cases @@ -41,7 +46,7 @@ get_enzyme_data/findECInDB.m Called from getECfromDatabase - tested when that fu get_enzyme_data/getECfromDatabase.m tc0007 - does not test download of the databases - this is tested in the manual workflows for Yeast-GEM and Human-GEM get_enzyme_data/getECfromGEM.m tc0006 get_enzyme_data/getECstring.m Called from findECInDB - tested when that function is tested. -download_databases/loadBRENDAdata.m Called from fuzzyKcatMatching - tested when that function is tested. +download_databases/loadBRENDAdata.m Called from fuzzyKcatMatching - tested when that function is tested; also tc0105 download_databases/fetchBrendaAsset.m tc0085 download_databases/loadDatabases.m Called from getECfromDatabase - tested when that function is tested. kcat_tuning/evotune/* Has only been tested by running it manually, not tested by test cases diff --git a/test/unit_tests/geckoCoreFunctionTests.m b/test/unit_tests/geckoCoreFunctionTests.m index c1b1e0f1..76b92d8f 100644 --- a/test/unit_tests/geckoCoreFunctionTests.m +++ b/test/unit_tests/geckoCoreFunctionTests.m @@ -2193,7 +2193,7 @@ function testGECKOInstallerIsOlderVersion_tc0059(testCase) verifyFalse(testCase, older('4.0.0b1', '4.0.0b2')) end -function testEcFSEOFSelectsTargetsByCorrelation_tc0099(testCase) +function testEcFSEOFSelectsTargetsByCorrelation_tc0100(testCase) % Verifies ecFSEOF classifies a reaction as a target when the absolute % flux correlates with the enforced production flux (default % corrThreshold 0.9), as RAVEN's FSEOF does, and does not cut the list to From 20da4e5680ad3ce3bbfb7df266cecee5bab4f56d Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 25 Oct 2026 07:46:26 +0100 Subject: [PATCH 06/13] Match geckopy on the BRENDA outlier filter, the fuzzy kcat matching and the custom kcats loadBRENDAdata drops a kcat more than four orders of magnitude below the median of its EC (filterOutlierKcats, default true). The reference of an EC is the median of the positive kcats of the max and median views of kcat.tsv and of the kcats derived from sa.tsv and mw.tsv, widened to the 3-, 2- and 1-level prefix while fewer than five values are available. On the BRENDA 2026.1 release it drops the same rows as geckopy's load_brenda_data. fuzzyKcatMatching escalates an EC number without a match to '-.-.-.-' (wildcard level 4), reports the wildcard level and origin of the selected EC when several EC numbers share the highest kcat, keeps the order and repeats of ecRxns given as indices, and returns no reactions for an empty selection other than []. kcatsFromBrenda refuses such a selection. applyCustomKcats skips an entry without a kcat, does not append an empty note, matches protein and reaction ids exactly and works on light ecModels. --- src/change_model/applyCustomKcats.m | 222 ++++++++++++---------- src/download_databases/loadBRENDAdata.m | 126 ++++++++++--- src/gather_kcats/fuzzyKcatMatching.m | 66 +++++-- src/gather_kcats/kcatsFromBrenda.m | 29 ++- test/VerificationMatrix.tsv | 8 +- test/unit_tests/geckoCoreFunctionTests.m | 14 +- test/unit_tests/kcatMatchTests.m | 224 +++++++++++++++++++++++ 7 files changed, 533 insertions(+), 156 deletions(-) create mode 100644 test/unit_tests/kcatMatchTests.m diff --git a/src/change_model/applyCustomKcats.m b/src/change_model/applyCustomKcats.m index 0cdd5e08..54f9ba67 100644 --- a/src/change_model/applyCustomKcats.m +++ b/src/change_model/applyCustomKcats.m @@ -43,9 +43,13 @@ % - gene_name : short gene name (optional, not used in matching). % - kcat : new kcat value (one per entry). % - rxns : reaction identifiers, multiple for the same kcat are separated by ',' (see further explanation below). -% - notes : will be appended to model.ec.notes (optional). +% - notes : will be appended to model.ec.notes, separated by ', ' from an existing note; an empty note changes nothing (optional). % - stoicho : complex stoichiometry, separated by ' + ' (examples: '1' or '3 + 1'), matching the order in proteins field. % +% An entry with an empty kcat is skipped. Protein and reaction identifiers +% are matched exactly, including their letter case. Light ecModels are +% supported: the rows of model.ec are matched by their ec reaction ids. +% % Matching order: % % 1. reactions are identified by .proteins and .rxns. @@ -121,125 +125,141 @@ enzInModel = cell(length(model.ec.rxns),1); ecRxnNoSuffix = regexprep(model.ec.rxns,'_EXP_\d+$',''); -% Implementation for full GECKO formulation -if ~model.ec.geckoLight - for i = 1:numel(customKcats.proteins) - if isempty(customKcats.proteins{i}) - %If only reaction ID(s) is/are specified (and no proteins), - %then apply the kcat to all isozymic reactions - rxns = strtrim(strsplit(customKcats.rxns{i}, ',')); - rxnIdxs = ismember(ecRxnNoSuffix,rxns); - rxnToUpdate(rxnIdxs) = 1; - model.ec.kcat(rxnIdxs) = customKcats.kcat(i); - - % Add note mentioning manual kcat change - model.ec.source(rxnIdxs) = {'custom'}; - if isfield(customKcats,'notes') - matchedIdxs = find(rxnIdxs); - for j = 1:numel(matchedIdxs) - if isempty(model.ec.notes{matchedIdxs(j), 1}) && ~isempty(customKcats.notes{i}) - model.ec.notes{matchedIdxs(j), 1} = customKcats.notes{i}; - else - model.ec.notes{matchedIdxs(j), 1} = [model.ec.notes{matchedIdxs(j), 1} ', ' customKcats.notes{i}]; - end - end +for i = 1:numel(customKcats.proteins) + if isnan(customKcats.kcat(i)) + printOrange(['WARNING: Entry ' num2str(i) ' has no kcat; it is skipped.']); + continue + end + if isempty(customKcats.proteins{i}) + %If only reaction ID(s) is/are specified (and no proteins), + %then apply the kcat to all isozymic reactions + rxns = strtrim(strsplit(customKcats.rxns{i}, ',')); + rxnIdxs = ismember(ecRxnNoSuffix,rxns); + rxnToUpdate(rxnIdxs) = 1; + model.ec.kcat(rxnIdxs) = customKcats.kcat(i); + + % Add note mentioning manual kcat change + model.ec.source(rxnIdxs) = {'custom'}; + if isfield(customKcats,'notes') + matchedIdxs = find(rxnIdxs); + for j = 1:numel(matchedIdxs) + model.ec.notes{matchedIdxs(j), 1} = appendNote(model.ec.notes{matchedIdxs(j), 1}, customKcats.notes{i}); end + end + else + prots = strtrim(strsplit(customKcats.proteins{i}, '+')); + + % Find the index for the enzymes which kcat will be changed. In + % case the protein is not in the model generate a warning. + try + enzIdx = cellfun(@(x) find(strcmp(model.ec.enzymes, x)), prots); + catch + enzIdx = []; + printOrange(['WARNING: Protein(s) ' customKcats.proteins{i} ' were not found in the model.']); + end + + % if not specific reactions are defined, find all the reaction + % index where the enzyme is used + if isempty(customKcats.rxns{i}) + temp_rxnIdxs = arrayfun(@(x) find(model.ec.rxnEnzMat(:, x)), enzIdx, 'UniformOutput', false); + % otherwhise, If a set of reactions if defined, only get the index for those + % but ignore any _EXP_ suffixes else - prots = strtrim(strsplit(customKcats.proteins{i}, '+')); - - % Find the index for the enzymes which kcat will be changed. In - % case the protein is not in the model generate a warning. - try - enzIdx = cellfun(@(x) find(strcmpi(model.ec.enzymes, x)), prots); - catch - enzIdx = []; - printOrange(['WARNING: Protein(s) ' customKcats.proteins{i} ' were not found in the model.']); - end + rxns = strtrim(strsplit(customKcats.rxns{i}, ',')); + temp_rxnIdxs = cellfun(@(x) find(strcmp(ecRxnNoSuffix, x)), rxns, 'UniformOutput', false); + end - % if not specific reactions are defined, find all the reaction - % index where the enzyme is used - if isempty(customKcats.rxns{i}) - temp_rxnIdxs = arrayfun(@(x) find(model.ec.rxnEnzMat(:, x)), enzIdx, 'UniformOutput', false); - % otherwhise, If a set of reactions if defined, only get the index for those - % but ignore any _EXP_ suffixes - else - rxns = strtrim(strsplit(customKcats.rxns{i}, ',')); - temp_rxnIdxs = arrayfun(@(x) find(strcmpi(ecRxnNoSuffix, x)), rxns, 'UniformOutput', false); + if ~isempty(temp_rxnIdxs) + rxnIdxs = []; + for j = 1:numel(temp_rxnIdxs) + rxnIdxs = [rxnIdxs; temp_rxnIdxs{j}]; end - if ~isempty(temp_rxnIdxs) - rxnIdxs = []; - for j = 1:numel(temp_rxnIdxs) - rxnIdxs = [rxnIdxs; temp_rxnIdxs{j}]; - end + % Check when multiple proteins are involved, since it can return same rxn n times + rxnIdxs = unique(rxnIdxs); %unique(rxnIdxs{1, :}); - % Check when multiple proteins are involved, since it can return same rxn n times - rxnIdxs = unique(rxnIdxs); %unique(rxnIdxs{1, :}); + for j = 1:numel(rxnIdxs) + % Get all the enzymes involved in the reaction + allEnzInRxn = find(model.ec.rxnEnzMat(rxnIdxs(j),:)); - for j = 1:numel(rxnIdxs) - % Get all the enzymes involved in the reaction - allEnzInRxn = find(model.ec.rxnEnzMat(rxnIdxs(j),:)); + C = intersect(enzIdx, allEnzInRxn); - C = intersect(enzIdx, allEnzInRxn); - - % Determine the match percentage bewteen the rules - if numel(C) == numel(enzIdx) && numel(C) == numel(allEnzInRxn) - match = 1; + % Determine the match percentage bewteen the rules + if numel(C) == numel(enzIdx) && numel(C) == numel(allEnzInRxn) + match = 1; + else + if numel(enzIdx) < numel(allEnzInRxn) + match = numel(C) / numel(allEnzInRxn); else - if numel(enzIdx) < numel(allEnzInRxn) - match = numel(C) / numel(allEnzInRxn); - else - match = numel(C) / numel(enzIdx); - end + match = numel(C) / numel(enzIdx); end + end - % Update the kcat value stored in the model, if match 100%, - % otherwhise if >= 0.5 inform to the user to be validated - if match == 1 - rxnToUpdate(rxnIdxs(j)) = 1; - model.ec.kcat(rxnIdxs(j)) = customKcats.kcat(i); - - % Add note mentioning manual kcat change - model.ec.source{rxnIdxs(j),1} = 'custom'; - if isfield(customKcats,'notes') - if isempty(model.ec.notes{rxnIdxs(j), 1}) && ~isempty(customKcats.notes{i}) - model.ec.notes{rxnIdxs(j), 1} = customKcats.notes{i}; - else - model.ec.notes{rxnIdxs(j), 1} = [model.ec.notes{rxnIdxs(j), 1} ', ' customKcats.notes{i}]; - end - end - elseif match >= 0.5 && match < 1 - rxnNotMatch(rxnIdxs(j)) = 1; - evaluatedRule{rxnIdxs(j), 1} = customKcats.proteins{i}; - enzInModel{rxnIdxs(j), 1} = strjoin(model.ec.enzymes(allEnzInRxn), ' + '); + % Update the kcat value stored in the model, if match 100%, + % otherwhise if >= 0.5 inform to the user to be validated + if match == 1 + rxnToUpdate(rxnIdxs(j)) = 1; + model.ec.kcat(rxnIdxs(j)) = customKcats.kcat(i); + + % Add note mentioning manual kcat change + model.ec.source{rxnIdxs(j),1} = 'custom'; + if isfield(customKcats,'notes') + model.ec.notes{rxnIdxs(j), 1} = appendNote(model.ec.notes{rxnIdxs(j), 1}, customKcats.notes{i}); end + elseif match >= 0.5 && match < 1 + rxnNotMatch(rxnIdxs(j)) = 1; + evaluatedRule{rxnIdxs(j), 1} = customKcats.proteins{i}; + enzInModel{rxnIdxs(j), 1} = strjoin(model.ec.enzymes(allEnzInRxn), ' + '); end end end end +end - % Apply the new kcat values to the model - if ~isempty(find(rxnToUpdate, 1)) - model = applyKcatConstraints(model, rxnToUpdate); +% Apply the new kcat values to the model +if ~isempty(find(rxnToUpdate, 1)) + if model.ec.geckoLight + %Light ecModels update by reaction of the model, not by ec row + updateRxns = ismember(model.rxns, extractAfter(model.ec.rxns(rxnToUpdate), 4)); else - printOrange('WARNING: No matches found. Consider checking the IDs or proteins in customKcats.'); + updateRxns = rxnToUpdate; end + model = applyKcatConstraints(model, updateRxns); +else + printOrange('WARNING: No matches found. Consider checking the IDs or proteins in customKcats.'); +end + +rxnUpdated = model.ec.rxns(find(rxnToUpdate)); - rxnUpdated = model.ec.rxns(find(rxnToUpdate)); - - % Remove from notMatch those reactions already updated - remove = and(rxnToUpdate, rxnNotMatch); - rxnNotMatch(remove) = 0; - evaluatedRule(remove) = ''; - enzInModel(remove) = ''; - - idRxns = model.ec.rxns(find(rxnNotMatch)); - fullIdx = cellfun(@(x) find(strcmpi(model.rxns, x)), idRxns); - rxnsNames = model.rxnNames(fullIdx); - evaluatedRule = evaluatedRule(~cellfun('isempty', evaluatedRule)); - enzInModel = enzInModel(~cellfun('isempty', enzInModel)); - rules = model.grRules(fullIdx); - notMatch = table(idRxns, rxnsNames, evaluatedRule, enzInModel, rules, ... - 'VariableNames',{'rxns', 'name', 'custom enzymes', 'enzymes in model', 'rules'}); +% Remove from notMatch those reactions already updated +remove = and(rxnToUpdate, rxnNotMatch); +rxnNotMatch(remove) = 0; +evaluatedRule(remove) = ''; +enzInModel(remove) = ''; + +idRxns = model.ec.rxns(find(rxnNotMatch)); +if model.ec.geckoLight + modelRxnIds = extractAfter(idRxns, 4); +else + modelRxnIds = idRxns; +end +fullIdx = cellfun(@(x) find(strcmpi(model.rxns, x), 1), modelRxnIds); +rxnsNames = model.rxnNames(fullIdx); +evaluatedRule = evaluatedRule(~cellfun('isempty', evaluatedRule)); +enzInModel = enzInModel(~cellfun('isempty', enzInModel)); +rules = model.grRules(fullIdx); +notMatch = table(idRxns, rxnsNames, evaluatedRule, enzInModel, rules, ... + 'VariableNames',{'rxns', 'name', 'custom enzymes', 'enzymes in model', 'rules'}); +end + +function notes = appendNote(notes, note) +% An empty note changes nothing; a note is separated by ', ' from an existing one. +note = strtrim(note); +if isempty(note) + return +elseif isempty(notes) + notes = note; +else + notes = [notes ', ' note]; end end diff --git a/src/download_databases/loadBRENDAdata.m b/src/download_databases/loadBRENDAdata.m index ec69d64c..07a96db4 100644 --- a/src/download_databases/loadBRENDAdata.m +++ b/src/download_databases/loadBRENDAdata.m @@ -21,6 +21,16 @@ % 'median' (default: params.kcatAggregation from the model adapter, % or 'max' if the adapter does not set it). A triple's max is the % highest turnover reported for it; its median is the middle one. +% filterOutlierKcats : logical +% whether to drop a row whose kcat is implausibly low for its EC +% (default: true). A row is dropped when its kcat is more than four +% orders of magnitude below the median kcat of its EC. The reference +% median of an EC is taken over the positive kcat values of the max and +% median views of kcat.tsv and of the kcat values derived from sa.tsv +% and mw.tsv, before any row is dropped, so it does not depend on +% kcatAggregation. Wildcard ECs (containing '-') do not contribute to a +% reference. The EC level widens from w.x.y.z to w.x.y, w.x and w while +% fewer than five values are available. Set to false to keep every row. % % Returns % ------- @@ -35,6 +45,7 @@ % [KCATcell, SAcell] = loadBRENDAdata(); % [KCATcell, SAcell] = loadBRENDAdata('modelAdapter', adapter); % [KCATcell, SAcell] = loadBRENDAdata('kcatAggregation', 'median'); +% [KCATcell, SAcell] = loadBRENDAdata('filterOutlierKcats', false); % % See also % -------- @@ -47,9 +58,14 @@ 'kcatAggregation', [], @(x) assert(isempty(x) || any(strcmpi(x, ... {'max','median'})), ... 'loadBRENDAdata:invalidKcatAggregation', ... - 'kcatAggregation must be ''max'' or ''median''.')}); -modelAdapter = p.modelAdapter; -kcatAggregation = p.kcatAggregation; + 'kcatAggregation must be ''max'' or ''median''.'); ... + 'filterOutlierKcats', true, @(x) assert((islogical(x) || isnumeric(x)) && ... + isscalar(x), ... + 'loadBRENDAdata:invalidFilterOutlierKcats', ... + 'filterOutlierKcats must be a logical scalar.')}); +modelAdapter = p.modelAdapter; +kcatAggregation = p.kcatAggregation; +filterOutlierKcats = logical(p.filterOutlierKcats); if isempty(modelAdapter) modelAdapter = ModelAdapterManager.getDefault(); @@ -59,7 +75,8 @@ end kcatAggregation = resolveKcatAggregation(kcatAggregation, modelAdapter); -%Column 4 holds the per-triple maximum and column 5 the median. +%Column 4 holds the per-triple maximum and column 5 the median; both are +%read, so that the outlier filter can use the max and the median views. valueCol = 4 + strcmpi(kcatAggregation,'median'); basePath = modelAdapter.getBrendaDBFolder(); @@ -69,20 +86,15 @@ %Extract BRENDA DATA from files information. kcat.tsv/sa.tsv are seven %columns (ec_code, substrate, organism, value_max, value_median, n, -%references); valueCol picks max or median. mw.tsv is six columns (no -%substrate-level aggregation choice): ec_code, substrate, organism, value, -%n, references. -KCATcell = openDataFile(KCAT_file,1,'%q %q %q %f %f %f %q',valueCol); +%references); the data cells hold ec_code, substrate, organism, the max and +%the median. mw.tsv is six columns (no substrate-level aggregation choice): +%ec_code, substrate, organism, value, n, references. +KCAT = openDataFile(KCAT_file,1,'%q %q %q %f %f %f %q',[4 5]); scalingFactor = 1/60; %[umol/min/mg] -> [mmol/s/g] -SA = openDataFile(SA_file,scalingFactor,'%q %q %q %f %f %f %q',valueCol); +SA = openDataFile(SA_file,scalingFactor,'%q %q %q %f %f %f %q',[4 5]); scalingFactor = 1/1000; %[g/mol] -> [g/mmol] MW = openDataFile(MW_file,scalingFactor,'%q %q %q %f %f %q',4); -for i=1:4 - SAcell{i} = []; -end -previousEC = []; EC_indexes = []; - %build an index on MW{1} to speed things up a bit %first just extract the genus (i.e. the first part of the name) MWECNum = unique(upper(MW{1})); @@ -93,7 +105,12 @@ MWECNumIndices{matchInd} = [MWECNumIndices{matchInd};i]; end - +%Match each specific activity to the molecular weight of its EC and +%organism. The first match is used. +previousEC = []; EC_indexes = []; +saRows = zeros(length(SA{1}),1); +mwVals = zeros(length(SA{1}),1); +nMatched = 0; for i=1:length(SA{1}) %Gets the indexes of the EC repetitions in the MW cell for every %new (different) EC @@ -111,18 +128,85 @@ % each EC# / Orgaism is reported on the file org_index = find(strcmpi(SA{3}(i),mwEC{1}),1); if ~isempty(org_index) - SAcell{1} = [SAcell{1};SA{1}(i)]; - SAcell{2} = [SAcell{2};SA{3}(i)]; - SAcell{3} = [SAcell{3}; SA{4}(i)*mwEC{2}(org_index)]; %[1/s] - SAcell{4} = [SAcell{4}; mwEC{2}(org_index)]; + nMatched = nMatched + 1; + saRows(nMatched) = i; + mwVals(nMatched) = mwEC{2}(org_index); end previousEC = SA{1}(i); end +saRows = saRows(1:nMatched); +mwVals = mwVals(1:nMatched); + +%Specific activities as kcat values [1/s], for the max and the median +saEC = SA{1}(saRows); +saMax = SA{4}(saRows).*mwVals; +saMedian = SA{5}(saRows).*mwVals; + +if valueCol == 4 + saSelected = saMax; +else + saSelected = saMedian; +end + +%Drop kcats that are implausibly low for their EC +keepKcat = true(size(KCAT{1})); +keepSA = true(size(saEC)); +if filterOutlierKcats + ecPool = [KCAT{1}; KCAT{1}; saEC; saEC]; + valPool = [KCAT{4}; KCAT{5}; saMax; saMedian]; + keepKcat = ~isOutlier(KCAT{1},KCAT{valueCol},ecPool,valPool); + keepSA = ~isOutlier(saEC,saSelected,ecPool,valPool); +end -function data_cell = openDataFile(fileName,scalingFactor,formatSpec,valueCol) +KCATcell = {KCAT{1}(keepKcat), KCAT{2}(keepKcat), KCAT{3}(keepKcat), ... + KCAT{valueCol}(keepKcat)}; +if nMatched > 0 + SAcell = {saEC(keepSA), SA{3}(saRows(keepSA)), saSelected(keepSA), ... + mwVals(keepSA)}; +else + SAcell = {[],[],[],[]}; +end +end + +function data_cell = openDataFile(fileName,scalingFactor,formatSpec,valueCols) fID = fopen(fileName); raw = textscan(fID,formatSpec,'delimiter','\t','HeaderLines',2); fclose(fID); -data_cell = {raw{1}, raw{2}, raw{3}, raw{valueCol}*scalingFactor}; +data_cell = [raw(1:3), cellfun(@(c) raw{c}*scalingFactor, num2cell(valueCols), ... + 'UniformOutput', false)]; end + +function outlier = isOutlier(ecs,vals,ecPool,valPool) +%A value is an outlier if it is more than 'decades' orders of magnitude below +%the median kcat of its own EC, taken over the pool of values. The EC widens +%from the full code to its 3-, 2- and 1-level prefix while fewer than minRef +%values are available. +decades = 4; +minRef = 5; +usePool = valPool > 0 & ~contains(ecPool,'-'); +ecPool = ecPool(usePool); +logPool = log10(valPool(usePool)); +refLog = NaN(size(ecs)); +if isempty(ecs) || isempty(logPool) + outlier = false(size(ecs)); + return +end +for level = 4:-1:1 + poolKey = ecPrefix(ecPool,level); + [grp,keys] = findgroups(poolKey); + nPerKey = accumarray(grp,1); + medPerKey = splitapply(@median,logPool,grp); + [found,loc] = ismember(ecPrefix(ecs,level),keys); + use = isnan(refLog) & found; + use(use) = nPerKey(loc(use)) >= minRef; + refLog(use) = medPerKey(loc(use)); +end +outlier = vals > 0 & ~isnan(refLog) & log10(vals) < refLog - decades; +end + +function prefix = ecPrefix(ecs,level) +%The first 'level' dot-separated parts of each EC code +prefix = regexp(ecs,['^[^.]*(?:\.[^.]*){0,' num2str(level-1) '}'],'match','once'); +prefix(cellfun(@isempty,prefix)) = {''}; +prefix = reshape(prefix,size(ecs)); end diff --git a/src/gather_kcats/fuzzyKcatMatching.m b/src/gather_kcats/fuzzyKcatMatching.m index 9956ee8e..ed3e1a05 100644 --- a/src/gather_kcats/fuzzyKcatMatching.m +++ b/src/gather_kcats/fuzzyKcatMatching.m @@ -18,10 +18,15 @@ % % Name-Value Arguments % -------------------- -% ecRxns : logical +% ecRxns : logical or double % for which reactions (from model.ec.rxns) kcat values should be found, -% provided as logical vector with same length as model.ec.rxns -% (default: all reactions). +% provided as a logical vector with the same length as model.ec.rxns +% (the reactions are returned in model order), or as indices into +% model.ec.rxns (the reactions are returned in the order given, and an +% index given twice is returned twice). The literal [] selects all +% reactions, which is the default; any other empty selection, such as +% find(mask) without true entries, selects none and returns an empty +% kcatList. % modelAdapter : ModelAdapter % a loaded model adapter (default: the current default model adapter). % forceWClvl : double @@ -43,7 +48,8 @@ % - kcats : proposed kcat value in /sec. % - eccodes : as used to query BRENDA. % - wildcardLvl : which level of EC wild-card was necessary to find a match -% (0: w.x.y.z; 1: w.x.y.-; 2: w.x.-.-; 3: w.-.-.-). +% (0: w.x.y.z; 1: w.x.y.-; 2: w.x.-.-; 3: w.-.-.-; 4: -.-.-.-, which +% matches every EC number in BRENDA). % - origin : which level of specificity was necessary to find a match % (1: correct organism, correct substrate, kcat; 2: any organism, correct % substrate, kcat; 3: correct organism, any substrate, kcat; 4: any @@ -52,7 +58,9 @@ % 1, 2, 3, 5, 4, 6: a specific activity from the correct organism is % preferred over a kcat from any organism with any substrate. When a % reaction has multiple EC numbers, the EC with the fewest wildcards and -% the earliest-tried level is selected, and among those the highest kcat. +% the earliest-tried level is selected, and among those the highest kcat +% (the first one when several EC numbers share it); the wildcardLvl and +% origin are those of the selected EC. % % If a wildcard is used, origin levels 1 and 2 are ignored. The last digits % in the E.C. number indicate the substrate specificity, so if this should @@ -93,14 +101,19 @@ forceWClvl = p.forceWClvl; kcatAggregation = p.kcatAggregation; -if isempty(ecRxns) - ecRxns = true(numel(model.ec.rxns),1); -elseif isnumeric(ecRxns) - ecRxnsVec = false(numel(model.ec.rxns),1); - ecRxnsVec(ecRxns) = true; - ecRxns = ecRxnsVec; +% Get indices: [] selects all reactions, a logical vector the true entries (in +% model order) and indices are used as given (order and repeats kept). +if isnumeric(ecRxns) && isequal(size(ecRxns),[0 0]) + ecRxns = (1:numel(model.ec.rxns)).'; +elseif islogical(ecRxns) + ecRxns = find(ecRxns(:)); +else + ecRxns = ecRxns(:); + if any(ecRxns < 1 | ecRxns > numel(model.ec.rxns) | ecRxns ~= round(ecRxns)) + error('fuzzyKcatMatching:invalidEcRxns', ... + 'ecRxns must hold indices of model.ec.rxns.'); + end end -ecRxns=find(ecRxns); % Get indices instead of logical if isempty(modelAdapter) modelAdapter = ModelAdapterManager.getDefault(); @@ -118,6 +131,17 @@ error('No EC codes defined in model.ec.eccodes. First run getECfromGEM() and/or getECfromDatabase().') end eccodes = model.ec.eccodes(ecRxns); +if isempty(ecRxns) + %No reactions selected: nothing to match, and BRENDA is not read + kcatList.source = 'brenda'; + kcatList.rxns = cell(0,1); + kcatList.substrates = cell(0,1); + kcatList.kcats = zeros(0,1); + kcatList.eccodes = cell(0,1); + kcatList.wildcardLvl = NaN(0,1); + kcatList.origin = NaN(0,1); + return +end substrates = cell(numel(ecRxns),1); substrCoeffs = cell(numel(ecRxns),1); @@ -264,6 +288,9 @@ %wild_num==4, an out-of-range index. kcat(k)/origin(k) keep %their zero-initialised "no match" values. success = true; + elseif wild_num == 3 + %Last escalation: '-.-.-.-' matches every EC number in BRENDA + EC{k} = '-.-.-.-'; else dot_pos = [2 strfind(EC{k},'.')]; wc_text = '-.-.-.-'; @@ -284,11 +311,14 @@ %(higher) wildcard count back in whenever it happens to share best_origin, %defeating the wildcard-count filter entirely. best_pos = best_pos & (origin == best_origin); - max_pos = find(kcat == max(kcat(best_pos))); - wc_num = wc_num(max_pos(1)); - origin = origin(max_pos(1)); - matches = matches(max_pos(1)); - kcat = kcat(max_pos(1)); + %Highest kcat among those, the first one if several share it + best_idx = find(best_pos); + [~,j] = max(kcat(best_idx)); + max_pos = best_idx(j); + wc_num = wc_num(max_pos); + origin = origin(max_pos); + matches = matches(max_pos); + kcat = kcat(max_pos); %Update dir and tot: dir.org_s(i) = matches*(origin == 1); @@ -465,6 +495,8 @@ for j = 1:length(X) EC_indexes = [EC_indexes,EcIndexIndices{X(j)}]; end + elseif isempty(EC) %Not optimized, fully wildcarded: every entry matches + EC_indexes = 1:length(EC_cell); else %Not optimized for j=1:length(EC_cell) if strfind(EC_cell{j},EC)==1 diff --git a/src/gather_kcats/kcatsFromBrenda.m b/src/gather_kcats/kcatsFromBrenda.m index 64d29428..9ced7d30 100644 --- a/src/gather_kcats/kcatsFromBrenda.m +++ b/src/gather_kcats/kcatsFromBrenda.m @@ -20,8 +20,10 @@ % 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). +% ecRxns : logical or double +% for which model.ec.rxns kcats are gathered (default: all). A logical +% vector or indices of model.ec.rxns, see fuzzyKcatMatching; [] selects +% all reactions and any other empty selection is an error. % save : logical % write the source file (default true). % saveAs : char @@ -60,6 +62,10 @@ end end aggregation = resolveKcatAggregation(p.aggregation, modelAdapter); +if isSubset(p.ecRxns) && countSelected(p.ecRxns) == 0 + error('kcatsFromBrenda:emptySubset', ... + 'ecRxns selects no reactions; leave it out to gather for all reactions.'); +end if p.save file = fullfile(modelAdapter.params.path, 'data', 'kcats', [char(p.saveAs) '.tsv']); warnSubset(p.ecRxns, numel(model.ec.rxns), file); @@ -70,8 +76,8 @@ if p.save parameters = struct('aggregation', aggregation, 'force_wildcard_level', p.forceWClvl); - if ~isempty(p.ecRxns) - parameters.ec_rxns = countRxns(p.ecRxns); + if isSubset(p.ecRxns) + parameters.ec_rxns = numel(kcatList.rxns); end saveKcatSource(model, kcatList, p.saveAs, 'modelAdapter', modelAdapter, ... 'generator', 'kcatsFromBrenda', 'parameters', parameters, ... @@ -96,6 +102,19 @@ end end +function tf = isSubset(ecRxns) +% [] selects all reactions; any other value selects the reactions it names. +tf = ~(isnumeric(ecRxns) && isequal(size(ecRxns),[0 0])); +end + +function n = countSelected(ecRxns) +if islogical(ecRxns) + n = sum(ecRxns); +else + n = numel(ecRxns); +end +end + function n = countRxns(ecRxns) if islogical(ecRxns) n = sum(ecRxns); @@ -106,7 +125,7 @@ function warnSubset(ecRxns, nRxns, file) % A subset written under the default name replaces the full file. -if ~isempty(ecRxns) && countRxns(ecRxns) < nRxns && isfile(file) +if isSubset(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 diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 91c07dce..6b036e07 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -7,8 +7,8 @@ change_model/getComplexData.m tc0003, tc0004 - File download not covered by test change_model/makeEcModel.m tc0001, tc0002, tc0210, tc0212 change_model/setProtPoolSize.m tc0005 change_model/getKcatAcrossIsoenzymes.m tc0011 -change_model/applyCustomKcats.m -gather_kcats/fuzzyKcatMatching.m tc0010 +change_model/applyCustomKcats.m tc0227, tc0228, tc0229, tc0230 +gather_kcats/fuzzyKcatMatching.m tc0010, tc0223, tc0224, tc0225, tc0226 gather_kcats/getStandardKcat.m gather_kcats/mergeDlkcatAndFuzzyKcats.m tc0011 gather_kcats/readDLKcatOutput.m tcm0001, tcm0002 - The test doesn't really check the function, only that some output is read @@ -18,7 +18,7 @@ gather_kcats/writeDLKcatInput.m partly by tc0011 gather_kcats/getEcRxnKeys.m tc0070, tc0076 gather_kcats/saveKcatSource.m tc0071, tc0073, tc0076, tc0078 gather_kcats/loadKcatSource.m tc0071, tc0072, tc0073, tc0076, tc0077, tc0078 -gather_kcats/kcatsFromBrenda.m tc0074, tc0078 +gather_kcats/kcatsFromBrenda.m tc0074, tc0078, tc0226 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, tc0087, tc0088, tc0089, tc0090, tc0098 @@ -46,7 +46,7 @@ get_enzyme_data/findECInDB.m Called from getECfromDatabase - tested when that fu get_enzyme_data/getECfromDatabase.m tc0007 - does not test download of the databases - this is tested in the manual workflows for Yeast-GEM and Human-GEM get_enzyme_data/getECfromGEM.m tc0006 get_enzyme_data/getECstring.m Called from findECInDB - tested when that function is tested. -download_databases/loadBRENDAdata.m Called from fuzzyKcatMatching - tested when that function is tested; also tc0105 +download_databases/loadBRENDAdata.m Called from fuzzyKcatMatching - tested when that function is tested; also tc0105, tc0220, tc0221, tc0222 download_databases/fetchBrendaAsset.m tc0085 download_databases/loadDatabases.m Called from getECfromDatabase - tested when that function is tested. kcat_tuning/evotune/* Has only been tested by running it manually, not tested by test cases diff --git a/test/unit_tests/geckoCoreFunctionTests.m b/test/unit_tests/geckoCoreFunctionTests.m index 76b92d8f..1ae0c61e 100644 --- a/test/unit_tests/geckoCoreFunctionTests.m +++ b/test/unit_tests/geckoCoreFunctionTests.m @@ -1270,28 +1270,26 @@ function testFindECInDBIntersectionDedupesWildcardPair_tc0032(testCase) function testFuzzyKcatMatchingWildcardExhaustionDoesNotCrash_tc0033(testCase) % Verifies that when an EC code has no BRENDA match at any wildcard - % level, including the fully wildcarded '-.-.-.-', fuzzyKcatMatching's - % escalation loop terminates cleanly and reports "no match" (kcat 0, - % origin NaN) instead of erroring. + % level, including the fully wildcarded '-.-.-.-' (which matches every + % EC number in BRENDA), fuzzyKcatMatching's escalation loop terminates + % cleanly and reports "no match" (kcat 0, origin NaN) instead of + % erroring. The fixture BRENDA kcat table has no rows, so nothing can match. geckoPath = findGECKOroot; adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM', 'TestGEMAdapter.m')); model = getGeckoTestModel(); ecModel = makeEcModel(model, false, adapter); ecModel = getECfromGEM(ecModel); - % An EC class unrelated, at every wildcard level, to anything in the - % fixture BRENDA file below. ecModel.ec.eccodes{strcmp(ecModel.ec.rxns, 'R2_EXP_1')} = '9.9.9.9'; scratchDir = fullfile(tempname); brendaDir = fullfile(scratchDir, 'data'); mkdir(brendaDir); cleanupDir = onCleanup(@() rmdir(scratchDir, 's')); - % Deliberately unrelated to '9.9.9.9' at every wildcard level, so the - % escalation loop runs out of levels without ever matching. + % No rows at all, so the escalation loop runs out of levels without + % ever matching. fid = fopen(fullfile(brendaDir, 'kcat.tsv'), 'w'); fprintf(fid, '# BRENDA release test-fixture - kcat in 1/s\n'); fprintf(fid, 'ec_code\tsubstrate\torganism\tkcat_max\tkcat_median\tn\treferences\n'); - fprintf(fid, '5.5.5.5\tm1\ttestus testus\t42\t42\t1\t*\n'); fclose(fid); fid = fopen(fullfile(brendaDir, 'mw.tsv'), 'w'); fprintf(fid, '# BRENDA release test-fixture - molecular weight in g/mol\n'); diff --git a/test/unit_tests/kcatMatchTests.m b/test/unit_tests/kcatMatchTests.m new file mode 100644 index 00000000..abc4471a --- /dev/null +++ b/test/unit_tests/kcatMatchTests.m @@ -0,0 +1,224 @@ +%run this test case with the command +%results = runtests('kcatMatchTests.m') + +function tests = kcatMatchTests + tests = functiontests(localfunctions); +end + +function testLoadBRENDAdataDropsKcatFarBelowItsEc_tc0220(testCase) + % A kcat more than four orders of magnitude below the median of its EC + % is dropped by default and kept when the filter is switched off. + adapter = brendaAdapter(testCase, ... + [kcatRows('1.1.1.1', {'o1','o2','o3','o4','o5','o6'}, 10, 10), ... + kcatRows('1.1.1.1', {'o7'}, 1e-4, 1e-4)]); + KCATcell = loadBRENDAdata('modelAdapter', adapter); + verifyEqual(testCase, KCATcell{4}, repmat(10,6,1)) + KCATcell = loadBRENDAdata('modelAdapter', adapter, 'filterOutlierKcats', false); + verifyEqual(testCase, sort(KCATcell{4}), sort([repmat(10,6,1); 1e-4])) +end + +function testLoadBRENDAdataReferenceUsesMaxAndMedian_tc0221(testCase) + % The reference of an EC is built from the max and the median values, + % whichever aggregation is read: a median that is far below its EC is an + % outlier even when the max of the same row is not. + adapter = brendaAdapter(testCase, ... + [kcatRows('2.1.1.1', {'o1','o2','o3','o4'}, 100, 100), ... + kcatRows('2.1.1.1', {'o5'}, 100, 0.001)]); + KCATcell = loadBRENDAdata('modelAdapter', adapter, 'kcatAggregation', 'median'); + verifyEqual(testCase, KCATcell{4}, repmat(100,4,1)) + KCATcell = loadBRENDAdata('modelAdapter', adapter, 'kcatAggregation', 'max'); + verifyEqual(testCase, KCATcell{4}, repmat(100,5,1)) +end + +function testLoadBRENDAdataReferenceNeedsFiveValues_tc0222(testCase) + % The values of the max and the median of a row both count towards the + % reference of its EC. With fewer than five of them at every EC level + % there is no reference and no row is dropped; with five or more at a + % shorter EC prefix, that prefix is the reference. + adapter = brendaAdapter(testCase, ... + [kcatRows('3.1.1.1', {'o1'}, 10, 10), kcatRows('3.1.1.1', {'o2'}, 1e-9, 1e-9)]); + KCATcell = loadBRENDAdata('modelAdapter', adapter); + verifyEqual(testCase, numel(KCATcell{4}), 2) + adapter = brendaAdapter(testCase, ... + [kcatRows('4.1.1.2', {'o1','o2','o3'}, 10, 10), kcatRows('4.1.1.1', {'o4'}, 1e-9, 1e-9)]); + KCATcell = loadBRENDAdata('modelAdapter', adapter); + verifyEqual(testCase, KCATcell{4}, repmat(10,3,1)) +end + +function testFuzzyKcatMatchingFullWildcard_tc0223(testCase) + % An EC number without any match at the wildcard levels 1 to 3 is matched + % against every EC number in BRENDA, at wildcard level 4. + [ecModel, fuzzyAdapter] = fuzzyModel(testCase, {'5.5.5.5\tm1\ttestus testus\t42\t42\t1\t*\n'}); + ecModel.ec.eccodes{strcmp(ecModel.ec.rxns, 'R2_EXP_1')} = '9.9.9.9'; + kcatList = fuzzyKcatMatching(ecModel, ismember(ecModel.ec.rxns,{'R2_EXP_1'}), fuzzyAdapter); + verifyEqual(testCase, kcatList.kcats, 42) + verifyEqual(testCase, kcatList.wildcardLvl, 4) +end + +function testFuzzyKcatMatchingTieReportsTheSelectedEc_tc0224(testCase) + % Two EC numbers of a reaction give the same kcat, one through a wildcard + % and one exactly: the wildcard level and origin are those of the exact one. + [ecModel, fuzzyAdapter] = fuzzyModel(testCase, ... + {'1.1.1.1\tm1\ttestus testus\t7\t7\t1\t*\n', ... + '1.1.1.9\tm1\ttestus testus\t7\t7\t1\t*\n'}); + ecModel.ec.eccodes{strcmp(ecModel.ec.rxns, 'R2_EXP_1')} = '1.1.1.5;1.1.1.1'; + kcatList = fuzzyKcatMatching(ecModel, ismember(ecModel.ec.rxns,{'R2_EXP_1'}), fuzzyAdapter); + verifyEqual(testCase, kcatList.kcats, 7) + verifyEqual(testCase, kcatList.wildcardLvl, 0) + verifyEqual(testCase, kcatList.origin, 1) +end + +function testFuzzyKcatMatchingSubsetKeepsOrderAndRepeats_tc0225(testCase) + % Indices select the reactions in the order given, a repeated index + % gives a repeated reaction, and a logical vector gives model order. + [ecModel, fuzzyAdapter] = fuzzyModel(testCase, {'1.1.1.1\tm1\ttestus testus\t7\t7\t1\t*\n'}); + ids = {'R3','R2_EXP_1','R2_EXP_1'}; + idx = cellfun(@(r) find(strcmp(ecModel.ec.rxns, r)), ids); + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', idx, 'modelAdapter', fuzzyAdapter); + verifyEqual(testCase, kcatList.rxns(:)', ids) + verifyEqual(testCase, numel(kcatList.kcats), 3) + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', ismember(ecModel.ec.rxns, ids), ... + 'modelAdapter', fuzzyAdapter); + verifyEqual(testCase, kcatList.rxns, ecModel.ec.rxns(ismember(ecModel.ec.rxns, ids))) +end + +function testFuzzyKcatMatchingEmptySubsetSelectsNone_tc0226(testCase) + % [] selects all reactions; any other empty selection selects none. + [ecModel, fuzzyAdapter] = fuzzyModel(testCase, {'1.1.1.1\tm1\ttestus testus\t7\t7\t1\t*\n'}); + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', [], 'modelAdapter', fuzzyAdapter); + verifyEqual(testCase, numel(kcatList.rxns), numel(ecModel.ec.rxns)) + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', zeros(0,1), 'modelAdapter', fuzzyAdapter); + verifyEmpty(testCase, kcatList.rxns) + verifyEmpty(testCase, kcatList.kcats) + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', false(numel(ecModel.ec.rxns),1), ... + 'modelAdapter', fuzzyAdapter); + verifyEmpty(testCase, kcatList.rxns) + adapter = getTestAdapter(); + verifyError(testCase, @() kcatsFromBrenda(ecModel, 'ecRxns', zeros(0,1), 'save', false, ... + 'modelAdapter', adapter), 'kcatsFromBrenda:emptySubset') +end + +function testApplyCustomKcatsSkipsABlankKcat_tc0227(testCase) + % An entry without a kcat leaves the reaction as it is. + [ecModel, adapter] = customModel(); + isR3 = strcmp(ecModel.ec.rxns,'R3'); + isR2 = strcmp(ecModel.ec.rxns,'R2_EXP_1'); + custom = customStruct({''; ''}, [NaN; 50], {'R3'; 'R2'}, {'';''}); + ecModel = applyCustomKcats(ecModel, custom, adapter); + verifyEqual(testCase, ecModel.ec.kcat(isR3), 10) + verifyFalse(testCase, strcmp(ecModel.ec.source{isR3}, 'custom')) + verifyEqual(testCase, ecModel.ec.kcat(isR2), 50) + verifyEqual(testCase, ecModel.ec.source{isR2}, 'custom') +end + +function testApplyCustomKcatsEmptyNoteIsNotAppended_tc0228(testCase) + % An empty note does not add a separator to the note of the reaction. + [ecModel, adapter] = customModel(); + custom = customStruct({''; ''}, [50; 60], {'R3'; 'R3'}, {'first'; ''}); + ecModel = applyCustomKcats(ecModel, custom, adapter); + verifyEqual(testCase, ecModel.ec.notes{strcmp(ecModel.ec.rxns,'R3')}, 'first') + custom = customStruct({''; ''}, [50; 60], {'R3'; 'R3'}, {''; 'second'}); + ecModel = applyCustomKcats(customModelAgain(), custom, adapter); + verifyEqual(testCase, ecModel.ec.notes{strcmp(ecModel.ec.rxns,'R3')}, 'second') +end + +function testApplyCustomKcatsMatchesIdsExactly_tc0229(testCase) + % Protein and reaction ids are case sensitive. + [ecModel, adapter] = customModel(); + before = ecModel.ec.kcat; + custom = customStruct({'p3'}, 123, {''}, {''}); + updated = applyCustomKcats(ecModel, custom, adapter); + verifyEqual(testCase, updated.ec.kcat, before) + custom = customStruct({'P3'}, 123, {'r3'}, {''}); + updated = applyCustomKcats(ecModel, custom, adapter); + verifyEqual(testCase, updated.ec.kcat, before) + custom = customStruct({'P3'}, 123, {''}, {''}); + updated = applyCustomKcats(ecModel, custom, adapter); + verifyTrue(testCase, any(updated.ec.kcat == 123)) +end + +function testApplyCustomKcatsLightModel_tc0230(testCase) + % A light ecModel is updated by the proteins of its rows, and all three + % outputs are returned. + adapter = getTestAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), true, adapter); + ecModel = getECfromGEM(ecModel); + custom = customStruct({'P3'}, 123, {''}, {'light'}); + [updated, rxnUpdated, notMatch] = applyCustomKcats(ecModel, custom, adapter); + verifyNotEmpty(testCase, rxnUpdated) + verifyTrue(testCase, all(updated.ec.kcat(ismember(updated.ec.rxns, rxnUpdated)) == 123)) + verifyTrue(testCase, istable(notMatch)) +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function adapter = getTestAdapter() + geckoPath = findGECKOroot; + adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); +end + +function rows = kcatRows(ec, organisms, kcatMax, kcatMedian) + rows = cell(1, numel(organisms)); + for i = 1:numel(organisms) + rows{i} = sprintf('%s\tm1\t%s\t%g\t%g\t1\t*\n', ec, organisms{i}, kcatMax, kcatMedian); + end +end + +function writeBrenda(brendaDir, kcatLines) + fid = fopen(fullfile(brendaDir, 'kcat.tsv'), 'w'); + fprintf(fid, '# BRENDA release test-fixture - kcat in 1/s\n'); + fprintf(fid, 'ec_code\tsubstrate\torganism\tkcat_max\tkcat_median\tn\treferences\n'); + for i = 1:numel(kcatLines) + fprintf(fid, kcatLines{i}); + end + fclose(fid); + fid = fopen(fullfile(brendaDir, 'mw.tsv'), 'w'); + fprintf(fid, '# BRENDA release test-fixture - molecular weight in g/mol\n'); + fprintf(fid, 'ec_code\tsubstrate\torganism\tmw\tn\treferences\n'); + fprintf(fid, '0.0.0.0\t*\tplaceholder\t1\t1\t*\n'); + fclose(fid); + fid = fopen(fullfile(brendaDir, 'sa.tsv'), 'w'); + fprintf(fid, '# BRENDA release test-fixture - specific activity in umol/min/mg\n'); + fprintf(fid, 'ec_code\tsubstrate\torganism\tsa_max\tsa_median\tn\treferences\n'); + fprintf(fid, '0.0.0.0\t*\tplaceholder\t1\t1\t1\t*\n'); + fclose(fid); +end + +function adapter = brendaAdapter(testCase, kcatLines) + scratchDir = tempname; + brendaDir = fullfile(scratchDir, 'data'); + mkdir(brendaDir); + testCase.addTeardown(@() rmdir(scratchDir, 's')); + writeBrenda(brendaDir, kcatLines); + adapter = getTestAdapter(); + adapter.params.path = scratchDir; +end + +function [ecModel, fuzzyAdapter] = fuzzyModel(testCase, kcatLines) + fuzzyAdapter = brendaAdapter(testCase, kcatLines); + geckoPath = findGECKOroot; + copyfile(fullfile(geckoPath,'test','unit_tests','ecTestGEM','data','PhylDist.mat'), ... + fullfile(fuzzyAdapter.params.path, 'data', 'PhylDist.mat')); + adapter = getTestAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + ecModel = getECfromGEM(ecModel); +end + +function [ecModel, adapter] = customModel() + adapter = getTestAdapter(); + ecModel = customModelAgain(); +end + +function ecModel = customModelAgain() + adapter = getTestAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + ecModel = getECfromGEM(ecModel); + ecModel.ec.kcat(:) = 10; +end + +function custom = customStruct(proteins, kcats, rxns, notes) + custom.proteins = proteins; + custom.kcat = kcats; + custom.rxns = rxns; + custom.notes = notes; +end From 71b73ce0b1fedce7ab4e3ceafb92eddd7e2ce80d Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 25 Oct 2026 08:35:14 +0100 Subject: [PATCH 07/13] Gene ids that contain 'and', case-insensitive gene rule operators, usage reactions out of the ecFSEOF gene targets, forced wildcards from the right, and usable kcats in light applyKcatConstraints addNewRxnsToEC couples a new reaction to the whole-word genes of its rule. makeEcModel (light) keeps one row for a reaction with an empty gene rule, and gene rule operators in capitals are accepted there and in addNewRxnsToEC. ecFSEOF leaves the usage reactions out before the gene slopes and actions are derived, reports an empty result for a corrThreshold above every correlation, and corrThreshold is tested. fuzzyKcatMatching forceWClvl replaces the rightmost level that is not a wildcard yet, for each EC number of a reaction. applyKcatConstraints (light) lets only isozymes with a positive kcat and a positive enzyme mass set the cost. --- src/change_model/addNewRxnsToEC.m | 8 ++-- src/change_model/applyKcatConstraints.m | 6 ++- src/change_model/makeEcModel.m | 4 ++ src/flux_analysis/ecFSEOF.m | 23 ++++++--- src/gather_kcats/fuzzyKcatMatching.m | 35 ++++++++++++-- test/VerificationMatrix.tsv | 10 ++-- test/unit_tests/ecModelBuildTests.m | 59 ++++++++++++++++++++++++ test/unit_tests/geckoCoreFunctionTests.m | 36 +++++++++++++++ test/unit_tests/kcatMatchTests.m | 13 ++++++ 9 files changed, 175 insertions(+), 19 deletions(-) diff --git a/src/change_model/addNewRxnsToEC.m b/src/change_model/addNewRxnsToEC.m index 7673c93b..fdab1b1b 100644 --- a/src/change_model/addNewRxnsToEC.m +++ b/src/change_model/addNewRxnsToEC.m @@ -140,7 +140,7 @@ % present in the model. genesInRules = strjoin(convertCharArray(newRxns.grRules)); genesInRules = regexp(genesInRules,'[^\s()]+','match'); % Gene ids and the words 'and' and 'or', without spaces and brackets - genesInRules = genesInRules(~ismember(genesInRules,{'and','or'})); % Only the gene ids + genesInRules = genesInRules(~ismember(lower(genesInRules),{'and','or'})); % Only the gene ids genesInRules = setdiff(unique(genesInRules),model.genes); % Only keep new genes % If there is any gene which is not present, as to be pass as @@ -250,8 +250,8 @@ model.ec.eccodes(end+1) = {''}; % Get rules - genes = newRxns.grRules(rxnIdx); - genes = strtrim(erase(split(genes,'and'), {'(' ')'})); + genes = regexp(newRxns.grRules{rxnIdx},'[^\s()]+','match'); % Gene ids and the word 'and', without spaces and brackets + genes = genes(~ismember(lower(genes),{'and','or'})); [~,enzIdx] = ismember(genes, model.ec.genes); % Update the enzyme rxns matrix @@ -269,6 +269,8 @@ %intact. An empty rule gives an empty cell array. function dnf = gprToDnf(rule) tokens = regexp(rule,'\(|\)|[^\s()]+','match'); +isOperator = strcmpi(tokens,'and') | strcmpi(tokens,'or'); +tokens(isOperator) = lower(tokens(isOperator)); % Operators are case insensitive if isempty(tokens) dnf = {}; return diff --git a/src/change_model/applyKcatConstraints.m b/src/change_model/applyKcatConstraints.m index d0536b7c..350ae1e7 100644 --- a/src/change_model/applyKcatConstraints.m +++ b/src/change_model/applyKcatConstraints.m @@ -140,8 +140,12 @@ % isozymes *before* selecting the cheapest one: including them in % the min() below would let an isozyme with no kcat or no enzyme % always win over a real kcat on another isozyme of the same reaction. + % A row with a kcat that is not positive or an enzyme mass that is not + % positive has no usable cost either. hasEnzyme = full(sum(model.ec.rxnEnzMat(ecIdx,:) ~= 0,2)) > 0; - validKcat = ~isinf(MWkcat) & ~isnan(MWkcat) & hasEnzyme; + enzMass = full(model.ec.rxnEnzMat(ecIdx,:) * model.ec.mw); + validKcat = ~isinf(MWkcat) & ~isnan(MWkcat) & hasEnzyme & ... + model.ec.kcat(ecIdx) > 0 & enzMass > 0; if any(validKcat) % Select the lowest MW/kcat (= most efficient), and convert to /hour model.S(prot_pool_idx, hasEc(i)) = -min(MWkcat(validKcat)/3600); diff --git a/src/change_model/makeEcModel.m b/src/change_model/makeEcModel.m index b8737793..313e1c0d 100644 --- a/src/change_model/makeEcModel.m +++ b/src/change_model/makeEcModel.m @@ -229,6 +229,8 @@ %the number of AND-clauses of the GPR in disjunctive normal form, so %'G1 and (G2 or G3)' gives the isozymes {G1,G2} and {G1,G3}. isozymes = cellfun(@gprToDnf, model.grRules(rxnWithGene), 'UniformOutput', false); + % A reaction with gene entries but an empty rule keeps one isozyme without enzymes + isozymes(cellfun(@isempty,isozymes)) = {{cell(1,0)}}; cpys = cellfun(@numel, isozymes); prevNumRxns = numel(isozymes); cpyIndices = repelem(rxnWithGene, cpys).'; %.' only matter when number of rxns is 1 @@ -459,6 +461,8 @@ %intact. An empty rule gives an empty cell array. function dnf = gprToDnf(rule) tokens = regexp(rule,'\(|\)|[^\s()]+','match'); +isOperator = strcmpi(tokens,'and') | strcmpi(tokens,'or'); +tokens(isOperator) = lower(tokens(isOperator)); % Operators are case insensitive if isempty(tokens) dnf = {}; return diff --git a/src/flux_analysis/ecFSEOF.m b/src/flux_analysis/ecFSEOF.m index d4b10c26..df24c4b7 100644 --- a/src/flux_analysis/ecFSEOF.m +++ b/src/flux_analysis/ecFSEOF.m @@ -202,6 +202,15 @@ slope_rxns = slope_rxns(order); target_type = target_type(order); +% Exclude enzyme usage reactions: they are not reported, and they would otherwise +% enter the gene slopes and the choice of the action per gene +notUsage = ~startsWith(target_rxns,'usage_prot_'); +target_rxns = target_rxns(notUsage); +v_matrix = v_matrix(notUsage,:); +rxnGeneM = rxnGeneM(notUsage,:); +slope_rxns = slope_rxns(notUsage); +target_type = target_type(notUsage); + % Create gene list of those connected to the remaining rxns genes = model.genes(sum(rxnGeneM,1) > 0); slope_genes = zeros(size(genes)); @@ -264,20 +273,22 @@ % Report values used to enforce objective target production fseof.alpha = alpha; -% Exclude enzyme usage reactions -toKeep = ~startsWith(target_rxns,'usage_prot_'); -rxnIdx = getIndexes(model,target_rxns(toKeep),'rxns'); +rxnIdx = getIndexes(model,target_rxns,'rxns'); % Report the fluxes at each alpha -fseof.v_matrix = array2table(v_matrix(toKeep,:), ... +fseof.v_matrix = array2table(v_matrix, ... 'VariableNames', string(alpha), 'RowNames', model.rxns(rxnIdx)); % Report reaction targets fseof.rxnTargets(:,1) = model.rxns(rxnIdx); fseof.rxnTargets(:,2) = model.rxnNames(rxnIdx); -fseof.rxnTargets(:,3) = num2cell(slope_rxns(toKeep)); +fseof.rxnTargets(:,3) = num2cell(slope_rxns); fseof.rxnTargets(:,4) = model.grRules(rxnIdx); -fseof.rxnTargets(:,5) = constructEquations(model,rxnIdx); +if isempty(rxnIdx) + fseof.rxnTargets(:,5) = cell(0,1); % constructEquations would give all reactions +else + fseof.rxnTargets(:,5) = constructEquations(model,rxnIdx); +end % Identify transport reactions; defined as reactions involving (at least) % one metabolite name in more than one compartment. diff --git a/src/gather_kcats/fuzzyKcatMatching.m b/src/gather_kcats/fuzzyKcatMatching.m index ed3e1a05..b8aa2017 100644 --- a/src/gather_kcats/fuzzyKcatMatching.m +++ b/src/gather_kcats/fuzzyKcatMatching.m @@ -222,10 +222,7 @@ end %Apply force wildcard level -while forceWClvl > 0 - eccodes=regexprep(eccodes,'(.)*(\.\d+)(\.-)*$','$1\.-$3'); - forceWClvl = forceWClvl - 1; -end +eccodes = forceWildcards(eccodes,forceWClvl); PB = progressReport(mM, 'Gathering kcat values by fuzzy matching to BRENDA database'); %Main loop: @@ -259,6 +256,36 @@ end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%Replaces the rightmost level of each EC code that is not a wildcard yet with +%a wildcard, forceWClvl times, so that level 4 gives -.-.-.- . The codes of one +%entry are separated by ';'. A code that does not have four levels is kept. +function eccodes = forceWildcards(eccodes,forceWClvl) +if forceWClvl <= 0 + return +end +for i = 1:numel(eccodes) + if isempty(eccodes{i}) + continue + end + codes = strsplit(eccodes{i},';'); + for j = 1:numel(codes) + parts = strsplit(codes{j},'.'); + if numel(parts) ~= 4 + continue + end + for k = 1:forceWClvl + last = find(~strcmp(parts,'-'),1,'last'); + if isempty(last) + break + end + parts{last} = '-'; + end + codes{j} = strjoin(parts,'.'); + end + eccodes{i} = strjoin(codes,';'); +end +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [kcat,dir,tot] =iterativeMatch(EC,subs,substrCoeff,i,KCATcell,dir,tot,... name,phylDist,org_index,SAcell,ECIndexIds,EcIndexIndices,kcatAggregation) %Will iteratively try to match the EC number to some registry in BRENDA, diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 6b036e07..4652cbda 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -1,14 +1,14 @@ Function/File (in src) Verification Comments BrendaParser/* Not explicitly tested - The code has existed for a long time change_model/applyComplexData.m tc0003, tc0004, tc0103 -change_model/applyKcatConstraints.m tc0011, tc0215 +change_model/applyKcatConstraints.m tc0011, tc0215, tc0235 change_model/findMetSmiles.m tc0012 change_model/getComplexData.m tc0003, tc0004 - File download not covered by test cases but is tested manually and works. -change_model/makeEcModel.m tc0001, tc0002, tc0210, tc0212 +change_model/makeEcModel.m tc0001, tc0002, tc0210, tc0212, tc0233, tc0234 change_model/setProtPoolSize.m tc0005 change_model/getKcatAcrossIsoenzymes.m tc0011 change_model/applyCustomKcats.m tc0227, tc0228, tc0229, tc0230 -gather_kcats/fuzzyKcatMatching.m tc0010, tc0223, tc0224, tc0225, tc0226 +gather_kcats/fuzzyKcatMatching.m tc0010, tc0223, tc0224, tc0225, tc0226, tc0231 gather_kcats/getStandardKcat.m gather_kcats/mergeDlkcatAndFuzzyKcats.m tc0011 gather_kcats/readDLKcatOutput.m tcm0001, tcm0002 - The test doesn't really check the function, only that some output is read @@ -33,10 +33,10 @@ gather_kcats/extractEnzymeSubstratePairs.m tc0096, tc0015 gather_kcats/convertCustomKcats.m tc0091, tc0092, tc0094 change_model/upgradeEcModel.m tc0093, tc0094, tc0099 io/loadProtData.m tc0101, tc0102 -change_model/addNewRxnsToEC.m tc0211 +change_model/addNewRxnsToEC.m tc0211, tc0232 change_model/mapRxnsToConv.m tc0213 flux_analysis/ecFVA.m tc0214 -flux_analysis/ecFSEOF.m tc0100 +flux_analysis/ecFSEOF.m tc0100, tc0236 utilities/upgradeGECKOproject.m tc0094, tc0099 gather_kcats/readOpenKineticsPredictorOutput.m tc0030, tc0069 gather_kcats/reviewKcatAssignment.m Has only been tested by running it manually, not tested by test cases diff --git a/test/unit_tests/ecModelBuildTests.m b/test/unit_tests/ecModelBuildTests.m index b1e39373..5d31ca9d 100644 --- a/test/unit_tests/ecModelBuildTests.m +++ b/test/unit_tests/ecModelBuildTests.m @@ -92,6 +92,65 @@ function testApplyKcatConstraintsLightSkipsIsozymeWithoutEnzymes_tc0215(testCase end +function testAddNewRxnsToECGeneIdsContainingAnd_tc0232(testCase) + % Gene ids that contain 'and' (Sand1, andA, GandH) are whole words: the + % new reaction is coupled to exactly the enzymes of the genes in its rule. + adapter = getAdapter(); + model = getGeckoTestModel(); + ecModel = makeEcModel(model, false, adapter); + + newRxns.rxns = {'RNEWB'}; + newRxns.rxnNames = {'RNEWB'}; + newRxns.equations = {'m1[c] => e2[e]'}; + newRxns.grRules = {'Sand1 and (andA or GandH)'}; + newEnzymes.enzymes = {'ENEWA';'ENEWB';'ENEWC'}; + newEnzymes.genes = {'Sand1';'andA';'GandH'}; + newEnzymes.mw = [15000;16000;17000]; + + ecModel = addNewRxnsToEC(ecModel, newRxns, newEnzymes, adapter); + verifyEqual(testCase, sort(coupledGenes(ecModel, 'RNEWB_EXP_1')), {'Sand1';'andA'}) + verifyEqual(testCase, sort(coupledGenes(ecModel, 'RNEWB_EXP_2')), {'GandH';'Sand1'}) +end + +function testLightGeneRuleOperatorsAreCaseInsensitive_tc0233(testCase) + % 'AND' and 'OR' are operators like 'and' and 'or'. + adapter = getAdapter(); + model = getGeckoTestModel(); + model.grRules{strcmp(model.rxns, 'R2')} = 'G1 AND (G2 OR G3)'; + ecModel = makeEcModel(model, true, adapter); + verifyEqual(testCase, sort(coupledGenes(ecModel, '001_R2')), {'G1';'G2'}) + verifyEqual(testCase, sort(coupledGenes(ecModel, '002_R2')), {'G1';'G3'}) +end + +function testLightReactionWithEmptyRuleKeepsOneRow_tc0234(testCase) + % A reaction that has gene entries but an empty gene rule keeps one row + % in ec.rxns, without enzymes. + adapter = getAdapter(); + model = getGeckoTestModel(); + model.grRules{strcmp(model.rxns, 'R2')} = ''; + ecModel = makeEcModel(model, true, adapter); + row = strcmp(ecModel.ec.rxns, '001_R2'); + verifyEqual(testCase, sum(row), 1) + verifyEmpty(testCase, coupledGenes(ecModel, '001_R2')) +end + +function testApplyKcatConstraintsLightSkipsNonPositiveKcat_tc0235(testCase) + % A negative kcat is not a usable kcat: the cost of R2 is the one of the + % isozyme with a positive kcat. + adapter = getAdapter(); + model = getGeckoTestModel(); + lecModel = makeEcModel(model, true, adapter); + lecModel.ec.kcat(:) = 10; + first = strcmp(lecModel.ec.rxns, '001_R2'); + second = strcmp(lecModel.ec.rxns, '002_R2'); + lecModel.ec.kcat(first) = -5; + lecModel = applyKcatConstraints(lecModel); + expected = -(lecModel.ec.rxnEnzMat(second,:) * lecModel.ec.mw) / 10 / 3600; + actual = full(lecModel.S(strcmp(lecModel.mets, 'prot_pool'), strcmp(lecModel.rxns, 'R2'))); + verifyGreaterThan(testCase, abs(expected), 0) + verifyEqual(testCase, actual, expected, 'AbsTol', 1e-10) +end + function adapter = getAdapter() geckoPath = findGECKOroot; adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); diff --git a/test/unit_tests/geckoCoreFunctionTests.m b/test/unit_tests/geckoCoreFunctionTests.m index 1ae0c61e..308b830f 100644 --- a/test/unit_tests/geckoCoreFunctionTests.m +++ b/test/unit_tests/geckoCoreFunctionTests.m @@ -2226,4 +2226,40 @@ function testEcFSEOFSelectsTargetsByCorrelation_tc0100(testCase) % The default corrThreshold is 0.9 sameAsDefault = ecFSEOF(ecModel, 'R3', 'S1', 'nSteps', 6, 'modelAdapter', adapter, 'corrThreshold', 0.9); verifyEqual(testCase, sameAsDefault.rxnTargets, fseof.rxnTargets) + + % The correlation of a reaction cannot exceed 1: a threshold above 1 gives no targets + none = ecFSEOF(ecModel, 'R3', 'S1', 'nSteps', 6, 'modelAdapter', adapter, 'corrThreshold', 1.5); + verifyEmpty(testCase, none.rxnTargets) + verifyEmpty(testCase, none.transportTargets) + verifyEmpty(testCase, none.geneTargets) +end + +function testEcFSEOFGeneTargetsComeFromReportedReactions_tc0236(testCase) + % The usage reactions of the enzymes are not reported, and they do not enter + % the gene targets either: the slope of a gene whose rule is the only gene of + % a reported reaction is the slope of that reaction. + geckoPath = findGECKOroot; + adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM', 'TestGEMAdapter.m')); + model = getGeckoTestModel(); + model.geneShortNames = model.genes; + ecModel = makeEcModel(model, false, adapter); + ecModel = getECfromGEM(ecModel); + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R2_EXP_1')) = 50; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R2_EXP_2')) = 10; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R2_REV_EXP_1')) = 50; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R2_REV_EXP_2')) = 10; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R3')) = 20; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R5')) = 100; + ecModel.ec.source(:) = {'manual'}; + ecModel = applyKcatConstraints(ecModel); + ecModel = setProtPoolSize(ecModel, [], [], [], adapter); + adapter.params.bioRxn = 'R5'; + + fseof = ecFSEOF(ecModel, 'R3', 'S1', 'nSteps', 6, 'modelAdapter', adapter); + verifyNotEmpty(testCase, fseof.geneTargets) + for i = 1:size(fseof.geneTargets,1) + rxnRow = strcmp(fseof.rxnTargets(:,4), fseof.geneTargets{i,1}); + verifyEqual(testCase, sum(rxnRow), 1) + verifyEqual(testCase, fseof.geneTargets{i,3}, fseof.rxnTargets{rxnRow,3}, 'AbsTol', 1e-9) + end end diff --git a/test/unit_tests/kcatMatchTests.m b/test/unit_tests/kcatMatchTests.m index abc4471a..891afa27 100644 --- a/test/unit_tests/kcatMatchTests.m +++ b/test/unit_tests/kcatMatchTests.m @@ -98,6 +98,19 @@ function testFuzzyKcatMatchingEmptySubsetSelectsNone_tc0226(testCase) 'modelAdapter', adapter), 'kcatsFromBrenda:emptySubset') end +function testFuzzyKcatMatchingForceWildcardEscalatesFromTheRight_tc0231(testCase) + % forceWClvl replaces the rightmost level that is not a wildcard yet, for + % every EC number of a reaction, also a provisional one: level 4 gives + % -.-.-.- and level 1 gives a wildcard at the last level. + [ecModel, fuzzyAdapter] = fuzzyModel(testCase, {'1.1.1.1\tm1\ttestus testus\t7\t7\t1\t*\n'}); + ecModel.ec.eccodes{strcmp(ecModel.ec.rxns, 'R2_EXP_1')} = '1.1.1.1;3.5.1.n3'; + sel = ismember(ecModel.ec.rxns,{'R2_EXP_1'}); + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', sel, 'modelAdapter', fuzzyAdapter, 'forceWClvl', 1); + verifyEqual(testCase, kcatList.eccodes, {'1.1.1.-;3.5.1.-'}) + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', sel, 'modelAdapter', fuzzyAdapter, 'forceWClvl', 4); + verifyEqual(testCase, kcatList.eccodes, {'-.-.-.-;-.-.-.-'}) +end + function testApplyCustomKcatsSkipsABlankKcat_tc0227(testCase) % An entry without a kcat leaves the reaction as it is. [ecModel, adapter] = customModel(); From 4b9d6d3cdc6b8d4b6b205ffefed34bbe76464546 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 25 Oct 2026 09:12:51 +0100 Subject: [PATCH 08/13] Gene rule operators in any case in addNewRxnsToEC; EC codes without four levels in fuzzyKcatMatching addNewRxnsToEC writes the operators AND and OR of every gene rule in lower case before the reversible reactions are split and the rules are passed to addRxns, so a rule with a single clause such as 'Sand1 AND andA' no longer reaches RAVEN with an operator it does not parse. Gene ids that contain the letters are left as they are. fuzzyKcatMatching stops escalating an EC code that does not have four levels (for example 1.1.1) and returns the no-match values, as geckopy does, instead of indexing a fourth level that is not there. The help text of forceWClvl in fuzzyKcatMatching and kcatsFromBrenda says that it replaces extra levels from the right, on top of any wildcard already in the code. Tests: upper-case operators on a reversible single-clause rule and an OR rule, an isozyme whose enzymes have no mass, trailing wildcards and a level above the remaining levels for forceWClvl, and a code without four levels. --- src/change_model/addNewRxnsToEC.m | 8 +++++- src/gather_kcats/fuzzyKcatMatching.m | 10 +++++-- src/gather_kcats/kcatsFromBrenda.m | 3 +- test/VerificationMatrix.tsv | 6 ++-- test/unit_tests/ecModelBuildTests.m | 43 ++++++++++++++++++++++++++++ test/unit_tests/kcatMatchTests.m | 23 +++++++++++++++ 6 files changed, 86 insertions(+), 7 deletions(-) diff --git a/src/change_model/addNewRxnsToEC.m b/src/change_model/addNewRxnsToEC.m index fdab1b1b..2ba53466 100644 --- a/src/change_model/addNewRxnsToEC.m +++ b/src/change_model/addNewRxnsToEC.m @@ -136,9 +136,15 @@ % 3. Add new reactions + % The operators of a gene rule are case insensitive: write them in lower + % case, as the rest of the function and RAVEN's addRxns expect. Gene ids + % that contain the letters (Sand1, andA, ORF1) are left as they are. + newRxns.grRules = regexprep(convertCharArray(newRxns.grRules), ... + '(?= 4 + if numel(strsplit(EC{k},'.')) ~= 4 + %An EC code without four levels (for example '1.1.1') cannot be + %escalated level by level: stop with the no-match values, as + %geckopy does. + success = true; + elseif wild_num >= 4 %Already fully wildcarded ('-.-.-.-') and still no match -- %stop here instead of computing dot_pos(4-wild_num) with %wild_num==4, an out-of-range index. kcat(k)/origin(k) keep diff --git a/src/gather_kcats/kcatsFromBrenda.m b/src/gather_kcats/kcatsFromBrenda.m index 9ced7d30..5df57419 100644 --- a/src/gather_kcats/kcatsFromBrenda.m +++ b/src/gather_kcats/kcatsFromBrenda.m @@ -19,7 +19,8 @@ % '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). +% number of extra levels, from the right, replaced by a wildcard in the +% EC numbers (on top of any wildcard already in the code; default 0). % ecRxns : logical or double % for which model.ec.rxns kcats are gathered (default: all). A logical % vector or indices of model.ec.rxns, see fuzzyKcatMatching; [] selects diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 4652cbda..61a26a63 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -1,14 +1,14 @@ Function/File (in src) Verification Comments BrendaParser/* Not explicitly tested - The code has existed for a long time change_model/applyComplexData.m tc0003, tc0004, tc0103 -change_model/applyKcatConstraints.m tc0011, tc0215, tc0235 +change_model/applyKcatConstraints.m tc0011, tc0215, tc0235, tc0238 change_model/findMetSmiles.m tc0012 change_model/getComplexData.m tc0003, tc0004 - File download not covered by test cases but is tested manually and works. change_model/makeEcModel.m tc0001, tc0002, tc0210, tc0212, tc0233, tc0234 change_model/setProtPoolSize.m tc0005 change_model/getKcatAcrossIsoenzymes.m tc0011 change_model/applyCustomKcats.m tc0227, tc0228, tc0229, tc0230 -gather_kcats/fuzzyKcatMatching.m tc0010, tc0223, tc0224, tc0225, tc0226, tc0231 +gather_kcats/fuzzyKcatMatching.m tc0010, tc0223, tc0224, tc0225, tc0226, tc0231, tc0239 gather_kcats/getStandardKcat.m gather_kcats/mergeDlkcatAndFuzzyKcats.m tc0011 gather_kcats/readDLKcatOutput.m tcm0001, tcm0002 - The test doesn't really check the function, only that some output is read @@ -33,7 +33,7 @@ gather_kcats/extractEnzymeSubstratePairs.m tc0096, tc0015 gather_kcats/convertCustomKcats.m tc0091, tc0092, tc0094 change_model/upgradeEcModel.m tc0093, tc0094, tc0099 io/loadProtData.m tc0101, tc0102 -change_model/addNewRxnsToEC.m tc0211, tc0232 +change_model/addNewRxnsToEC.m tc0211, tc0232, tc0237 change_model/mapRxnsToConv.m tc0213 flux_analysis/ecFVA.m tc0214 flux_analysis/ecFSEOF.m tc0100, tc0236 diff --git a/test/unit_tests/ecModelBuildTests.m b/test/unit_tests/ecModelBuildTests.m index 5d31ca9d..523ecd7b 100644 --- a/test/unit_tests/ecModelBuildTests.m +++ b/test/unit_tests/ecModelBuildTests.m @@ -112,6 +112,30 @@ function testAddNewRxnsToECGeneIdsContainingAnd_tc0232(testCase) verifyEqual(testCase, sort(coupledGenes(ecModel, 'RNEWB_EXP_2')), {'GandH';'Sand1'}) end +function testAddNewRxnsToECUpperCaseOperators_tc0237(testCase) + % 'AND' and 'OR' are operators like 'and' and 'or', also in a rule with a + % single clause and on a reversible reaction, where the rule is copied as + % it is to the reverse reaction. + adapter = getAdapter(); + model = getGeckoTestModel(); + ecModel = makeEcModel(model, false, adapter); + + newRxns.rxns = {'RNEWC';'RNEWD'}; + newRxns.rxnNames = {'RNEWC';'RNEWD'}; + newRxns.equations = {'m1[c] <=> e2[e]';'m1[c] => e2[e]'}; + newRxns.grRules = {'Sand1 AND andA';'Sand1 OR andA'}; + newEnzymes.enzymes = {'ENEWA';'ENEWB'}; + newEnzymes.genes = {'Sand1';'andA'}; + newEnzymes.mw = [15000;16000]; + + ecModel = addNewRxnsToEC(ecModel, newRxns, newEnzymes, adapter); + verifyEqual(testCase, ecModel.grRules(strcmp(ecModel.rxns, 'RNEWC')), {'Sand1 and andA'}) + verifyEqual(testCase, sort(coupledGenes(ecModel, 'RNEWC')), {'Sand1';'andA'}) + verifyEqual(testCase, sort(coupledGenes(ecModel, 'RNEWC_REV')), {'Sand1';'andA'}) + verifyEqual(testCase, coupledGenes(ecModel, 'RNEWD_EXP_1'), {'Sand1'}) + verifyEqual(testCase, coupledGenes(ecModel, 'RNEWD_EXP_2'), {'andA'}) +end + function testLightGeneRuleOperatorsAreCaseInsensitive_tc0233(testCase) % 'AND' and 'OR' are operators like 'and' and 'or'. adapter = getAdapter(); @@ -151,6 +175,25 @@ function testApplyKcatConstraintsLightSkipsNonPositiveKcat_tc0235(testCase) verifyEqual(testCase, actual, expected, 'AbsTol', 1e-10) end +function testApplyKcatConstraintsLightSkipsEnzymeWithoutMass_tc0238(testCase) + % An isozyme whose enzymes have no mass gives a cost of 0, which is not a + % usable cost: the cost of R2 is the one of the isozyme with a mass. + adapter = getAdapter(); + model = getGeckoTestModel(); + lecModel = makeEcModel(model, true, adapter); + lecModel.ec.kcat(:) = 10; + first = strcmp(lecModel.ec.rxns, '001_R2'); + second = strcmp(lecModel.ec.rxns, '002_R2'); + onlyFirst = find(full(lecModel.ec.rxnEnzMat(first,:)) > 0 & full(lecModel.ec.rxnEnzMat(second,:)) == 0); + verifyNotEmpty(testCase, onlyFirst) + lecModel.ec.mw(onlyFirst) = 0; + lecModel = applyKcatConstraints(lecModel); + expected = -full(lecModel.ec.rxnEnzMat(second,:) * lecModel.ec.mw) / 10 / 3600; + actual = full(lecModel.S(strcmp(lecModel.mets, 'prot_pool'), strcmp(lecModel.rxns, 'R2'))); + verifyGreaterThan(testCase, abs(expected), 0) + verifyEqual(testCase, actual, expected, 'AbsTol', 1e-10) +end + function adapter = getAdapter() geckoPath = findGECKOroot; adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); diff --git a/test/unit_tests/kcatMatchTests.m b/test/unit_tests/kcatMatchTests.m index 891afa27..3ec73a03 100644 --- a/test/unit_tests/kcatMatchTests.m +++ b/test/unit_tests/kcatMatchTests.m @@ -109,6 +109,29 @@ function testFuzzyKcatMatchingForceWildcardEscalatesFromTheRight_tc0231(testCase verifyEqual(testCase, kcatList.eccodes, {'1.1.1.-;3.5.1.-'}) kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', sel, 'modelAdapter', fuzzyAdapter, 'forceWClvl', 4); verifyEqual(testCase, kcatList.eccodes, {'-.-.-.-;-.-.-.-'}) + + % A code with a trailing wildcard escalates from the last level that is not + % a wildcard yet, and a level above the number of levels left stops at -.-.-.- + ecModel.ec.eccodes{strcmp(ecModel.ec.rxns, 'R2_EXP_1')} = '1.1.-.-;3.5.1.n3'; + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', sel, 'modelAdapter', fuzzyAdapter, 'forceWClvl', 1); + verifyEqual(testCase, kcatList.eccodes, {'1.-.-.-;3.5.1.-'}) + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', sel, 'modelAdapter', fuzzyAdapter, 'forceWClvl', 6); + verifyEqual(testCase, kcatList.eccodes, {'-.-.-.-;-.-.-.-'}) +end + +function testFuzzyKcatMatchingEcCodeWithoutFourLevels_tc0239(testCase) + % An EC code without four levels is not escalated level by level: it is + % left as it is by forceWClvl and, when BRENDA has no match, ends with the + % no-match values instead of an error. + [ecModel, fuzzyAdapter] = fuzzyModel(testCase, {'1.1.1.1\tm1\ttestus testus\t7\t7\t1\t*\n'}); + ecModel.ec.eccodes{strcmp(ecModel.ec.rxns, 'R2_EXP_1')} = '9.9.9'; + sel = ismember(ecModel.ec.rxns,{'R2_EXP_1'}); + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', sel, 'modelAdapter', fuzzyAdapter); + verifyEqual(testCase, kcatList.eccodes, {'9.9.9'}) + verifyEqual(testCase, kcatList.kcats, 0) + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', sel, 'modelAdapter', fuzzyAdapter, 'forceWClvl', 1); + verifyEqual(testCase, kcatList.eccodes, {'9.9.9'}) + verifyEqual(testCase, kcatList.kcats, 0) end function testApplyCustomKcatsSkipsABlankKcat_tc0227(testCase) From a2bf2317995b806d9d443a6560f0de594bd864cf Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 25 Oct 2026 09:16:11 +0100 Subject: [PATCH 09/13] mergeKcats and mergeDLKcatAndFuzzyKcats default to topOriginLimit 4; findECInDB drops a wildcard code next to a specific code for a single gene Origin 5 and 6 matches come from specific activities times molecular weights and do not always beat a prediction, so they rank after DLKcat by default; topOriginLimit 6 ranks them above, as the gecko_classic preset does. A gene whose database cell holds '1.1.1.4 1.1.1.- 1.1.1.303' now gives 1.1.1.4;1.1.1.303, as an enzyme complex already did. --- src/gather_kcats/mergeDLKcatAndFuzzyKcats.m | 5 +- src/gather_kcats/mergeKcats.m | 7 ++- src/get_enzyme_data/findECInDB.m | 15 ++++-- test/VerificationMatrix.tsv | 5 +- test/unit_tests/kcatAssignmentTests.m | 3 +- test/unit_tests/parityDecisionTests.m | 58 +++++++++++++++++++++ 6 files changed, 82 insertions(+), 11 deletions(-) create mode 100644 test/unit_tests/parityDecisionTests.m diff --git a/src/gather_kcats/mergeDLKcatAndFuzzyKcats.m b/src/gather_kcats/mergeDLKcatAndFuzzyKcats.m index 4bb1e0dd..364efeac 100644 --- a/src/gather_kcats/mergeDLKcatAndFuzzyKcats.m +++ b/src/gather_kcats/mergeDLKcatAndFuzzyKcats.m @@ -28,7 +28,8 @@ % -------------------- % topOriginLimit : double % origin limit for prioritized BRENDA matches. Origin is explained in -% more detail below (default 6). +% more detail below (default 4: BRENDA matches of origin 5 and 6, from +% specific activities, rank after DLKcat; use 6 to rank them above). % bottomOriginLimit : double % origin limit for low priority BRENDA matches. Origin is explained in % more detail below (default 6). @@ -74,7 +75,7 @@ % mergeKcats, readDLKcatOutput, fuzzyKcatMatching, assignKcatValues p = parseGECKOargs(varargin, { ... - 'topOriginLimit', 6; ... + 'topOriginLimit', 4; ... 'bottomOriginLimit', 6; ... 'wildcardLimit', 3}); topOriginLimit = p.topOriginLimit; diff --git a/src/gather_kcats/mergeKcats.m b/src/gather_kcats/mergeKcats.m index ad0f65c5..483c51f3 100644 --- a/src/gather_kcats/mergeKcats.m +++ b/src/gather_kcats/mergeKcats.m @@ -33,7 +33,10 @@ % Name-Value Arguments % -------------------- % topOriginLimit : double -% origin limit for database_top (default 6). +% origin limit for database_top (default 4). Origins 5 and 6 are specific +% activity times molecular weight kcats, which do not always beat a +% prediction; they rank after the other sources and are used where +% nothing else exists, or above DLKcat with topOriginLimit 6. % bottomOriginLimit : double % origin limit for database_bottom (default 6). % wildcardLimit : double @@ -87,7 +90,7 @@ % mergeDLKcatAndFuzzyKcats, readDLKcatOutput, fuzzyKcatMatching, assignKcatValues p = parseGECKOargs(varargin, { ... - 'topOriginLimit', 6; ... + 'topOriginLimit', 4; ... 'bottomOriginLimit', 6; ... 'wildcardLimit', 3; ... 'databaseSources', {'brenda','sabio_rk'}}); diff --git a/src/get_enzyme_data/findECInDB.m b/src/get_enzyme_data/findECInDB.m index f23c8dbb..fa59305b 100644 --- a/src/get_enzyme_data/findECInDB.m +++ b/src/get_enzyme_data/findECInDB.m @@ -4,10 +4,12 @@ % For each gene, finds its matching database entries and takes the EC number % of the lightest-weight matching protein when a single EC number applies to % that gene; genes with several distinct EC numbers among their matches are -% also reported as conflicts (the first match is still used for EC). When -% gene_set has more than one gene (an enzyme complex), the numbers found for -% each gene are combined by intersection where they share a match (wildcards -% such as "-" are accounted for), otherwise by union. +% also reported as conflicts (the first match is still used for EC). A +% wildcard code such as "1.1.1.-" is dropped when a more specific code of +% the same prefix ("1.1.1.4") is present. When gene_set has more than one +% gene (an enzyme complex), the numbers found for each gene are combined by +% intersection where they share a match (wildcards such as "-" are +% accounted for), otherwise by union. % % Parameters % ---------- @@ -102,6 +104,11 @@ int_EC = intersection(int_EC,other_EC); end end +%A single gene keeps no wildcard code next to a more specific code of the +%same prefix (e.g. 'EC1.1.1.-' next to 'EC1.1.1.4'), as in a complex: +if isscalar(EC_set) + int_EC = compare_wild(int_EC); +end %Use the intersection found, if any. If not, use the union: if isempty(int_EC) EC_set = uni_EC; diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 61a26a63..6df447cd 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -10,7 +10,7 @@ change_model/getKcatAcrossIsoenzymes.m tc0011 change_model/applyCustomKcats.m tc0227, tc0228, tc0229, tc0230 gather_kcats/fuzzyKcatMatching.m tc0010, tc0223, tc0224, tc0225, tc0226, tc0231, tc0239 gather_kcats/getStandardKcat.m -gather_kcats/mergeDlkcatAndFuzzyKcats.m tc0011 +gather_kcats/mergeDlkcatAndFuzzyKcats.m tc0011; tc0241 gather_kcats/readDLKcatOutput.m tcm0001, tcm0002 - The test doesn't really check the function, only that some output is read 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; also tc0104 @@ -41,8 +41,9 @@ utilities/upgradeGECKOproject.m tc0094, tc0099 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 +gather_kcats/mergeKcats.m tc0240 get_enzyme_data/calculateMW.m Not explicitly tested - The code has existed for a long time -get_enzyme_data/findECInDB.m Called from getECfromDatabase - tested when that function is tested. +get_enzyme_data/findECInDB.m Called from getECfromDatabase - tested when that function is tested; also tc0242 get_enzyme_data/getECfromDatabase.m tc0007 - does not test download of the databases - this is tested in the manual workflows for Yeast-GEM and Human-GEM get_enzyme_data/getECfromGEM.m tc0006 get_enzyme_data/getECstring.m Called from findECInDB - tested when that function is tested. diff --git a/test/unit_tests/kcatAssignmentTests.m b/test/unit_tests/kcatAssignmentTests.m index 4f612967..68312c88 100644 --- a/test/unit_tests/kcatAssignmentTests.m +++ b/test/unit_tests/kcatAssignmentTests.m @@ -14,7 +14,8 @@ function testAssignKcatsReproducesChain_tc0079(testCase) 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'}); + % gecko_classic ranks the specific-activity matches (origin 5 and 6) above DLKcat + merged = mergeKcats({fuzzy, dlkcat}, {'database_top','dlkcat','database_bottom'}, 'topOriginLimit', 6); A = assignKcatValues(base, merged); A = applyCustomKcats(A, 'modelAdapter', adapter); A = getKcatAcrossIsozymes(A); diff --git a/test/unit_tests/parityDecisionTests.m b/test/unit_tests/parityDecisionTests.m new file mode 100644 index 00000000..e9f97dd5 --- /dev/null +++ b/test/unit_tests/parityDecisionTests.m @@ -0,0 +1,58 @@ +%run this test case with the command +%results = runtests('parityDecisionTests.m') + +function tests = parityDecisionTests + tests = functiontests(localfunctions); +end + +function testMergeKcatsDefaultTopOriginLimitIs4_tc0240(testCase) + % A BRENDA match of origin 5 (specific activity, correct organism) ranks + % after DLKcat by default and above it with topOriginLimit 6. + [fuzzy, dlkcat] = originFiveLists(); + tiers = {'database_top', 'dlkcat', 'database_bottom'}; + merged = mergeKcats({fuzzy, dlkcat}, tiers); + verifyEqual(testCase, merged.kcats(strcmp(merged.rxns, 'R1')), 50) + merged = mergeKcats({fuzzy, dlkcat}, tiers, 'topOriginLimit', 4); + verifyEqual(testCase, merged.kcats(strcmp(merged.rxns, 'R1')), 50) + merged = mergeKcats({fuzzy, dlkcat}, tiers, 'topOriginLimit', 6); + verifyEqual(testCase, merged.kcats(strcmp(merged.rxns, 'R1')), 5) +end + +function testMergeDLKcatAndFuzzyKcatsDefaultTopOriginLimitIs4_tc0241(testCase) + [fuzzy, dlkcat] = originFiveLists(); + merged = mergeDLKcatAndFuzzyKcats(dlkcat, fuzzy); + verifyEqual(testCase, merged.kcats(strcmp(merged.rxns, 'R1')), 50) + merged = mergeDLKcatAndFuzzyKcats(dlkcat, fuzzy, 'topOriginLimit', 6); + verifyEqual(testCase, merged.kcats(strcmp(merged.rxns, 'R1')), 5) +end + +function testFindECInDBSingleGeneDropsWildcardParent_tc0242(testCase) + % One gene whose protein carries '1.1.1.4 1.1.1.- 1.1.1.303' keeps the + % two specific codes and loses the wildcard of the same prefix. + geneHashMap = containers.Map({'G1'}, {1}); + geneIndex = {1}; + DBecNum = {'1.1.1.4 1.1.1.- 1.1.1.303'}; + DBMW = 10000; + EC = findECInDB({'G1'}, DBecNum, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.4;1.1.1.303') + % A wildcard without a specific code of its prefix stays + DBecNum = {'1.1.1.- 2.2.2.2'}; + EC = findECInDB({'G1'}, DBecNum, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.-;2.2.2.2') +end + +function [fuzzy, dlkcat] = originFiveLists() + fuzzy.rxns = {'R1'}; + fuzzy.kcats = 5; + fuzzy.genes = {{}}; + fuzzy.substrates = {{}}; + fuzzy.eccodes = {'1.1.1.1'}; + fuzzy.wildcardLvl = 0; + fuzzy.origin = 5; + fuzzy.source = 'brenda'; + dlkcat.rxns = {'R1'}; + dlkcat.kcats = 50; + dlkcat.genes = {'G1'}; + dlkcat.substrates = {'m1'}; + dlkcat.source = 'DLKcat'; +end From 9465dd5b57f2522ffe798913d0b89e0f7d85a832 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 25 Oct 2026 09:36:37 +0100 Subject: [PATCH 10/13] applyKcatConstraints clears the constraints of a full ecModel when every kcat is 0, sensitivityTuning stops when growth gains less than 1e-9, getKcatAcrossIsozymes averages the isozymes with a kcat above 0 --- src/change_model/applyKcatConstraints.m | 13 +++- src/change_model/getKcatAcrossIsozymes.m | 9 ++- src/kcat_tuning/sensitivityTuning.m | 12 ++- test/VerificationMatrix.tsv | 6 +- test/unit_tests/parityRound3Tests.m | 94 ++++++++++++++++++++++++ 5 files changed, 121 insertions(+), 13 deletions(-) create mode 100644 test/unit_tests/parityRound3Tests.m diff --git a/src/change_model/applyKcatConstraints.m b/src/change_model/applyKcatConstraints.m index 350ae1e7..e87f7ce2 100644 --- a/src/change_model/applyKcatConstraints.m +++ b/src/change_model/applyKcatConstraints.m @@ -3,7 +3,9 @@ % % Existing enzyme constraints are first removed (unless updateRxns is % provided), and new constraints are defined based on the content of -% model.ec.kcat. +% model.ec.kcat. When every kcat is 0, the existing constraints of a full +% ecModel are removed and no new ones are written; a light ecModel is +% returned unchanged. % % Parameters % ---------- @@ -66,8 +68,13 @@ ' not a valid GECKO3 ecModel. First run makeEcModel(model).']) end if all(model.ec.kcat==0) - printOrange('WARNING: No kcat values are provided in model.ec.kcat, model remains unchanged.\n'); - return + if model.ec.geckoLight + printOrange('WARNING: No kcat values are provided in model.ec.kcat, model remains unchanged.\n'); + return + end + % In a full model the existing enzyme constraints of the selected + % reactions are still cleared below, and no new ones are written. + printOrange('WARNING: No kcat values are provided in model.ec.kcat, existing enzyme constraints are cleared and no new ones are written.\n'); end %Clear existing incorporation of enzyme usage diff --git a/src/change_model/getKcatAcrossIsozymes.m b/src/change_model/getKcatAcrossIsozymes.m index 47c03664..71dfe485 100644 --- a/src/change_model/getKcatAcrossIsozymes.m +++ b/src/change_model/getKcatAcrossIsozymes.m @@ -4,7 +4,9 @@ % For reactions without kcat value (0 in model.ec.kcat), isozymes are found % (being based on the same reaction in the conventional GEM), that do have a % kcat value assigned. The mean kcat value of these isozymes is then used to -% fill in model.ec.kcat. +% fill in model.ec.kcat. Only isozymes with a kcat above 0 enter the mean: a +% NaN or negative kcat of an isozyme is left out, and an entry whose kcat is +% NaN is not filled. % % Parameters % ---------- @@ -34,10 +36,11 @@ end noKcats = model.ec.kcat==0; +hasKcat = model.ec.kcat>0; % NaN and negative kcats are not isozyme values rxnIDs = regexprep(model.ec.rxns,'_EXP_\d+',''); noKcatID = rxnIDs(noKcats); -yesKcatID = rxnIDs(~noKcats); -yesKcatVal = model.ec.kcat(~noKcats); +yesKcatID = rxnIDs(hasKcat); +yesKcatVal = model.ec.kcat(hasKcat); noKcatVal = cellfun(@(x) strcmp(x, yesKcatID), noKcatID, 'UniformOutput', false); noKcatVal = cell2mat(cellfun(@(x) mean(yesKcatVal(x)), noKcatVal, 'UniformOutput', false)); diff --git a/src/kcat_tuning/sensitivityTuning.m b/src/kcat_tuning/sensitivityTuning.m index 9e0facc2..600f6531 100644 --- a/src/kcat_tuning/sensitivityTuning.m +++ b/src/kcat_tuning/sensitivityTuning.m @@ -3,7 +3,8 @@ % % Iteratively increases the kcat of whichever enzyme currently limits growth % the most, until the model reaches the desired growth rate, updating the -% tuned values in model.ec.kcat. +% tuned values in model.ec.kcat. The loop also stops, with a warning, when an +% iteration raises growth by less than 1e-9. % % Parameters % ---------- @@ -87,7 +88,10 @@ if isempty(res.x) error('FBA of input model gives no valid result. Reduce protein pool constraint with setProtPoolSize and/or check if exchange constraints are correctly defined.') end -lastGrowth = 0; +% The loop stops when growth rises by less than growthStallTol over the +% previous iteration (the first iteration has nothing to compare with). +lastGrowth = NaN; +growthStallTol = 1e-9; if ~m.ec.geckoLight %for the full model, we first find the draw reaction with the most flux drawRxns = startsWith(m.rxns, 'usage_prot_'); @@ -95,7 +99,7 @@ iteration = 1; while true [res,hs] = solveLP(m,0,[],hs); %skip parsimonius, only takes time - if (lastGrowth == res.f) + if ~isnan(lastGrowth) && res.f <= lastGrowth + growthStallTol printOrange('WARNING: No growth increase from increased kcats - check if the constraints on the uptake reactions are too tight.\n'); break; end @@ -146,7 +150,7 @@ iteration = 1; while true res = solveLP(m,0); %skip parsimonius, only takes time - if (lastGrowth == res.f) + if ~isnan(lastGrowth) && res.f <= lastGrowth + growthStallTol printOrange('No growth increase from increased kcats - check if the constraints on the uptake reactions are too tight.\n'); break; end diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 6df447cd..35e78574 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -1,12 +1,12 @@ Function/File (in src) Verification Comments BrendaParser/* Not explicitly tested - The code has existed for a long time change_model/applyComplexData.m tc0003, tc0004, tc0103 -change_model/applyKcatConstraints.m tc0011, tc0215, tc0235, tc0238 +change_model/applyKcatConstraints.m tc0011, tc0215, tc0235, tc0238, tc0250 change_model/findMetSmiles.m tc0012 change_model/getComplexData.m tc0003, tc0004 - File download not covered by test cases but is tested manually and works. change_model/makeEcModel.m tc0001, tc0002, tc0210, tc0212, tc0233, tc0234 change_model/setProtPoolSize.m tc0005 -change_model/getKcatAcrossIsoenzymes.m tc0011 +change_model/getKcatAcrossIsoenzymes.m tc0011, tc0252 change_model/applyCustomKcats.m tc0227, tc0228, tc0229, tc0230 gather_kcats/fuzzyKcatMatching.m tc0010, tc0223, tc0224, tc0225, tc0226, tc0231, tc0239 gather_kcats/getStandardKcat.m @@ -53,7 +53,7 @@ download_databases/loadDatabases.m Called from getECfromDatabase - tested when t kcat_tuning/evotune/* Has only been tested by running it manually, not tested by test cases kcat_tuning/fillCarbonNum.m tc0055 (deprecated-alias check only, not a real usage test) kcat_tuning/findTopLimitations.m Not explicitly tested - The code has existed for a long time -kcat_tuning/sensitivityTuning.m tcm0001 +kcat_tuning/sensitivityTuning.m tcm0001, tc0251 kcat_tuning/sigmaFitter.m Not explicitly tested - The code has existed for a long time kcat_tuning/topUsedEnzymes.m Not explicitly tested - The code has existed for a long time Not sure this works anymore - written for the old model structure? kcat_tuning/truncateValues.m Not explicitly tested - The code has existed for a long time diff --git a/test/unit_tests/parityRound3Tests.m b/test/unit_tests/parityRound3Tests.m new file mode 100644 index 00000000..82bf394f --- /dev/null +++ b/test/unit_tests/parityRound3Tests.m @@ -0,0 +1,94 @@ +%run this test case with the command +%results = runtests('parityRound3Tests.m') + +function tests = parityRound3Tests + tests = functiontests(localfunctions); +end + +function testApplyKcatConstraintsAllZeroClearsFullModel_tc0250(testCase) + % With every kcat 0 the enzyme constraints of a full ecModel are removed + % and no new ones are written. A light ecModel is returned unchanged. + adapter = getAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + ecModel.ec.kcat(:) = 10; + ecModel = applyKcatConstraints(ecModel); + protMets = startsWith(ecModel.mets, 'prot_') & ~strcmp(ecModel.mets, 'prot_pool'); + metabolic = ~startsWith(ecModel.rxns, 'usage_prot_') & ~strcmp(ecModel.rxns, 'prot_pool_exchange'); + verifyGreaterThan(testCase, nnz(ecModel.S(protMets, metabolic)), 0) + ecModel.ec.kcat(:) = 0; + cleared = applyKcatConstraints(ecModel); + verifyEqual(testCase, nnz(cleared.S(protMets, metabolic)), 0) + % the usage reactions and the pool exchange keep their fixed entries + verifyEqual(testCase, cleared.S(:, ~metabolic), ecModel.S(:, ~metabolic)) +end + +function testApplyKcatConstraintsAllZeroLightIsUnchanged_tc0250light(testCase) + % A light ecModel with every kcat 0 is returned unchanged (control: the + % behaviour is the same before and after the full-model change). + adapter = getAdapter(); + lecModel = makeEcModel(getGeckoTestModel(), true, adapter); + lecModel.ec.kcat(:) = 10; + lecModel = applyKcatConstraints(lecModel); + poolRow = strcmp(lecModel.mets, 'prot_pool'); + verifyGreaterThan(testCase, nnz(lecModel.S(poolRow, :)), 1) + lecModel.ec.kcat(:) = 0; + same = applyKcatConstraints(lecModel); + verifyEqual(testCase, same.S, lecModel.S) +end + +function testSensitivityTuningStopsOnGrowthBelowTolerance_tc0251(testCase) + % R3 limits growth and its enzyme P4 is ignored, so tuning only raises + % R5. Growth stops rising once R3 alone limits it: the loop ends on the + % first iteration whose growth gain is below 1e-9, instead of raising + % the kcat of R5 until growth is exactly unchanged. + adapter = getAdapter(); + model = loadConventionalGEM('modelAdapter', adapter); + for id = {'R2', 'R4'} + idx = strcmp(model.rxns, id{1}); + model.lb(idx) = 0; model.ub(idx) = 0; + end + ecModel = makeEcModel(model, false, adapter); + ecModel = getECfromGEM(ecModel); + kcats = struct('R3', 0.05, 'R5', 200); + for id = fieldnames(kcats)' + j = strcmp(ecModel.ec.rxns, id{1}); + ecModel.ec.kcat(j) = kcats.(id{1}); + ecModel.ec.source{j} = 'manual'; + end + ecModel = applyKcatConstraints(ecModel); + ecModel = setProtPoolSize(ecModel, [], [], [], adapter); + tuneAdapter = adapter; + tuneAdapter.params.bioRxn = 'R5'; + tunedModel = sensitivityTuning(ecModel, 'desiredGrowthRate', 50, ... + 'modelAdapter', tuneAdapter, 'protToIgnore', {'P4'}, 'verbose', false); + r5 = tunedModel.ec.kcat(strcmp(tunedModel.ec.rxns, 'R5')); + verifyGreaterThan(testCase, r5, 200) + verifyLessThan(testCase, r5, 1e12) +end + +function testGetKcatAcrossIsozymesUsesPositiveSiblingsOnly_tc0252(testCase) + % Only isozymes with a kcat above 0 enter the mean: a negative or NaN + % sibling does not fill a kcat of 0. + adapter = getAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + ecModel.ec.kcat(:) = 0; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R3')) = 8; + % R2: a zero next to a negative sibling is not filled + ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R2_EXP_2')) = -3; + % R2_REV: a zero next to a NaN sibling is not filled + ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R2_REV_EXP_2')) = NaN; + out = getKcatAcrossIsozymes(ecModel); + verifyEqual(testCase, out.ec.kcat(strcmp(out.ec.rxns, 'R2_EXP_1')), 0) + verifyEqual(testCase, out.ec.kcat(strcmp(out.ec.rxns, 'R2_REV_EXP_1')), 0) + % a positive sibling fills the zero (control: the same before and after) + ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R2_EXP_2')) = 6; + out = getKcatAcrossIsozymes(ecModel); + j = strcmp(out.ec.rxns, 'R2_EXP_1'); + verifyEqual(testCase, out.ec.kcat(j), 6) + verifyEqual(testCase, out.ec.source{j}, 'isozymes') +end + +function adapter = getAdapter() +geckoPath = findGECKOroot; +adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); +end From 25fcbf47fb6b8b92af242cd7accf963d33516447 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 25 Oct 2026 09:43:52 +0100 Subject: [PATCH 11/13] sensitivityTuning ends with a warning when no protein usage reaction carries flux --- src/kcat_tuning/sensitivityTuning.m | 8 ++++++++ test/VerificationMatrix.tsv | 2 +- test/unit_tests/parityRound3Tests.m | 23 +++++++++++++++++++++++ 3 files changed, 32 insertions(+), 1 deletion(-) diff --git a/src/kcat_tuning/sensitivityTuning.m b/src/kcat_tuning/sensitivityTuning.m index 600f6531..dba95750 100644 --- a/src/kcat_tuning/sensitivityTuning.m +++ b/src/kcat_tuning/sensitivityTuning.m @@ -117,6 +117,10 @@ drawFluxes(drawRxns) = res.x(drawRxns); % Remove from the list user defined proteins drawFluxes(idxToIgnore) = 0; + if ~any(drawFluxes > 0) + printOrange('WARNING: No protein usage reaction carries flux, tuning cannot continue.\n'); + break; + end [~,sel] = max(drawFluxes); % since bounds 0 to 1000 %Now get the metabolite metSel = m.S(:,sel) > 0; % positive coeff (forward usage) @@ -166,6 +170,10 @@ %find the highest protein usage flux protPoolStoich = m.S(strcmp(m.mets, 'prot_pool'),:).'; protPoolStoich(idxToIgnore) = 0; + if ~any(res.x .* protPoolStoich < 0) + printOrange('WARNING: No reaction draws on the protein pool, tuning cannot continue.\n'); + break; + end [~,sel] = min(res.x .* protPoolStoich); %max consumption kcatList = [kcatList, sel]; rxn = m.rxns(sel.'); diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 35e78574..3686ea47 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -53,7 +53,7 @@ download_databases/loadDatabases.m Called from getECfromDatabase - tested when t kcat_tuning/evotune/* Has only been tested by running it manually, not tested by test cases kcat_tuning/fillCarbonNum.m tc0055 (deprecated-alias check only, not a real usage test) kcat_tuning/findTopLimitations.m Not explicitly tested - The code has existed for a long time -kcat_tuning/sensitivityTuning.m tcm0001, tc0251 +kcat_tuning/sensitivityTuning.m tcm0001, tc0251, tc0253 kcat_tuning/sigmaFitter.m Not explicitly tested - The code has existed for a long time kcat_tuning/topUsedEnzymes.m Not explicitly tested - The code has existed for a long time Not sure this works anymore - written for the old model structure? kcat_tuning/truncateValues.m Not explicitly tested - The code has existed for a long time diff --git a/test/unit_tests/parityRound3Tests.m b/test/unit_tests/parityRound3Tests.m index 82bf394f..fc99e96d 100644 --- a/test/unit_tests/parityRound3Tests.m +++ b/test/unit_tests/parityRound3Tests.m @@ -66,6 +66,29 @@ function testSensitivityTuningStopsOnGrowthBelowTolerance_tc0251(testCase) verifyLessThan(testCase, r5, 1e12) end +function testSensitivityTuningStopsWhenNothingCarriesFlux_tc0253(testCase) + % R5, the growth reaction, is blocked: growth is 0 from the start and no + % usage reaction carries flux. The loop ends with a warning instead of + % tuning an arbitrary reaction or failing, and no kcat changes. + adapter = getAdapter(); + model = loadConventionalGEM('modelAdapter', adapter); + for id = {'R2', 'R4', 'R5'} + idx = strcmp(model.rxns, id{1}); + model.lb(idx) = 0; model.ub(idx) = 0; + end + ecModel = makeEcModel(model, false, adapter); + ecModel = getECfromGEM(ecModel); + ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R3')) = 0.05; + ecModel = applyKcatConstraints(ecModel); + ecModel = setProtPoolSize(ecModel, [], [], [], adapter); + tuneAdapter = adapter; + tuneAdapter.params.bioRxn = 'R5'; + [tunedModel, tunedKcats] = sensitivityTuning(ecModel, 'desiredGrowthRate', 50, ... + 'modelAdapter', tuneAdapter, 'verbose', false); + verifyEqual(testCase, tunedModel.ec.kcat, ecModel.ec.kcat) + verifyEmpty(testCase, tunedKcats.rxns) +end + function testGetKcatAcrossIsozymesUsesPositiveSiblingsOnly_tc0252(testCase) % Only isozymes with a kcat above 0 enter the mean: a negative or NaN % sibling does not fill a kcat of 0. From 8e2f63da5e5fac46b523262195d98956c93da8c1 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 25 Oct 2026 16:46:13 +0100 Subject: [PATCH 12/13] sensitivityTuning ends with a warning when the LP has no solution; findECInDB leaves codes with fewer than four levels alone sensitivityTuning stops with a warning when solveLP returns no objective value inside the tuning loop (full and light branches), as geckopy does for a solve that is not optimal. For a single gene, findECInDB removes a wildcard code only next to a four-level code of the same prefix and removes empty codes, so '1.1.1 1.1.1.4' keeps both codes and a trailing space leaves no empty code. The help text states that the combined EC set of a complex is not cleaned of a wildcard next to a specific code. --- src/get_enzyme_data/findECInDB.m | 69 +++++++++++++-- src/kcat_tuning/sensitivityTuning.m | 8 ++ test/VerificationMatrix.tsv | 4 +- test/unit_tests/parityDecisionTests.m | 12 +++ test/unit_tests/parityRound3Tests.m | 123 +++++++++++++++++++++----- 5 files changed, 183 insertions(+), 33 deletions(-) diff --git a/src/get_enzyme_data/findECInDB.m b/src/get_enzyme_data/findECInDB.m index fa59305b..85277fec 100644 --- a/src/get_enzyme_data/findECInDB.m +++ b/src/get_enzyme_data/findECInDB.m @@ -4,12 +4,14 @@ % For each gene, finds its matching database entries and takes the EC number % of the lightest-weight matching protein when a single EC number applies to % that gene; genes with several distinct EC numbers among their matches are -% also reported as conflicts (the first match is still used for EC). A -% wildcard code such as "1.1.1.-" is dropped when a more specific code of -% the same prefix ("1.1.1.4") is present. When gene_set has more than one -% gene (an enzyme complex), the numbers found for each gene are combined by -% intersection where they share a match (wildcards such as "-" are -% accounted for), otherwise by union. +% also reported as conflicts (the first match is still used for EC). For a +% single gene, a wildcard code such as "1.1.1.-" is dropped when a more +% specific code of the same prefix ("1.1.1.4") is present; only codes with +% four levels take part, and empty codes are removed. When gene_set has more +% than one gene (an enzyme complex), the numbers found for each gene are +% combined by intersection where they share a match (wildcards such as "-" +% are accounted for), otherwise by union, and the combined set is not +% cleaned of a wildcard next to a specific code. % % Parameters % ---------- @@ -105,9 +107,9 @@ end end %A single gene keeps no wildcard code next to a more specific code of the -%same prefix (e.g. 'EC1.1.1.-' next to 'EC1.1.1.4'), as in a complex: +%same prefix (e.g. 'EC1.1.1.-' next to 'EC1.1.1.4'): if isscalar(EC_set) - int_EC = compare_wild(int_EC); + int_EC = dropWildcardParents(int_EC); end %Use the intersection found, if any. If not, use the union: if isempty(int_EC) @@ -153,6 +155,57 @@ int_EC = unique(int_EC,'stable'); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function int_EC = dropWildcardParents(EC) +%Removes empty codes and every code that is a wildcard parent of another +%code in the cell array (e.g. 'EC1.2.3.-' when 'EC1.2.3.4' is present). Only +%codes with four levels are compared, so 'EC1.1.1' is neither a parent nor a +%child of 'EC1.1.1.4'. Equal codes collapse to the first. +EC = EC(~cellfun(@isempty,EC) & ~strcmp(EC,'EC')); +keep = true(1,numel(EC)); +for i = 1:numel(EC) + if ~keep(i) + continue; + end + for j = i+1:numel(EC) + if ~keep(j) + continue; + end + iInJ = subsumes(EC{i},EC{j}); + jInI = subsumes(EC{j},EC{i}); + if iInJ && jInI + keep(j) = false; + elseif iInJ + keep(i) = false; + break; + elseif jInI + keep(j) = false; + end + end +end +int_EC = EC(keep); +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function tf = subsumes(general,specific) +%True if the code general is equal to or a wildcard parent of the code +%specific; both carry the "EC" prefix and need four levels. +general = strsplit(extractAfter(general,2),'.'); +specific = strsplit(extractAfter(specific,2),'.'); +tf = false; +if numel(general) ~= 4 || numel(specific) ~= 4 + return; +end +for k = 1:4 + if strcmp(general{k},'-') + tf = true; + return; + end + if strcmp(specific{k},'-') || ~strcmp(general{k},specific{k}) + return; + end +end +tf = true; +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function int_EC = compare_wild(EC) %Goes through a cell array of EC numbers, and erases any repetitions, %considering also wildcards (e.g. will erase 'EC1.2.3.-' if 'EC1.2.3.4' is diff --git a/src/kcat_tuning/sensitivityTuning.m b/src/kcat_tuning/sensitivityTuning.m index dba95750..96d1efce 100644 --- a/src/kcat_tuning/sensitivityTuning.m +++ b/src/kcat_tuning/sensitivityTuning.m @@ -99,6 +99,10 @@ iteration = 1; while true [res,hs] = solveLP(m,0,[],hs); %skip parsimonius, only takes time + if isempty(res.f) + printOrange('WARNING: The model has no feasible solution with the increased kcats, tuning stops.\n'); + break; + end if ~isnan(lastGrowth) && res.f <= lastGrowth + growthStallTol printOrange('WARNING: No growth increase from increased kcats - check if the constraints on the uptake reactions are too tight.\n'); break; @@ -154,6 +158,10 @@ iteration = 1; while true res = solveLP(m,0); %skip parsimonius, only takes time + if isempty(res.f) + printOrange('WARNING: The model has no feasible solution with the increased kcats, tuning stops.\n'); + break; + end if ~isnan(lastGrowth) && res.f <= lastGrowth + growthStallTol printOrange('No growth increase from increased kcats - check if the constraints on the uptake reactions are too tight.\n'); break; diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 3686ea47..024a5070 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -43,7 +43,7 @@ gather_kcats/reviewKcatAssignment.m Has only been tested by running it manually, gather_kcats/ecStatsFromBrenda.m Not verified against real BRENDA data in this environment -- see its own TODO comment gather_kcats/mergeKcats.m tc0240 get_enzyme_data/calculateMW.m Not explicitly tested - The code has existed for a long time -get_enzyme_data/findECInDB.m Called from getECfromDatabase - tested when that function is tested; also tc0242 +get_enzyme_data/findECInDB.m Called from getECfromDatabase - tested when that function is tested; also tc0242, tc0243 get_enzyme_data/getECfromDatabase.m tc0007 - does not test download of the databases - this is tested in the manual workflows for Yeast-GEM and Human-GEM get_enzyme_data/getECfromGEM.m tc0006 get_enzyme_data/getECstring.m Called from findECInDB - tested when that function is tested. @@ -53,7 +53,7 @@ download_databases/loadDatabases.m Called from getECfromDatabase - tested when t kcat_tuning/evotune/* Has only been tested by running it manually, not tested by test cases kcat_tuning/fillCarbonNum.m tc0055 (deprecated-alias check only, not a real usage test) kcat_tuning/findTopLimitations.m Not explicitly tested - The code has existed for a long time -kcat_tuning/sensitivityTuning.m tcm0001, tc0251, tc0253 +kcat_tuning/sensitivityTuning.m tcm0001, tc0251, tc0253, tc0254 kcat_tuning/sigmaFitter.m Not explicitly tested - The code has existed for a long time kcat_tuning/topUsedEnzymes.m Not explicitly tested - The code has existed for a long time Not sure this works anymore - written for the old model structure? kcat_tuning/truncateValues.m Not explicitly tested - The code has existed for a long time diff --git a/test/unit_tests/parityDecisionTests.m b/test/unit_tests/parityDecisionTests.m index e9f97dd5..de3f7277 100644 --- a/test/unit_tests/parityDecisionTests.m +++ b/test/unit_tests/parityDecisionTests.m @@ -41,6 +41,18 @@ function testFindECInDBSingleGeneDropsWildcardParent_tc0242(testCase) verifyEqual(testCase, EC, '1.1.1.-;2.2.2.2') end +function testFindECInDBSingleGeneShortCodesAndTrailingSpace_tc0243(testCase) + % A code with fewer than four levels is neither the wildcard parent nor + % the child of a four-level code, and a trailing space leaves no empty code. + geneHashMap = containers.Map({'G1'}, {1}); + geneIndex = {1}; + DBMW = 10000; + EC = findECInDB({'G1'}, {'1.1.1 1.1.1.4'}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1;1.1.1.4') + EC = findECInDB({'G1'}, {'1.1.1.4 '}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.4') +end + function [fuzzy, dlkcat] = originFiveLists() fuzzy.rxns = {'R1'}; fuzzy.kcats = 5; diff --git a/test/unit_tests/parityRound3Tests.m b/test/unit_tests/parityRound3Tests.m index fc99e96d..d230faf9 100644 --- a/test/unit_tests/parityRound3Tests.m +++ b/test/unit_tests/parityRound3Tests.m @@ -37,33 +37,31 @@ function testApplyKcatConstraintsAllZeroLightIsUnchanged_tc0250light(testCase) end function testSensitivityTuningStopsOnGrowthBelowTolerance_tc0251(testCase) - % R3 limits growth and its enzyme P4 is ignored, so tuning only raises - % R5. Growth stops rising once R3 alone limits it: the loop ends on the - % first iteration whose growth gain is below 1e-9, instead of raising - % the kcat of R5 until growth is exactly unchanged. - adapter = getAdapter(); - model = loadConventionalGEM('modelAdapter', adapter); - for id = {'R2', 'R4'} - idx = strcmp(model.rxns, id{1}); - model.lb(idx) = 0; model.ub(idx) = 0; - end - ecModel = makeEcModel(model, false, adapter); - ecModel = getECfromGEM(ecModel); - kcats = struct('R3', 0.05, 'R5', 200); - for id = fieldnames(kcats)' - j = strcmp(ecModel.ec.rxns, id{1}); - ecModel.ec.kcat(j) = kcats.(id{1}); - ecModel.ec.source{j} = 'manual'; - end - ecModel = applyKcatConstraints(ecModel); - ecModel = setProtPoolSize(ecModel, [], [], [], adapter); - tuneAdapter = adapter; - tuneAdapter.params.bioRxn = 'R5'; + % The solver reports a growth of 1 in the first iteration and 1 + 5e-10 in + % the second: the gain is below 1e-9, so the loop stops after the second + % iteration (3 solves with the one before the loop). An exact-equality + % loop would go on until two growths are identical (4 solves). + [ecModel, tuneAdapter] = tuningModel(); + cleanup = shadowSolveLP('growth', [1 1 1+5e-10]); %#ok tunedModel = sensitivityTuning(ecModel, 'desiredGrowthRate', 50, ... 'modelAdapter', tuneAdapter, 'protToIgnore', {'P4'}, 'verbose', false); + verifyEqual(testCase, getappdata(0, 'gecko_test_solvelp_count'), 3) + r5 = tunedModel.ec.kcat(strcmp(tunedModel.ec.rxns, 'R5')); + verifyGreaterThan(testCase, r5, 200) +end + +function testSensitivityTuningEndsWhenTheLPBecomesInfeasible_tc0254(testCase) + % The third solve (the second iteration) has no solution: RAVEN returns + % an empty objective value. The loop ends with a warning and returns + % what was tuned in the first iteration. + [ecModel, tuneAdapter] = tuningModel(); + cleanup = shadowSolveLP('infeasibleFrom', 3); %#ok + out = evalc(['tunedModel = sensitivityTuning(ecModel, ''desiredGrowthRate'', 50, ' ... + '''modelAdapter'', tuneAdapter, ''protToIgnore'', {''P4''}, ''verbose'', false);']); + verifyTrue(testCase, contains(out, 'no feasible solution')) + verifyEqual(testCase, getappdata(0, 'gecko_test_solvelp_count'), 3) r5 = tunedModel.ec.kcat(strcmp(tunedModel.ec.rxns, 'R5')); verifyGreaterThan(testCase, r5, 200) - verifyLessThan(testCase, r5, 1e12) end function testSensitivityTuningStopsWhenNothingCarriesFlux_tc0253(testCase) @@ -109,6 +107,85 @@ function testGetKcatAcrossIsozymesUsesPositiveSiblingsOnly_tc0252(testCase) j = strcmp(out.ec.rxns, 'R2_EXP_1'); verifyEqual(testCase, out.ec.kcat(j), 6) verifyEqual(testCase, out.ec.source{j}, 'isozymes') + % a NaN sibling next to a positive one is left out of the mean: three + % isozymes of R2 with kcats 0, NaN and 6 give 6 to the first + ecModel.ec.kcat(:) = 0; + ecModel.ec.rxns{strcmp(ecModel.ec.rxns, 'R3')} = 'R2_EXP_3'; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R2_EXP_2')) = NaN; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R2_EXP_3')) = 6; + out = getKcatAcrossIsozymes(ecModel); + j = strcmp(out.ec.rxns, 'R2_EXP_1'); + verifyEqual(testCase, out.ec.kcat(j), 6) + verifyEqual(testCase, out.ec.source{j}, 'isozymes') + verifyTrue(testCase, isnan(out.ec.kcat(strcmp(out.ec.rxns, 'R2_EXP_2')))) +end + +function [ecModel, tuneAdapter] = tuningModel() +% R3 limits growth, its enzyme P4 is ignored by the tests, so tuning raises R5. +adapter = getAdapter(); +model = loadConventionalGEM('modelAdapter', adapter); +for id = {'R2', 'R4'} + idx = strcmp(model.rxns, id{1}); + model.lb(idx) = 0; model.ub(idx) = 0; +end +ecModel = makeEcModel(model, false, adapter); +ecModel = getECfromGEM(ecModel); +kcats = struct('R3', 0.05, 'R5', 200); +for id = fieldnames(kcats)' + j = strcmp(ecModel.ec.rxns, id{1}); + ecModel.ec.kcat(j) = kcats.(id{1}); + ecModel.ec.source{j} = 'manual'; +end +ecModel = applyKcatConstraints(ecModel); +ecModel = setProtPoolSize(ecModel, [], [], [], adapter); +tuneAdapter = adapter; +tuneAdapter.params.bioRxn = 'R5'; +end + +function cleanup = shadowSolveLP(mode, values) +% Puts a solveLP on the path that calls RAVEN's and then reports the growth +% of solve number n as values(min(n, end)) ('growth'), or no solution from +% solve number values on ('infeasibleFrom'). The number of solves is kept in +% the root appdata gecko_test_solvelp_count. The path is restored when the +% returned object is cleared. +setappdata(0, 'gecko_test_solvelp_real', @solveLP); +setappdata(0, 'gecko_test_solvelp_count', 0); +setappdata(0, 'gecko_test_solvelp_mode', mode); +setappdata(0, 'gecko_test_solvelp_values', values); +stubDir = tempname; +mkdir(stubDir); +fid = fopen(fullfile(stubDir, 'solveLP.m'), 'w'); +fprintf(fid, '%s\n', ... + 'function [res,hs] = solveLP(varargin)', ... + 'n = getappdata(0,''gecko_test_solvelp_count'') + 1;', ... + 'setappdata(0,''gecko_test_solvelp_count'',n);', ... + 'real = getappdata(0,''gecko_test_solvelp_real'');', ... + 'mode = getappdata(0,''gecko_test_solvelp_mode'');', ... + 'values = getappdata(0,''gecko_test_solvelp_values'');', ... + '[res,hs] = real(varargin{:});', ... + 'switch mode', ... + ' case ''growth''', ... + ' res.f = values(min(n,numel(values)));', ... + ' case ''infeasibleFrom''', ... + ' if n >= values', ... + ' res.f = [];', ... + ' res.x = [];', ... + ' end', ... + 'end', ... + 'end'); +fclose(fid); +addpath(stubDir, '-begin'); +cleanup = onCleanup(@() removeShadow(stubDir)); +end + +function removeShadow(stubDir) +rmpath(stubDir); +rmdir(stubDir, 's'); +for name = {'real', 'mode', 'values'} + if isappdata(0, ['gecko_test_solvelp_' name{1}]) + rmappdata(0, ['gecko_test_solvelp_' name{1}]); + end +end end function adapter = getAdapter() From 0b6060fd55e23172e004a492d9b769561ba79c15 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 25 Oct 2026 17:29:03 +0100 Subject: [PATCH 13/13] findECInDB cleans the EC set of an enzyme complex of a wildcard next to a specific code, as the set of a single gene, and leaves out genes without an EC number The cross-gene union and intersection use the same four-level wildcard rule as the single-gene clean-up, replacing the trimmed-substring comparison, and the combined set is cleaned before the codes are joined. --- src/get_enzyme_data/findECInDB.m | 122 ++++++++++---------------- test/VerificationMatrix.tsv | 2 +- test/unit_tests/parityDecisionTests.m | 36 ++++++++ 3 files changed, 82 insertions(+), 78 deletions(-) diff --git a/src/get_enzyme_data/findECInDB.m b/src/get_enzyme_data/findECInDB.m index 85277fec..50abdd05 100644 --- a/src/get_enzyme_data/findECInDB.m +++ b/src/get_enzyme_data/findECInDB.m @@ -4,14 +4,14 @@ % For each gene, finds its matching database entries and takes the EC number % of the lightest-weight matching protein when a single EC number applies to % that gene; genes with several distinct EC numbers among their matches are -% also reported as conflicts (the first match is still used for EC). For a -% single gene, a wildcard code such as "1.1.1.-" is dropped when a more -% specific code of the same prefix ("1.1.1.4") is present; only codes with -% four levels take part, and empty codes are removed. When gene_set has more -% than one gene (an enzyme complex), the numbers found for each gene are -% combined by intersection where they share a match (wildcards such as "-" -% are accounted for), otherwise by union, and the combined set is not -% cleaned of a wildcard next to a specific code. +% also reported as conflicts (the first match is still used for EC). Genes +% without an EC number are left out. When gene_set has more than one gene +% (an enzyme complex), the numbers found for each gene are combined by +% intersection where they share a match (a wildcard such as "1.1.1.-" matches +% the specific codes of its prefix), otherwise by union. The EC set of every +% gene set, one gene or a complex, is cleaned of a wildcard code such as +% "1.1.1.-" when a more specific code of the same prefix ("1.1.1.4") is +% present; only codes with four levels take part, and empty codes are removed. % % Parameters % ---------- @@ -92,30 +92,33 @@ end end -%EC: Find union and intersection between all units (only applies for -%complexes, i.e. length(EC_set) > 1): -uni_EC = strsplit(EC_set{1},' '); -int_EC = uni_EC; -for j = 2:length(EC_set) - other_EC = strsplit(EC_set{j},' '); - if isempty(uni_EC) - uni_EC = other_EC; - int_EC = other_EC; - elseif ~isempty(other_EC) - uni_EC = compare_wild([uni_EC other_EC]); - int_EC = intersection(int_EC,other_EC); - end -end -%A single gene keeps no wildcard code next to a more specific code of the -%same prefix (e.g. 'EC1.1.1.-' next to 'EC1.1.1.4'): -if isscalar(EC_set) - int_EC = dropWildcardParents(int_EC); +%EC: Split the codes of each gene, leave out the genes without any, and +%combine the rest (union and intersection only apply for complexes, i.e. more +%than one gene with codes): +tokens = cell(size(EC_set)); +for j = 1:length(EC_set) + codes = strsplit(EC_set{j},' '); + tokens{j} = codes(~cellfun(@isempty,codes)); end -%Use the intersection found, if any. If not, use the union: -if isempty(int_EC) - EC_set = uni_EC; +tokens = tokens(~cellfun(@isempty,tokens)); +if isempty(tokens) + EC_set = {}; +elseif isscalar(tokens) + EC_set = dropWildcardParents(tokens{1}); else - EC_set = int_EC; + uni_EC = tokens{1}; + int_EC = tokens{1}; + for j = 2:length(tokens) + uni_EC = dropWildcardParents([uni_EC tokens{j}]); + int_EC = wildcardIntersection(int_EC,tokens{j}); + end + %Use the intersection found, if any. If not, use the union: + if isempty(int_EC) + EC_set = uni_EC; + else + EC_set = int_EC; + end + EC_set = dropWildcardParents(EC_set); end %Add new codes as new possible isozymes and remove the initial "EC" text: @@ -127,7 +130,11 @@ end %EC_set = cellfun(@(x) extractAfter(x,2),EC_set,'UniformOutput',false);%remove 2 first chars -EC = strjoin(EC_set,';'); +if isempty(EC_set) + EC = ''; +else + EC = strjoin(EC_set,';'); +end %Delete repeated Uniprots (and the corresponding ECs): @@ -136,23 +143,21 @@ end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -function int_EC = intersection(prev_EC,new_EC) +function int_EC = wildcardIntersection(prev_EC,new_EC) %Finds the common elements between two cell arrays, if any. Also considers %wildcards (e.g. if 'EC1.1.1.1' is in one array and 'EC1.1.1.-' is in the -%other one, then 'EC1.1.1.1' is added to the intersection). +%other one, then 'EC1.1.1.1' is added to the intersection). A pair adds its +%code only when one code of the pair is the wildcard parent of the other, or +%the two are equal four-level codes. int_EC = {}; for i = 1:length(prev_EC) for j = 1:length(new_EC) - new_int = compare_wild({prev_EC{i} new_EC{j}}); - if length(new_int) == 1 - int_EC = [int_EC new_int]; + survivors = dropWildcardParents({prev_EC{i} new_EC{j}}); + if length(survivors) == 1 + int_EC = [int_EC survivors]; end end end -%When prev_EC carries both a code and its own wildcard parent (e.g. -%'EC1.1.1.1' and 'EC1.1.1.-'), both independently pair-match the same -%new_EC entry above and the loop appends the resolved code twice; dedupe. -int_EC = unique(int_EC,'stable'); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function int_EC = dropWildcardParents(EC) @@ -205,40 +210,3 @@ end tf = true; end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -function int_EC = compare_wild(EC) -%Goes through a cell array of EC numbers, and erases any repetitions, -%considering also wildcards (e.g. will erase 'EC1.2.3.-' if 'EC1.2.3.4' is -%already present). - -%Trim all EC numbers of wild cards (e.g. 'EC1.2.-.-' -> 'EC1.2.'): -EC_trimmed = EC; -for i = 1:length(EC) - %Add a final dot for avoiding issues (e.g. 'EC1.1.1.1' & 'EC1.1.1.12'): - ECi = [EC{i} '.']; - pos = strfind(ECi,'-'); - if ~isempty(pos) - ECi = ECi(1:pos(1)-1); - end - EC_trimmed{i} = ECi; -end - -%Compare all EC numbers between them to find if 1 fits in the other: -non_repeated = true(1,length(EC)); -for i = 1:length(EC)-1 - for j = i+1:length(EC) - ECi = EC_trimmed{i}; - ECj = EC_trimmed{j}; - %If ECj fits in ECi then ECj can be disregarded: - if strfind(ECi,ECj) - non_repeated(j) = false; - %Else, if ECi fits in ECj then ECi can be disregarded: - elseif strfind(ECj,ECi) - non_repeated(i) = false; - end - end -end - -int_EC = EC(non_repeated); - -end diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 024a5070..d14c6743 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -43,7 +43,7 @@ gather_kcats/reviewKcatAssignment.m Has only been tested by running it manually, gather_kcats/ecStatsFromBrenda.m Not verified against real BRENDA data in this environment -- see its own TODO comment gather_kcats/mergeKcats.m tc0240 get_enzyme_data/calculateMW.m Not explicitly tested - The code has existed for a long time -get_enzyme_data/findECInDB.m Called from getECfromDatabase - tested when that function is tested; also tc0242, tc0243 +get_enzyme_data/findECInDB.m Called from getECfromDatabase - tested when that function is tested; also tc0242, tc0243, tc0244, tc0245 get_enzyme_data/getECfromDatabase.m tc0007 - does not test download of the databases - this is tested in the manual workflows for Yeast-GEM and Human-GEM get_enzyme_data/getECfromGEM.m tc0006 get_enzyme_data/getECstring.m Called from findECInDB - tested when that function is tested. diff --git a/test/unit_tests/parityDecisionTests.m b/test/unit_tests/parityDecisionTests.m index de3f7277..0c11fee0 100644 --- a/test/unit_tests/parityDecisionTests.m +++ b/test/unit_tests/parityDecisionTests.m @@ -53,6 +53,42 @@ function testFindECInDBSingleGeneShortCodesAndTrailingSpace_tc0243(testCase) verifyEqual(testCase, EC, '1.1.1.4') end +function testFindECInDBComplexDropsWildcardParent_tc0244(testCase) + % The shared EC set of an enzyme complex loses a wildcard code next to a + % specific code of its prefix, as the set of a single gene does. + DBMW = [10000; 20000; 30000]; + geneHashMap = containers.Map({'G1','G2','G3'}, {1,2,3}); + geneIndex = {1; 2; 3}; + EC = findECInDB({'G1','G2'}, {'1.1.1.4 1.1.1.-'; '1.1.1.4 1.1.1.-'; ''}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.4') + % A wildcard that only the union carries is dropped too + EC = findECInDB({'G1','G2','G3'}, {'1.1.1.4'; '1.1.1.- 2.2.2.2'; '5.5.5.5'}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.4;2.2.2.2;5.5.5.5') + % A wildcard shared by every gene stays, as nothing more specific exists + EC = findECInDB({'G1','G2'}, {'1.1.1.-'; '1.1.1.-'; ''}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.-') +end + +function testFindECInDBGeneWithoutECIsLeftOut_tc0245(testCase) + % A gene without an EC number does not change the codes of the others: + % the complex is treated as the genes that have codes, in either order. + DBMW = [10000; 20000]; + geneHashMap = containers.Map({'G1','G2'}, {1,2}); + geneIndex = {1; 2}; + EC = findECInDB({'G1','G2'}, {'1.1.1.4 1.1.1.-'; ''}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.4') + EC = findECInDB({'G1','G2'}, {''; '1.1.1.4 1.1.1.-'}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.4') + % No gene has a code + EC = findECInDB({'G1','G2'}, {''; ''}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '') + % A gene missing from the database is left out as well + geneHashMap = containers.Map({'G2'}, {2}); + geneIndex = {[]; 2}; + EC = findECInDB({'G1','G2'}, {''; '1.1.1.4 1.1.1.-'}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.4') +end + function [fuzzy, dlkcat] = originFiveLists() fuzzy.rxns = {'R1'}; fuzzy.kcats = 5;