From 708774181a1c990d66e2680a5feab8e9d5f44ec5 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 00:32:32 +0100 Subject: [PATCH 1/2] sigmaFitter: return the start of the growth plateau on every solver Growth that is flat from some sigma on is returned with rounding noise by some solvers. The smallest sigma whose growth is within 1e-9 (relative) of the best growth is now chosen instead of the exact minimum error. --- src/kcat_tuning/sigmaFitter.m | 9 +++++- test/VerificationMatrix.tsv | 2 +- test/unit_tests/parityRound3Tests.m | 45 +++++++++++++++++++++++++++++ 3 files changed, 54 insertions(+), 2 deletions(-) diff --git a/src/kcat_tuning/sigmaFitter.m b/src/kcat_tuning/sigmaFitter.m index 4bff81bc..9ffbfe0c 100644 --- a/src/kcat_tuning/sigmaFitter.m +++ b/src/kcat_tuning/sigmaFitter.m @@ -95,7 +95,14 @@ error = num2str(((growthRate-solution.x(objPos))/growthRate)*100); sigParam = [sigParam; sigma]; end -[~, minIndx] = min(errors); +% Growth that is flat from some sigma on (a target above what the model can +% reach) is returned by the solver with rounding noise. The smallest sigma +% whose growth is within 1e-9 (relative) of the best growth is chosen, so the +% result is the start of the plateau on every solver. The tolerance is +% expressed on the percent error. +[minError, minIndx] = min(errors); +errTol = 1e-9*100*abs(objValues(minIndx))/abs(growthRate); +minIndx = find(errors <= minError + errTol, 1); sigma = sigParam(minIndx); % The loop above leaves the protein pool sized for the *last* sigma tried % (i=100, i.e. sigma=1), not the best-fitting one found by the search -- diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index f673afb8..8ab10779 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -54,7 +54,7 @@ kcat_tuning/evotune/* Has only been tested by running it manually, not tested by 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, tc0254 -kcat_tuning/sigmaFitter.m Not explicitly tested - The code has existed for a long time +kcat_tuning/sigmaFitter.m tc0020, tc0255 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 limit_proteins/constrainProtConcs.m tc0013 diff --git a/test/unit_tests/parityRound3Tests.m b/test/unit_tests/parityRound3Tests.m index d230faf9..28191148 100644 --- a/test/unit_tests/parityRound3Tests.m +++ b/test/unit_tests/parityRound3Tests.m @@ -120,6 +120,51 @@ function testGetKcatAcrossIsozymesUsesPositiveSiblingsOnly_tc0252(testCase) verifyTrue(testCase, isnan(out.ec.kcat(strcmp(out.ec.rxns, 'R2_EXP_2')))) end +function testSigmaFitterReturnsStartOfGrowthPlateau_tc0255(testCase) +% The growth target is above what the uptake-limited model reaches, so growth +% is flat from sigma 0.06 on. GLPK returns the plateau growth with rounding +% noise (49.999999999999986 at sigma 0.06, 50 at sigma 1). The fitted sigma is +% the start of the plateau, 0.06, whatever the noise, and the same on Gurobi. +previous = getpref('RAVEN', 'solver', 'glpk'); +restore = onCleanup(@() setRavenSolver(previous)); %#ok +setRavenSolver('glpk'); +sigmaGlpk = plateauSigma(); +verifyEqual(testCase, sigmaGlpk, 0.06, 'AbsTol', 1e-12) +hasGurobi = false; +if exist('gurobi', 'file') + try + setRavenSolver('gurobi'); + hasGurobi = true; + catch + end +end +assumeTrue(testCase, hasGurobi, 'Gurobi is not available') +verifyEqual(testCase, plateauSigma(), sigmaGlpk, 'AbsTol', 1e-12) +end + +function sigma = plateauSigma() +adapter = getAdapter(); +model = loadConventionalGEM('modelAdapter', adapter); +model.lb(strcmp(model.rxns, 'S1')) = -50; +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', 5, '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.c = double(strcmp(ecModel.rxns, 'R5')); +ecModel = applyKcatConstraints(ecModel); +ecModel = setProtPoolSize(ecModel, [], [], [], adapter); +[~, sigma] = sigmaFitter(ecModel, 'growthRate', 100, 'Ptot', 0.5, 'f', 4, ... + 'makePlot', false, 'modelAdapter', adapter); +end + function [ecModel, tuneAdapter] = tuningModel() % R3 limits growth, its enzyme P4 is ignored by the tests, so tuning raises R5. adapter = getAdapter(); From 793593cb2c29d33c5fb0e7845382c9ffc7687c51 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 08:47:41 +0100 Subject: [PATCH 2/2] getStandardKcat, getSubsetEcModel: handle subSystem case and empty subSystems, models without genes to remove, light ecModels getStandardKcat treats subSystem names that differ only in letter case as one subSystem, as its lookup already did, instead of raising when two list entries match. A reaction with an empty subSystem gets the model-wide standard kcat instead of raising an index error. getSubsetEcModel returns the model when no gene has to be removed, instead of raising in removeGenes. In a light ecModel the reaction identifiers of the ec structure carry the prefix 001_ and were all compared to the model reactions with it, so every entry was removed. The prefix is now removed before the comparison. --- src/change_model/getSubsetEcModel.m | 14 +++++- src/gather_kcats/getStandardKcat.m | 59 +++++++--------------- test/VerificationMatrix.tsv | 3 +- test/unit_tests/parityRound3Tests.m | 77 +++++++++++++++++++++++++++++ 4 files changed, 109 insertions(+), 44 deletions(-) diff --git a/src/change_model/getSubsetEcModel.m b/src/change_model/getSubsetEcModel.m index d7f484ed..ae09990c 100644 --- a/src/change_model/getSubsetEcModel.m +++ b/src/change_model/getSubsetEcModel.m @@ -50,7 +50,11 @@ if ismember('standard', genesToRemove) genesToRemove(strcmp(genesToRemove, 'standard')) = []; end -smallEcModel = removeGenes(bigEcModel,genesToRemove,true,true,false); +if isempty(genesToRemove) + smallEcModel = bigEcModel; +else + smallEcModel = removeGenes(bigEcModel,genesToRemove,true,true,false); +end % Remove genes from ec-structure enzToRemove = ismember(smallEcModel.ec.genes,genesToRemove); @@ -68,7 +72,13 @@ smallEcModel = removeReactions(smallEcModel,~keepRxns, true, true, true); % Remove reactions from ec-structure -ecRxnsToRemove = ~ismember(smallEcModel.ec.rxns,smallEcModel.rxns); +% The reaction identifiers of a light ecModel carry a 4 character prefix +if smallEcModel.ec.geckoLight + ecModelRxns = extractAfter(smallEcModel.ec.rxns,4); +else + ecModelRxns = smallEcModel.ec.rxns; +end +ecRxnsToRemove = ~ismember(ecModelRxns,smallEcModel.rxns); smallEcModel.ec.rxns(ecRxnsToRemove) = []; smallEcModel.ec.kcat(ecRxnsToRemove) = []; diff --git a/src/gather_kcats/getStandardKcat.m b/src/gather_kcats/getStandardKcat.m index 7ebe2345..6cad6083 100644 --- a/src/gather_kcats/getStandardKcat.m +++ b/src/gather_kcats/getStandardKcat.m @@ -133,13 +133,14 @@ end % Map ec-rxns to model.rxns [~,rxnIdx] = ismember(modRxns,model.rxns); - % Choose first subSystem - enzSubSystems = flattenCell(model.subSystems(rxnIdx)); - enzSubSystems = enzSubSystems(:,1); + % Choose first subSystem ('' when a reaction has none) + enzSubSystems = cellfun(@firstSubSystem, model.subSystems(rxnIdx), 'UniformOutput', false); if ~all(cellfun(@isempty, enzSubSystems)) - % Make list of unique subsystems, and which rxns are linked to them - [enzSubSystem_names, ~, rxnToSub] = unique(enzSubSystems); + % Make list of unique subsystems, and which rxns are linked to them. + % Names that differ only in letter case are one subSystem. + [~, firstOccurrence, rxnToSub] = unique(lower(enzSubSystems)); + enzSubSystem_names = enzSubSystems(firstOccurrence); % Calculate the median kcat per subSystem, counting only reactions % with a real (nonzero, non-NaN) kcat toward both the median itself % and the threshold that decides whether a subSystem-specific value @@ -289,8 +290,9 @@ % all(), is the right test -- all() only ever succeeds when the model % has exactly one subSystem in total, silently falling back to the % model-wide standardKcat for every reaction otherwise. - kcatSubSystemIdx = strcmpi(enzSubSystem_names, model.subSystems{rxnIdx}(1)); - if any(kcatSubSystemIdx) + rxnSubSystem = firstSubSystem(model.subSystems{rxnIdx}); + kcatSubSystemIdx = strcmpi(enzSubSystem_names, rxnSubSystem); + if ~isempty(rxnSubSystem) && any(kcatSubSystemIdx) model.ec.kcat(end+1) = kcatSubSystem(kcatSubSystemIdx); else model.ec.kcat(end+1) = standardKcat; @@ -319,39 +321,14 @@ end end -function Cflat = flattenCell(C,strFlag) -%FLATTENCELL Flatten a nested column cell array into a matrix cell array. -% -% CFLAT = FLATTENCELL(C) takes a column cell array in which one or more -% entries is a nested cell array, and flattens it into a 2D matrix cell -% array, where the nested entries are spread across new columns. -% -% CFLAT = FLATTENCELL(C,STRFLAG) if STRFLAG is TRUE, empty entries in the -% resulting CFLAT will be replaced with empty strings {''}. Default = FALSE -if nargin < 2 - strFlag = false; -end - -% determine which entries are cells -cells = cellfun(@iscell,C); - -% determine number of elements in each nested cell -cellsizes = cellfun(@numel,C); -cellsizes(~cells) = 1; % ignore non-cell entries - -% flatten single-entry cells -Cflat = C; -Cflat(cells & (cellsizes == 1)) = cellfun(@(x) x{1},Cflat(cells & (cellsizes == 1)),'UniformOutput',false); - -% iterate through multi-entry cells -multiCells = find(cellsizes > 1); -for i = 1:length(multiCells) - cellContents = Cflat{multiCells(i)}; - Cflat(multiCells(i),1:length(cellContents)) = cellContents; -end - -% change empty elements to strings, if specified -if ( strFlag ) - Cflat(cellfun(@isempty,Cflat)) = {''}; +function name = firstSubSystem(subSystems) +% First subSystem of a reaction as a character vector; '' when the reaction +% has none. +if ischar(subSystems) + name = subSystems; +elseif iscell(subSystems) && ~isempty(subSystems) && ischar(subSystems{1}) + name = subSystems{1}; +else + name = ''; end end diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 8ab10779..0ec808fd 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -7,9 +7,10 @@ change_model/getComplexData.m tc0003, tc0004 - File download not covered by test change_model/makeEcModel.m tc0001, tc0002, tc0210, tc0212, tc0233, tc0234 change_model/setProtPoolSize.m tc0005 change_model/getKcatAcrossIsoenzymes.m tc0011, tc0252 +change_model/getSubsetEcModel.m tc0258, tc0259 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/getStandardKcat.m tc0256, tc0257 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 diff --git a/test/unit_tests/parityRound3Tests.m b/test/unit_tests/parityRound3Tests.m index 28191148..b28b74e8 100644 --- a/test/unit_tests/parityRound3Tests.m +++ b/test/unit_tests/parityRound3Tests.m @@ -142,6 +142,83 @@ function testSigmaFitterReturnsStartOfGrowthPlateau_tc0255(testCase) verifyEqual(testCase, plateauSigma(), sigmaGlpk, 'AbsTol', 1e-12) end +function testGetStandardKcatSubSystemNamesIgnoreCase_tc0256(testCase) +% R3 is in 'Alpha', R5 in 'alpha', the entries of R2 in 'Beta'. The names +% 'Alpha' and 'alpha' are one subSystem, so R6 ('ALPHA') gets the median of +% R3 and R5 (5) and not the model-wide median (100). +[ecModel, adapter] = standardKcatModel(struct('R5', 'alpha', 'R2', 'Beta', 'R6', 'ALPHA')); +ecModel = getStandardKcat(ecModel, 'modelAdapter', adapter, 'threshold', 2); +verifyEqual(testCase, ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R6')), 5) +end + +function testGetStandardKcatReactionWithoutSubSystem_tc0257(testCase) +% R7 has an empty subSystem while the other reactions have one: it gets the +% model-wide median kcat (100), R6 the one of its subSystem 'Alpha' (5). +[ecModel, adapter] = standardKcatModel(struct('R2', 'Beta', 'R7', '')); +ecModel = getStandardKcat(ecModel, 'modelAdapter', adapter, 'threshold', 2); +verifyEqual(testCase, ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R7')), 100) +verifyEqual(testCase, ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R6')), 5) +end + +function testGetSubsetEcModelKeepsEverythingForTheFullModel_tc0258(testCase) +% The context model is the conventional model, so no gene is removed and the +% ecModel is returned as it is, for the full and the light ecModel. The +% reaction identifiers of the light ecModel have a prefix; they are matched +% by the reaction they stand for, so they are not all removed. +adapter = getAdapter(); +model = loadConventionalGEM('modelAdapter', adapter); +for light = [false true] + bigEcModel = makeEcModel(model, light, adapter); + smallEcModel = getSubsetEcModel(bigEcModel, model); + verifyEqual(testCase, smallEcModel, bigEcModel) + verifyNotEmpty(testCase, smallEcModel.ec.rxns) +end +end + +function testGetSubsetEcModelLightDropsEcEntriesOfRemovedRxns_tc0259(testCase) +% The context model lacks R2: the light ecModel loses the entry 001_R2 of +% its ec structure and keeps the others. +adapter = getAdapter(); +model = loadConventionalGEM('modelAdapter', adapter); +small = removeReactions(model, {'R2'}, 'removeUnusedMets', true, 'removeUnusedGenes', true); +bigEcModel = makeEcModel(model, true, adapter); +smallEcModel = getSubsetEcModel(bigEcModel, small); +verifyFalse(testCase, ismember('001_R2', smallEcModel.ec.rxns)) +verifyTrue(testCase, all(ismember({'001_R3', '001_R5'}, smallEcModel.ec.rxns))) +verifyEqual(testCase, size(smallEcModel.ec.rxnEnzMat, 1), numel(smallEcModel.ec.rxns)) +end + +function [ecModel, adapter] = standardKcatModel(subSystems) +% ecTestGEM with R6 and R7 added (no genes) and the subSystems given as a +% struct of reaction id to name (an empty name is no subSystem); reactions +% not listed are in 'Alpha'. The kcats +% are R3 2, R5 8, the entries of R2 100. +adapter = getAdapter(); +model = loadConventionalGEM('modelAdapter', adapter); +rxnsToAdd.rxns = {'R6'; 'R7'}; +rxnsToAdd.equations = {'m1[c] => m2[c]'; 'm1[c] => m2[c]'}; +rxnsToAdd.grRules = {''; ''}; +model = addRxns(model, rxnsToAdd, 3); +model.subSystems = cell(numel(model.rxns), 1); +for i = 1:numel(model.rxns) + name = 'Alpha'; + if isfield(subSystems, model.rxns{i}) + name = subSystems.(model.rxns{i}); + end + if isempty(name) + model.subSystems{i} = {}; + else + model.subSystems{i} = {name}; + end +end +ecModel = makeEcModel(model, false, adapter); +ecModel = getECfromGEM(ecModel); +ecModel.ec.kcat(:) = 100; +ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R3')) = 2; +ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R5')) = 8; +ecModel.ec.source(:) = {'manual'}; +end + function sigma = plateauSigma() adapter = getAdapter(); model = loadConventionalGEM('modelAdapter', adapter);