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/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..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 @@ -54,7 +55,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..b28b74e8 100644 --- a/test/unit_tests/parityRound3Tests.m +++ b/test/unit_tests/parityRound3Tests.m @@ -120,6 +120,128 @@ 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 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); +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();