From a923c89625f264691b427f845d1abca0f4544ec4 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:40:54 +0200 Subject: [PATCH 01/14] fix: use the population std for the reporterMetabolites background The closed-form background correction assumes genes drawn with replacement from the scored-gene pool, which needs the population standard deviation. Scores now match raven-toolbox. The help states that results are sorted by Z-score and that P-values are one-sided. --- src/analysis/reporterMetabolites.m | 16 ++++++++++------ tests/function_tests/tAnalysis.m | 18 ++++++++++++++++++ 2 files changed, 28 insertions(+), 6 deletions(-) diff --git a/src/analysis/reporterMetabolites.m b/src/analysis/reporterMetabolites.m index 46ee4eb3..d2a2ae83 100755 --- a/src/analysis/reporterMetabolites.m +++ b/src/analysis/reporterMetabolites.m @@ -36,10 +36,10 @@ % supplied. % - mets : a cell array of metabolite IDs for the metabolites for which % a score could be calculated. -% - metZScores : Z-scores for differential expression around each -% metabolite in "mets". -% - metPValues : P-values for differential expression around each -% metabolite in "mets". +% - metZScores : background-corrected Z-scores for differential +% expression around each metabolite in "mets". The entries of all +% fields are sorted by this score, highest first. +% - metPValues : one-sided (upper-tail) P-values of "metZScores". % - metNGenes : number of neighbouring genes for each metabolite in % "mets". % - meanZ : average Z-scores for the genes around each metabolite in @@ -57,6 +57,10 @@ % For details about the algorithm, see Patil KR, Nielsen J, Uncovering % transcriptional regulation of metabolism by using metabolic network % topology. Proc. Natl Acad. Sci. USA 2005;102:2685-2689. +% +% The background correction uses the closed-form mean and standard +% deviation of the aggregate Z-score of n genes drawn at random from all +% scored genes, so the results are deterministic. genes=convertCharArray(genes); p=parseRAVENargs(varargin, {'printResults',false; 'outputFile',[]; 'geneFoldChanges',[]}); @@ -136,11 +140,11 @@ %Correct for background using the closed-form moments of the sampling %distribution. For n genes drawn with replacement from the scored-gene pool -%(mean mu, std sigma), the aggregate Z = sum(z)/sqrt(n) has: +%(mean mu, population std sigma), the aggregate Z = sum(z)/sqrt(n) has: % E[Z] = sqrt(n)*mu and Std[Z] = sigma (independent of n) %so the corrected score is (metZ - sqrt(n)*mu) / sigma. mu = mean(geneZScores); -sigma = std(geneZScores); +sigma = std(geneZScores,1); metZScores = (metZScores - sqrt(metNGenes) .* mu) ./ sigma; %Calculate the P-values (upper-tail of the standard normal). Uses the base diff --git a/tests/function_tests/tAnalysis.m b/tests/function_tests/tAnalysis.m index d19f8000..a4e825cf 100644 --- a/tests/function_tests/tAnalysis.m +++ b/tests/function_tests/tAnalysis.m @@ -274,6 +274,24 @@ function reporterMetabolitesIsDeterministic(testCase) testCase.verifyEqual(rm1.metZScores, rm2.metZScores); end + function reporterMetabolitesUsesPopulationStd(testCase) + % The closed-form background assumes genes drawn with + % replacement from the scored pool, so it divides by the + % population standard deviation of the gene Z-scores. + m = struct(); + m.mets = {'A';'B'}; + m.metNames = {'A';'B'}; + m.genes = {'g1';'g2';'g3'}; + m.S = sparse([1 0; 0 1]); + m.rxnGeneMat = sparse([1 1 0; 0 0 1]); + pvals = [0.01; 0.2; 0.5]; + z = sqrt(2)*erfcinv(2*pvals); + expectedA = ((z(1)+z(2))/sqrt(2) - sqrt(2)*mean(z)) / std(z,1); + rm = reporterMetabolites(m, m.genes, pvals); + testCase.verifyEqual(rm.metZScores(strcmp(rm.mets,'A')), expectedA, 'AbsTol', 1e-12); + testCase.verifyEqual(rm.metPValues(strcmp(rm.mets,'A')), 0.5*erfc(expectedA/sqrt(2)), 'AbsTol', 1e-12); + end + function FSEOFRuns(testCase) biomassRxn = testCase.model.rxns{find(testCase.model.c == 1, 1)}; [~, exchIdx] = getExchangeRxns(testCase.model); From 52a8b0cd8cd9ef6a48db809e38641e44157d1f3d Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:45:08 +0200 Subject: [PATCH 02/14] fix: pass a sparse constraint matrix from sampleChebyshevCenter Gurobi rejects a dense model.A, so randomSampling with method 'chrr' failed whenever Gurobi was the RAVEN solver. --- src/analysis/sampleChebyshevCenter.m | 2 +- tests/function_tests/tSampling.m | 11 +++++++++++ 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/src/analysis/sampleChebyshevCenter.m b/src/analysis/sampleChebyshevCenter.m index 7204a639..8c8198cd 100644 --- a/src/analysis/sampleChebyshevCenter.m +++ b/src/analysis/sampleChebyshevCenter.m @@ -26,7 +26,7 @@ norms = sqrt(sum(A.^2, 2)); % Variables [x (n, free); r (1, >=0)]; minimise -r (maximises the ball radius). -prob.A = [A, norms]; +prob.A = sparse([A, norms]); % Gurobi and SCIP require a sparse A prob.a = prob.A; prob.b = b; prob.csense = repmat('L', 1, m); diff --git a/tests/function_tests/tSampling.m b/tests/function_tests/tSampling.m index c365addb..797ee052 100644 --- a/tests/function_tests/tSampling.m +++ b/tests/function_tests/tSampling.m @@ -109,6 +109,17 @@ function randomSamplingCHRRDispatch(testCase) testCase.verifyLessThan(max(abs(resid(:))), 1e-6); end + function randomSamplingCHRRWithGurobi(testCase) + % sampleChebyshevCenter must hand Gurobi a sparse constraint + % matrix, which Gurobi requires. + testCase.assumeSolver('gurobi'); + evalc('setRavenSolver(''gurobi'')'); + testCase.addTeardown(@() evalc('setRavenSolver(''glpk'')')); + evalc(['sols = randomSampling(testCase.model, 10, ''method'', ''chrr'', ' ... + '''thinning'', 5, ''nBurnin'', 20, ''seed'', 1);']); + testCase.verifySize(sols, [numel(testCase.model.rxns), 10]); + end + function randomSamplingCHRRExposesInfo(testCase) % The CHRR convergence diagnostics must be reachable through the % documented randomSampling entry point, not only sampleCHRR. From 25f07c939ff59c79fef16b8273db2d63a18730b0 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:46:01 +0200 Subject: [PATCH 03/14] fix: return all-blocked from haveFlux when no reaction survives simplification When every reaction is a dead end, simplifyModel removes them all and the LP has an empty constraint matrix, which GLPK rejects. --- src/analysis/haveFlux.m | 7 +++++++ tests/function_tests/tAnalysis.m | 19 +++++++++++++++++++ 2 files changed, 26 insertions(+) diff --git a/src/analysis/haveFlux.m b/src/analysis/haveFlux.m index 8c005dca..22fec7e6 100755 --- a/src/analysis/haveFlux.m +++ b/src/analysis/haveFlux.m @@ -66,6 +66,13 @@ %First remove all dead-end reactions smallModel=simplifyModel(model,false,false,true,true); +%Nothing is left to solve when every reaction is a dead end, and reactions +%without metabolites are limited by their bounds only +if isempty(smallModel.rxns) || isempty(smallModel.mets) + I=ismember(rxns,smallModel.rxns(smallModel.ub>cutOff | smallModel.lb<-cutOff)); + return +end + %Get the indexes of the reactions in the connected model indexes=getIndexes(smallModel,intersect(smallModel.rxns,rxns),'rxns'); J=false(numel(indexes),1); diff --git a/tests/function_tests/tAnalysis.m b/tests/function_tests/tAnalysis.m index a4e825cf..fe54823a 100644 --- a/tests/function_tests/tAnalysis.m +++ b/tests/function_tests/tAnalysis.m @@ -29,6 +29,25 @@ function haveFluxSeedIsReproducible(testCase) testCase.verifyEqual(I1, I2); end + function haveFluxAllDeadEnds(testCase) + % A model in which every reaction is a dead end has nothing + % left to solve after simplification; all reactions are + % blocked and no solver error is raised. + m = struct(); + m.id = 'deadEnds'; m.rxns = {}; m.S = []; m.rev = []; + m.mets = {'A';'B';'C'}; m.metNames = m.mets; + m.comps = {'c'}; m.compNames = m.comps; m.metComps = [1;1;1]; + m.genes = {}; m.grRules = {}; m.rxnGeneMat = []; + r = struct(); + r.rxns = {'R1';'R2'}; + r.equations = {'A[c] => B[c]';'B[c] => C[c]'}; + evalc('m = addRxns(m, r, 3);'); + m.lb = [0;0]; m.ub = [1000;1000]; + m.c = [0;0]; m.b = zeros(3,1); + I = haveFlux(m); + testCase.verifyEqual(I, [false; false]); + end + function getMinNrFluxesReturnsFlux(testCase) testCase.assumeMILPSolver(); evalc('[x, I, exitFlag] = getMinNrFluxes(testCase.model, testCase.model.rxns);'); From 4878db0de81dc5d35a7f9c6604e2b6dd8940e144 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:46:46 +0200 Subject: [PATCH 04/14] fix: refuse a biomass balance that would make the balancing component negative scaleBiomassFraction rescales the balancing component so the biomass total becomes 1 g/gDW. A required fraction of zero or less inverted the balancing pseudoreaction without notice; it now raises an error, as does a balancing component with no current fraction. The help states that the total becomes 1 g/gDW. --- src/biomass/scaleBiomassFraction.m | 26 +++++++++++++++---- tests/function_tests/tBiomass.m | 41 ++++++++++++++++++++++++++++++ 2 files changed, 62 insertions(+), 5 deletions(-) diff --git a/src/biomass/scaleBiomassFraction.m b/src/biomass/scaleBiomassFraction.m index ba360243..6c66f864 100644 --- a/src/biomass/scaleBiomassFraction.m +++ b/src/biomass/scaleBiomassFraction.m @@ -2,7 +2,8 @@ % scaleBiomassFraction Rescale a biomass component to a target value. % % Rescale a biomass component to a target g/gDW value, optionally -% balancing a second component so the total biomass mass stays at 1 g/gDW. +% rescaling a second component so that the biomass total becomes 1 g/gDW, +% whatever the total was before. % The MATLAB counterpart of yeast-GEM's legacy scaleBioMass. % % Parameters @@ -19,8 +20,10 @@ % Name-Value Arguments % -------------------- % balanceOut : char -% Second component name to adjust so the biomass total remains 1 -% g/gDW. Empty / omit to skip balancing. +% Second component name to adjust so the biomass total becomes 1 +% g/gDW. It must be large enough to absorb the change: an error is +% raised if its fraction would become zero or negative, or if it is +% 0 to begin with. Empty / omit to skip balancing. % % Returns % ------- @@ -55,11 +58,24 @@ if ~isempty(balanceOut) fractions = getBiomassFractions(model, biomassConfig); total = fractions.total; + if ~isfield(fractions, balanceOut) + error('scaleBiomassFraction:unknownComponent', ... + 'biomassConfig has no component named %s', balanceOut); + end balanceCurrent = fractions.(balanceOut); if balanceCurrent == 0 - return; + error('scaleBiomassFraction:zeroCurrent', ... + ['Cannot balance with %s: current fraction is 0 ' ... + '(pseudoreaction missing or empty).'], balanceOut); + end + balanceNew = balanceCurrent + (1 - total); + if balanceNew <= 0 + error('scaleBiomassFraction:negativeBalance', ... + ['Cannot balance with %s: its fraction would become %g g/gDW. ' ... + 'It is %g g/gDW, but the other components total %g g/gDW.'], ... + balanceOut, balanceNew, balanceCurrent, total - balanceCurrent); end - balanceFactor = (balanceCurrent + (1 - total)) / balanceCurrent; + balanceFactor = balanceNew / balanceCurrent; model = scaleBiomassPseudoreaction(model, biomassConfig, balanceOut, balanceFactor); end end diff --git a/tests/function_tests/tBiomass.m b/tests/function_tests/tBiomass.m index dd9c7276..7b7c0faa 100644 --- a/tests/function_tests/tBiomass.m +++ b/tests/function_tests/tBiomass.m @@ -26,6 +26,22 @@ function scaleBiomassFractionErrorsOnMissingComponent(testCase) 'protein', 0.5), ?MException); end + function scaleBiomassFractionBalancesToOne(testCase) + [m, cfg] = tBiomass.twoComponentToyModel(); + out = scaleBiomassFraction(m, cfg, 'protein', 0.7, 'balanceOut', 'DNA'); + f = getBiomassFractions(out, cfg); + testCase.verifyEqual(f.protein, 0.7, 'AbsTol', 1e-9); + testCase.verifyEqual(f.DNA, 0.3, 'AbsTol', 1e-9); + end + + function scaleBiomassFractionRefusesNegativeBalance(testCase) + % Balancing protein 1.2 against DNA would need DNA at -0.2, + % which would turn the DNA pseudoreaction around. + [m, cfg] = tBiomass.twoComponentToyModel(); + testCase.verifyError(@() scaleBiomassFraction(m, cfg, 'protein', 1.2, ... + 'balanceOut', 'DNA'), 'scaleBiomassFraction:negativeBalance'); + end + function scaleBiomassPseudoreactionErrorsOnMissingComponent(testCase) cfg = testCase.biomassConfig(); testCase.verifyError(@() scaleBiomassPseudoreaction(testCase.model, cfg, ... @@ -93,6 +109,31 @@ function fitParametersRuns(testCase) end methods (Static, Access = private) + function [m, cfg] = twoComponentToyModel() + % Protein and DNA pseudoreactions in grams, 0.5 g/gDW each. + m = struct(); + m.id = 'toy'; + m.rxns = {'R1';'R2'}; + m.rxnNames = {'protein pseudoreaction';'DNA pseudoreaction'}; + m.mets = {'x_c';'y_c';'h_c';'prot_c';'dna_c'}; + m.metNames = {'x';'y';'H+';'protein';'DNA'}; + m.metComps = ones(5,1); + m.metCharges= zeros(5,1); + m.comps = {'c'}; + m.compNames = {'cytosol'}; + m.S = sparse([-0.5 0; 0 -0.5; 0 0; 1 0; 0 1]); + m.lb = [0;0]; m.ub = [1000;1000]; m.rev = [0;0]; m.c = [0;0]; m.b = zeros(5,1); + m.grRules = {'';''}; m.genes = {}; m.rxnGeneMat = sparse(2,0); + + cfg = struct(); + cfg.biomass_rxn = 'R1'; + cfg.proton_met = 'h_c'; + cfg.components = {struct('name','protein', ... + 'pseudoreaction_name','protein pseudoreaction','mass_strategy','grams'), ... + struct('name','DNA', ... + 'pseudoreaction_name','DNA pseudoreaction','mass_strategy','grams')}; + end + function [m, cfg] = protonToyModel() % One pseudoreaction: x_c -> protein, with h_c available to carry % the charge balance. From da9c35a126d8f29078d6bf494113696d6107f850 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:48:58 +0200 Subject: [PATCH 05/14] feat: solve MILPs with the bundled GLPK optimizeProb no longer refuses a MILP when 'glpk' is the RAVEN solver. The bundled GLPK MEX supports integer variables and gives the same answers as Gurobi on the gap-filling and minimal-medium problems of the guide models. A RAVEN:glpkMILP warning recommends Gurobi for larger problems. The GLPK branch uses the default time limit for MILPs and reports a mipgap, as the Gurobi and SCIP branches do. MILPs through the COBRA Toolbox's GLPK interface are still refused. The help of setRavenSolver, gapFillMILP, fillGaps and getMinimalMedium states the solver behaviour; getMinimalMedium also states that the returned medium need not be unique. --- src/analysis/getMinimalMedium.m | 10 ++++++++ src/gapfilling/fillGaps.m | 7 ++++++ src/gapfilling/gapFillMILP.m | 6 +++++ src/solver/optimizeProb.m | 36 +++++++++++++++++++++------- src/solver/setRavenSolver.m | 9 +++++-- tests/function_tests/RavenTestCase.m | 4 ++-- tests/function_tests/tSolver.m | 36 ++++++++++++++++++++++++++++ 7 files changed, 95 insertions(+), 13 deletions(-) diff --git a/src/analysis/getMinimalMedium.m b/src/analysis/getMinimalMedium.m index 754c8c30..186f044b 100644 --- a/src/analysis/getMinimalMedium.m +++ b/src/analysis/getMinimalMedium.m @@ -50,6 +50,16 @@ % % A y_i = 0 forces reaction i to carry non-negative flux (no net uptake); % y_i = 1 allows it to carry flux as low as lb_i. +% +% The minimal number of uptakes is guaranteed, but the medium itself need +% not be unique: when several media of the same size support the growth +% rate, the solver returns one of them, and another solver or solver +% version can return a different one. +% +% The MILP is solved with the solver set by setRavenSolver. RAVEN's bundled +% GLPK solves it for small models (with a RAVEN:glpkMILP warning); for +% larger models Gurobi is much faster and more reliable +% (setRavenSolver('gurobi')). p = parseRAVENargs(varargin, {'minGrowth',[]; 'verbose',true}); minGrowth = p.minGrowth; diff --git a/src/gapfilling/fillGaps.m b/src/gapfilling/fillGaps.m index 56480cf1..ee9f13bb 100755 --- a/src/gapfilling/fillGaps.m +++ b/src/gapfilling/fillGaps.m @@ -93,6 +93,13 @@ % fillGaps(model,models,allowNetProduction,useModelConstraints,... % supressWarnings,rxnScores,params); % +% Notes +% ----- +% The 'reference' and 'gapfillMILP' algorithms solve a MILP with the +% solver set by setRavenSolver. RAVEN's bundled GLPK solves it for small +% models (with a RAVEN:glpkMILP warning); for larger models Gurobi is much +% faster and more reliable (setRavenSolver('gurobi')). +% % See Also % -------- % ftINITFillGaps : the ftINIT-side counterpart, used when the reference diff --git a/src/gapfilling/gapFillMILP.m b/src/gapfilling/gapFillMILP.m index 0127fe1a..e66d228c 100644 --- a/src/gapfilling/gapFillMILP.m +++ b/src/gapfilling/gapFillMILP.m @@ -52,6 +52,12 @@ % exitFlag : double % 1 = optimal, -1 = infeasible, -2 = timeout or other failure. % +% Notes +% ----- +% The MILP is solved with the solver set by setRavenSolver. RAVEN's bundled GLPK solves it for small +% models (with a RAVEN:glpkMILP warning); for larger models Gurobi is much +% faster and more reliable (setRavenSolver('gurobi')). +% % See also % -------- % gapFillFastLP, gapFillTopological, fillGaps diff --git a/src/solver/optimizeProb.m b/src/solver/optimizeProb.m index de278ad1..31991dad 100755 --- a/src/solver/optimizeProb.m +++ b/src/solver/optimizeProb.m @@ -27,6 +27,15 @@ % time limit. Whatever solution is in res is then either absent or not % proven optimal. % +% Notes +% ----- +% With 'glpk' selected, MILPs are solved by RAVEN's bundled GLPK, and a +% RAVEN:glpkMILP warning recommends Gurobi, which is much faster and more +% reliable on larger MILPs (such as ftINIT or getMinNrFluxes on +% genome-scale models). +% The GLPK interface of the COBRA Toolbox is not used for MILPs: with +% 'cobra' selected, set a different MILP solver with changeCobraSolver. +% % See also % -------- % splitProbForConditioning @@ -53,13 +62,16 @@ solver=RAVENSOLVER; if ~all(lower(prob.vartype) == 'c') milp=true; - errorText = 'glpk is not suitable for solving MILPs, '; switch solver case 'glpk' - error([errorText 'select a different solver with setRavenSolver().']) + warning('RAVEN:glpkMILP', ['Solving a MILP with GLPK. This works for small ' ... + 'problems, but for larger problems Gurobi is much faster and more reliable; ' ... + 'select it with setRavenSolver(''gurobi''). Silence this warning with ' ... + 'warning(''off'',''RAVEN:glpkMILP'').']) case 'cobra' if strcmp(CBT_MILP_SOLVER,'glpk') - error([errorText 'select a different solver with changeCobraSolver() or setRavenSolver().']) + error(['The GLPK interface of the COBRA Toolbox is not suitable for solving MILPs, ' ... + 'select a different solver with changeCobraSolver() or setRavenSolver().']) end end else @@ -211,16 +223,13 @@ solverparams.tolint = defaultparams.intTol; solverparams.tolobj = defaultparams.objTol; solverparams.msglev = 0; % Level of verbosity + if milp + solverparams.tmlim = defaultparams.timeLimit; % seconds + end solverparams = structUpdate(solverparams,params); prob.csense = renameparams(prob.csense, {'L','G','E'}, {'U','L','S'}); - if milp - solverparams.tmlim = solverparams.tmlim*10; - solverparams.msglev = 1; % Level of verbosity - disp('Issues have been observed when using GLPK for MILP solving. Be advised to carefully observe the results, or us another solver.') - end - % Ensure that RAVEN glpk binary is used, return to original % directory afterwards [ravenDir,currDir]=findRAVENroot(); @@ -241,6 +250,15 @@ res.obj = fmin; res.dual = -extra.lambda*prob.osense; res.rcost = -extra.redcosts*prob.osense; + if milp + % glpk reports no gap; a proven optimum has none, anything else + % is not accepted by callers that test the gap + if res.stat == 1 + res.mipgap = 0; + else + res.mipgap = inf; + end + end %% Use scip case {'soplex','scip'} % Old 'soplex' option also allowed % scip takes a row range [rl,ru] rather than a right-hand side plus diff --git a/src/solver/setRavenSolver.m b/src/solver/setRavenSolver.m index 0adf3663..dc534e91 100755 --- a/src/solver/setRavenSolver.m +++ b/src/solver/setRavenSolver.m @@ -9,11 +9,16 @@ function setRavenSolver(solver) % solver : char % string with solver name. The following options are available: % -% - glpk : uses RAVEN-provided binaries. +% - glpk : uses RAVEN-provided binaries. Solves LPs and MILPs; for +% MILPs a RAVEN:glpkMILP warning recommends Gurobi, which is much +% faster and more reliable on larger problems. % - gurobi : requires a working Gurobi installation. +% - scip : requires a SCIP MEX interface. RAVEN provides a +% precompiled binary for Windows only. % - cobra : requires a working COBRA Toolbox installation. RAVEN % parses the problem to COBRA and uses the solver that is set by -% changeCobraSolver. +% changeCobraSolver. MILPs are refused when COBRA's MILP solver is +% GLPK. % % Examples % -------- diff --git a/tests/function_tests/RavenTestCase.m b/tests/function_tests/RavenTestCase.m index f6af4ee7..245e896d 100644 --- a/tests/function_tests/RavenTestCase.m +++ b/tests/function_tests/RavenTestCase.m @@ -71,8 +71,8 @@ function assumeDependency(testCase, isAvailable, label) end function assumeMILPSolver(testCase) - % Switch to a MILP-capable solver for this test, or skip if none. - % GLPK (the default) cannot solve MILPs in RAVEN. + % Switch to Gurobi or SCIP for this test, or skip if neither is + % available. if ~isempty(which('gurobi')) evalc('setRavenSolver(''gurobi'')'); else diff --git a/tests/function_tests/tSolver.m b/tests/function_tests/tSolver.m index b86eb529..dbe99a0b 100644 --- a/tests/function_tests/tSolver.m +++ b/tests/function_tests/tSolver.m @@ -37,6 +37,42 @@ function optimizeProbReportsTimeLimit(testCase) testCase.verifyFalse(res.hitTimeLimit); end + function optimizeProbSolvesMILPWithGlpk(testCase) + % A binary knapsack whose LP relaxation (9.33) differs from the + % integer optimum (8, items 1 and 3). GLPK solves it, with a + % warning that recommends Gurobi. + evalc('setRavenSolver(''glpk'')'); + prob = struct(); + prob.c = [-5; -4; -3]; + prob.A = sparse([2 3 1]); + prob.a = prob.A; + prob.b = 4; + prob.csense = 'L'; + prob.lb = zeros(3,1); + prob.ub = ones(3,1); + prob.osense = 1; + prob.vartype = 'BBB'; + res = testCase.verifyWarning(@() optimizeProb(prob,[],false), 'RAVEN:glpkMILP'); + testCase.verifyEqual(res.stat, 1); + testCase.verifyEqual(res.obj, -8, 'AbsTol', 1e-9); + testCase.verifyEqual(res.full, [1; 0; 1], 'AbsTol', 1e-9); + testCase.verifyEqual(res.mipgap, 0); + end + + function getMinimalMediumRunsOnGlpk(testCase) + evalc('setRavenSolver(''glpk'')'); + warnState = warning('off','RAVEN:glpkMILP'); + testCase.addTeardown(@() warning(warnState)); + [medium, idx] = getMinimalMedium(testCase.model, 'verbose', false); + testCase.verifyNotEmpty(medium); + m = testCase.model; + inMedium = false(numel(m.rxns),1); inMedium(idx) = true; + [~, exchIdx] = getExchangeRxns(m); + m.lb(intersect(exchIdx, find(~inMedium & m.lb < 0))) = 0; + sol = solveLP(m); + testCase.verifyGreaterThan(sol.f, 0); + end + function checkSolutionOnFeasibleProblem(testCase) res = optimizeProb(testCase.modelLP()); [isFeasible, isOptimal] = checkSolution(res); From b347457dc09189a2cf1738e310c3a3989f21e1a4 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:49:21 +0200 Subject: [PATCH 06/14] fix: raise a clear error for an ambiguous name[comp] in getIndexes A metaboliteName[comp] that matches more than one metabolite ended in an unexplained assignment error. The error now names the matching metabolites. --- src/queries/getIndexes.m | 10 ++++++++-- tests/function_tests/tQueries.m | 23 +++++++++++++++++++++++ 2 files changed, 31 insertions(+), 2 deletions(-) diff --git a/src/queries/getIndexes.m b/src/queries/getIndexes.m index 87770bbc..1c64efa5 100755 --- a/src/queries/getIndexes.m +++ b/src/queries/getIndexes.m @@ -12,6 +12,8 @@ % 'rxns', 'mets', or 'genes' depending on what to retrieve. 'metnames' % queries metabolite names, while 'metcomps' allows providing specific % metabolites and their compartments in the format metaboliteName[comp]. +% An error is raised if such a name matches no metabolite, or more than +% one. % If a model.ec structure exists (GECKO 3), then also 'ecenzymes', % 'ecrxns' and 'ecgenes' are allowed. % @@ -65,10 +67,14 @@ for i=1:numel(objects) index=find(strcmp(mets(i),model.metNames)); index=index(strcmp(comps(i),model.comps(model.metComps(index)))); - if ~isempty(index) + if isscalar(index) indexes(i)=index; + elseif isempty(index) + error(['Could not find object ''' objects{i} ''' in the model']); else - error(['Could not find object ' objects{i} ' in the model']); + error(['''' objects{i} ''' is ambiguous: it matches the ' ... + 'metabolites ' strjoin(model.mets(index),', ') ... + '. Use metabolite ids (type ''mets'') instead.']); end end indexes=indexes(:); diff --git a/tests/function_tests/tQueries.m b/tests/function_tests/tQueries.m index 1fc2c5a6..7e918050 100644 --- a/tests/function_tests/tQueries.m +++ b/tests/function_tests/tQueries.m @@ -293,6 +293,29 @@ function getIndexesByName(testCase) testCase.verifyEqual(idx, [1;2;3]); end + function getIndexesMetcompsFindsName(testCase) + m = testCase.model; + name = [m.metNames{1} '[' m.comps{m.metComps(1)} ']']; + testCase.verifyEqual(getIndexes(m, name, 'metcomps'), 1); + end + + function getIndexesMetcompsAmbiguousNameErrors(testCase) + % Two metabolites with the same name in the same compartment + % cannot be told apart by name[comp]; the error names both. + m = testCase.model; + m.metNames{2} = m.metNames{1}; + m.metComps(2) = m.metComps(1); + name = [m.metNames{1} '[' m.comps{m.metComps(1)} ']']; + try + getIndexes(m, name, 'metcomps'); + testCase.verifyFail('No error for an ambiguous name[comp]'); + catch e + testCase.verifySubstring(e.message, 'ambiguous'); + testCase.verifySubstring(e.message, m.mets{1}); + testCase.verifySubstring(e.message, m.mets{2}); + end + end + function getIndexesLogical(testCase) lg = getIndexes(testCase.model, testCase.model.rxns(1), 'rxns', true); testCase.verifyClass(lg, 'logical'); From fef2abaf87e9f140f36a157e27f9c3981d9de57e Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:49:59 +0200 Subject: [PATCH 07/14] docs: correct the help of six functions - randomSampling: the random-objective description applies to 'randomObjective' only, not to the default ACHR. - checkTasks: the model's own exchange reactions are not closed. - findGeneDeletions: details 2 means the solve failed, not that the deletion is lethal. - checkSolution: it reads an optimizeProb result, not a Mosek one. - assignSBOterms: where the terms are stored and how to read them. - haveFlux: the test uses the current bounds. --- src/analysis/findGeneDeletions.m | 6 ++++-- src/analysis/haveFlux.m | 9 +++++---- src/analysis/randomSampling.m | 9 ++++++--- src/annotation/assignSBOterms.m | 5 ++++- src/solver/checkSolution.m | 9 +++------ src/tasks/checkTasks.m | 6 ++++++ 6 files changed, 28 insertions(+), 16 deletions(-) diff --git a/src/analysis/findGeneDeletions.m b/src/analysis/findGeneDeletions.m index dd0e6602..14f4c0fa 100755 --- a/src/analysis/findGeneDeletions.m +++ b/src/analysis/findGeneDeletions.m @@ -36,8 +36,10 @@ % with details about each gene in originalGenes and why or why not it % was deleted: % -% - 1 : was deleted -% - 2 : proved lethal in sgd (single gene deletion) +% - 1 : was deleted, and the problem was solved (a lethal deletion +% with zero growth also gets 1) +% - 2 : the single-deletion problem had no optimal solution +% (infeasible or failed); fluxes and grRatioMuts are 0 % - 3 : redundant, no longer used % - 4 : involved in dead-end reaction % grRatioMuts : double diff --git a/src/analysis/haveFlux.m b/src/analysis/haveFlux.m index 22fec7e6..aa0ae4f0 100755 --- a/src/analysis/haveFlux.m +++ b/src/analysis/haveFlux.m @@ -1,9 +1,10 @@ function I=haveFlux(model,varargin) % haveFlux Check which reactions can carry a flux. % -% Checks which reactions can carry a (positive or negative) flux. Is used as -% a faster version of getAllowedBounds if it is only interesting whether the -% reactions can carry a flux or not. +% Checks which reactions can carry a (positive or negative) flux under the +% model's current bounds; other media or bounds are not explored. Is used +% as a faster version of getAllowedBounds if it is only interesting whether +% the reactions can carry a flux or not. % % Parameters % ---------- @@ -28,7 +29,7 @@ % ------- % I : logical % logical array with true if the corresponding reaction can carry a -% flux. +% flux under the tested bounds. % % Examples % -------- diff --git a/src/analysis/randomSampling.m b/src/analysis/randomSampling.m index 95f674e6..8d3347fe 100755 --- a/src/analysis/randomSampling.m +++ b/src/analysis/randomSampling.m @@ -81,9 +81,12 @@ % % Notes % ----- -% The solutions are generated by maximizing (with random weights) for a -% random set of nObjectives reactions (default 2). For reversible reactions -% it randomly chooses between maximizing and minimizing. +% With 'randomObjective', the solutions are generated by maximizing (with +% random weights) for a random set of nObjectives reactions (default 2). +% For reversible reactions it randomly chooses between maximizing and +% minimizing. The default 'achr' and the 'chrr' methods use no random +% objectives; they take hit-and-run steps through the interior of the +% solution space. % % If the model is a GECKO v3+ ecModel, then usage_prot reactions are not % selected for sampling, instead focusing on sampling from the metabolic diff --git a/src/annotation/assignSBOterms.m b/src/annotation/assignSBOterms.m index 6aa0a983..d2a42190 100644 --- a/src/annotation/assignSBOterms.m +++ b/src/annotation/assignSBOterms.m @@ -6,7 +6,10 @@ % organism-agnostic, parameterised entirely by opts. % % SBO is written via editMiriam(..., 'fill') so pre-existing SBO -% annotations are preserved. +% annotations are preserved. The terms are stored in model.metMiriams and +% model.rxnMiriams under the namespace 'sbo' (there is no separate SBO +% field); read them with extractMiriam(model.rxnMiriams, 'sbo'). +% exportModel writes them as the SBML sboTerm attribute. % % Parameters % ---------- diff --git a/src/solver/checkSolution.m b/src/solver/checkSolution.m index 3fb59b49..e6b3478b 100755 --- a/src/solver/checkSolution.m +++ b/src/solver/checkSolution.m @@ -1,10 +1,10 @@ function [isFeasible, isOptimal]=checkSolution(res) -% checkSolution Check if a Mosek solution is feasible and optimal. +% checkSolution Check if a solver result is feasible and optimal. % % Parameters % ---------- % res : struct -% the output structure from mosekopt. +% the output structure from optimizeProb. % % Returns % ------- @@ -16,10 +16,7 @@ % Examples % -------- % [isFeasible, isOptimal] = checkSolution(res); -% -% Notes -% ----- -% This function also throws an error if the license has expired. + isFeasible=false; isOptimal=false; diff --git a/src/tasks/checkTasks.m b/src/tasks/checkTasks.m index 3f8a3804..e9e309bc 100755 --- a/src/tasks/checkTasks.m +++ b/src/tasks/checkTasks.m @@ -8,6 +8,12 @@ % one new equation, and one change of reaction bounds. If more bounds are % needed to define the task, then several rows can be used for each task. % +% The task inputs and outputs are applied to the boundary metabolites of +% the model. The model's own exchange reactions are not closed: they keep +% their bounds, so in a model without boundary metabolites they stay open +% during every task. To let the task define the whole boundary, run +% checkTasks on closeModel(model). +% % Parameters % ---------- % model : struct From 2a16c3b5d66de314a34e714e63266395a341a73f Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:50:18 +0200 Subject: [PATCH 08/14] fix: read folded YAML scalars whose continuation starts with a dash readYAMLmodel treated any indented line starting with '-' as a list item, so a wrapped name continuing with '-epoxidizing)' was reported as an unknown entry. Only a dash followed by a space, or a lone dash, starts a list item. This made the kegg118 reference model unreadable and with it getKEGGModelForOrganism. --- src/io/readYAMLmodel.m | 8 +++++--- tests/function_tests/tIO.m | 19 +++++++++++++++++++ 2 files changed, 24 insertions(+), 3 deletions(-) diff --git a/src/io/readYAMLmodel.m b/src/io/readYAMLmodel.m index 9d1dcb43..90a3333c 100755 --- a/src/io/readYAMLmodel.m +++ b/src/io/readYAMLmodel.m @@ -58,9 +58,11 @@ % A continuation line is recognised by NOT starting with '- ': every real % entry in this format, at every nesting depth (including a nested list % under a key, e.g. an annotation sub-entry indented 6+ spaces), is itself -% introduced with a literal leading dash, while a folded scalar's -% continuation is raw text with no list-item marker of its own. -newLine=regexp(line_raw,'^ {6,}([^-\s].*)','tokenExtents'); +% introduced with a list-item marker (a dash followed by a space, or a dash +% ending the line), while a folded scalar's continuation is raw text with no +% list-item marker of its own. Such text can itself start with a dash, as in +% a KEGG reaction name wrapped before "-epoxidizing)". +newLine=regexp(line_raw,'^ {6,}((?!-(\s|$))\S.*)','tokenExtents'); brokenLine=find(~cellfun('isempty',newLine)); for i=flip(1:numel(brokenLine)) extraLine = char(line_raw(brokenLine(i))); diff --git a/tests/function_tests/tIO.m b/tests/function_tests/tIO.m index f26ace24..86a06ae6 100644 --- a/tests/function_tests/tIO.m +++ b/tests/function_tests/tIO.m @@ -207,6 +207,25 @@ function writeReadYAMLRoundTrip(testCase) testCase.verifyEqual(numel(m2.rxns), numel(testCase.model.rxns)); end + function readYAMLmodelJoinsContinuationStartingWithDash(testCase) + % A writer that wraps long scalars can break a name right + % before a dash, so the continuation line starts with "-" + % without being a list item (as in the KEGG reference model). + m = testCase.model; + m.rxnNames{1} = 'PLACEHOLDER'; + f = [tempname '.yml']; + testCase.addTeardown(@() delete(f)); + evalc('writeYAMLmodel(m, f);'); + txt = fileread(f); + txt = strrep(txt, 'name: PLACEHOLDER', sprintf(['name: naphthalene,NADPH:oxygen ' ... + 'oxidoreductase (RH-hydroxylating or \n -epoxidizing)'])); + fid = fopen(f, 'w'); fwrite(fid, txt); fclose(fid); + evalc('m2 = readYAMLmodel(f);'); + testCase.verifyEqual(m2.rxnNames{1}, ... + 'naphthalene,NADPH:oxygen oxidoreductase (RH-hydroxylating or -epoxidizing)'); + testCase.verifyEqual(numel(m2.rxns), numel(m.rxns)); + end + function writeReadYAMLNestsEcRxnEnzymes(testCase) % Each ec-rxns entry lists its enzymes as a mapping under the % "enzymes" key, so the enzyme lines are indented deeper than From 33d5272b62860eafea43fcd1acd26d9892a2ab88 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:50:18 +0200 Subject: [PATCH 09/14] fix: read curation tables by header and curate existing entities fully curateModelFromTables read its tables by column position, so a metabolite table without an inchi column wrote its notes over the InChIs. Columns are now found by header name, in any order, and optional columns can be left out. Every supplied field is written to existing metabolites, genes and reactions as well as to new ones, including metNotes and rxnReferences, which were never written to existing entities; a field the model lacks is created. MIRIAM columns of the reaction table are taken from the columns after the core fields rather than from a fixed position, and grRules are matched to the right existing reaction when existing and new reactions are mixed. --- src/curation/curateModelFromTables.m | 383 ++++++++++++++++----------- tests/function_tests/tCuration.m | 86 ++++++ 2 files changed, 308 insertions(+), 161 deletions(-) diff --git a/src/curation/curateModelFromTables.m b/src/curation/curateModelFromTables.m index 2cda31fe..47c6d976 100644 --- a/src/curation/curateModelFromTables.m +++ b/src/curation/curateModelFromTables.m @@ -12,28 +12,39 @@ % This function can therefore be used to add new entities in the model, or % curate those already existing in the model. % +% Columns are identified by their header (case-insensitive), so their order +% does not matter and optional columns can be left out. Every field that a +% table supplies is written to new and existing entities alike; a field the +% model does not have yet is created, with its default value for the other +% entities, unless all supplied values for it are empty. An empty grRules +% cell leaves the rule of an existing reaction unchanged. +% % Parameters % ---------- % model : struct % RAVEN model structure to be curated. % metsInfo : char % Path to a *.tsv file with metabolite information, or 'none' to skip -% metabolite curation. Columns: metNames, comps, formula, charge, -% inchi, metNotes, then any number of MIRIAM-namespace columns. +% metabolite curation. Columns: metNames, comps (both required), +% formula, charge, inchi, metNotes, and any number of MIRIAM-namespace +% columns. Singular forms and the RAVEN field names (metFormulas, +% metCharges, inchis) are accepted as headers too. % % Name-Value Arguments % -------------------- % genesInfo : char % Path to a *.tsv file with gene information, or 'none'. Columns: -% genes, geneShortNames, then MIRIAM. +% genes (required), geneShortNames (or geneShortName), and MIRIAM. % rxnsCoeffs : char % Path to a *.tsv file with reaction stoichiometric coefficients, or -% 'none'. Columns: rxnIdx, rxnNames, metNames, comps, coefficient. One -% row per (reaction, metabolite) pair. +% 'none'. Columns (all required): rxnIdx (or index), rxnNames, +% metNames, comps, coefficient. One row per (reaction, metabolite) +% pair. % rxnsInfo : char % Path to a *.tsv file with reaction information, or 'none'. Columns: -% rxnIdx, rxnNames, grRules, lb, ub, rev, subSystems, eccodes, -% rxnNotes, rxnReferences, rxnConfidenceScores, then MIRIAM. +% rxnIdx (or index), rxnNames (both required), grRules, lb, ub, rev, +% subSystems, eccodes, rxnNotes, rxnReferences, rxnConfidenceScores, +% and MIRIAM. % metPrefix : char % Prefix used to mint fresh metabolite ids (e.g. 's_' for yeast-GEM, % 'M_' for the BiGG default) (default 'M_'). @@ -52,10 +63,9 @@ % % Notes % ----- -% The 'everything after the core columns is MIRIAM' convention applies to -% all three info tables: any column whose header is not one of the listed -% core fields is treated as a MIRIAM annotation namespace and stored on -% the matching entity. +% In all three info tables, any column whose header is not one of the +% listed core fields is treated as a MIRIAM annotation namespace and +% stored on the matching entity. p=parseRAVENargs(varargin, {'genesInfo','none'; 'rxnsCoeffs','none'; 'rxnsInfo','none'; 'metPrefix',[]; 'rxnPrefix',[]}); genesInfo=p.genesInfo; @@ -76,25 +86,32 @@ %% Metabolites if ~strcmp(metsInfo,'none') - fid = fopen(metsInfo); - raw = textscan(fid,['%q' repmat(' %q',1,17)],'Delimiter','\t'); - fclose(fid); - metsToAdd.metNames = raw{1}(2:end); - metsToAdd.compartments = raw{2}(2:end); - metsToAdd.metFormulas = raw{3}(2:end); - metsToAdd.metCharges = cellfun(@str2num,raw{4}(2:end),'UniformOutput',false); - emptyEntry = cellfun(@isempty,metsToAdd.metCharges); - if all(emptyEntry) - metsToAdd = rmfield(metsToAdd,'metCharges'); - elseif any(emptyEntry) % If some charges are given, assume 0 for those without charges specified - metsToAdd.metCharges(emptyEntry) = {0}; - metsToAdd.metCharges = cell2mat(metsToAdd.metCharges); - else - metsToAdd.metCharges = cell2mat(metsToAdd.metCharges); + [header,data] = readCurationTable(metsInfo); + [col,miriamCols] = matchColumns(header,metsInfo,{'metNames','comps'},{ ... + 'metNames', {'metNames','metName'}; ... + 'comps', {'comps','comp'}; ... + 'metFormulas', {'formula','formulas','metFormula','metFormulas'}; ... + 'metCharges', {'charge','charges','metCharge','metCharges'}; ... + 'inchis', {'inchi','inchis'}; ... + 'metNotes', {'metNotes','metNote'}}); + metsToAdd.metNames = data(:,col.metNames); + metsToAdd.compartments = data(:,col.comps); + textFields = {'metFormulas','inchis','metNotes'}; + for f = textFields(isfield(col,textFields)) + metsToAdd.(f{1}) = data(:,col.(f{1})); + end + if isfield(col,'metCharges') + metsToAdd.metCharges = cellfun(@str2num,data(:,col.metCharges),'UniformOutput',false); + emptyEntry = cellfun(@isempty,metsToAdd.metCharges); + if all(emptyEntry) + metsToAdd = rmfield(metsToAdd,'metCharges'); + else % If some charges are given, assume 0 for those without charges specified + metsToAdd.metCharges(emptyEntry) = {0}; + metsToAdd.metCharges = cell2mat(metsToAdd.metCharges); + end end - metsToAdd.inchis = raw{5}(2:end); - metsToAdd.metNotes = raw{6}(2:end); - + optFields = intersect({'metFormulas','metCharges','inchis','metNotes'},fieldnames(metsToAdd)','stable'); + % Check if metabolite already exists (check by metName[comp]) existingMets = []; existingMetsIdx = []; @@ -109,133 +126,135 @@ newMets = [newMets, i]; end end - + % Overwrite annotation of existing entries if any(existingMets) - if isfield(newModel,'metFormulas') - newModel.metFormulas(existingMetsIdx) = metsToAdd.metFormulas(existingMets); - end - if isfield(newModel,'metCharges') && isfield(metsToAdd,'metCharges') - newModel.metCharges(existingMetsIdx) = metsToAdd.metCharges(existingMets); + for f = optFields + newModel = writeExistingField(newModel,f{1},existingMetsIdx,metsToAdd.(f{1})(existingMets)); end - if isfield(newModel,'inchis') - newModel.inchis(existingMetsIdx) = metsToAdd.inchis(existingMets); - end - newModel = extractAndAddMiriam(newModel,raw(7:end),existingMets,existingMetsIdx,'met'); + newModel = extractAndAddMiriam(newModel,header(miriamCols),data(:,miriamCols),existingMets,existingMetsIdx,'met'); warning(['The following metabolites are already present in the model, '... 'their annotation will be overwritten to match the metsInfo file. '... 'If you do not particularly want to curate their annotations, it '... 'would be better to removes these metabolites from metsInfo:\n\t\t%s'],... strjoin(strcat(metsToAdd.metNames(existingMets),'[',metsToAdd.compartments(existingMets),']'),'\n\t\t')); end - + % Continue with new metabolites if numel(newMets)>0 - metsToAdd.metNames(existingMets) = []; - metsToAdd.compartments(existingMets) = []; - metsToAdd.metFormulas(existingMets) = []; - if all(cellfun(@isempty,metsToAdd.metFormulas)) - metsToAdd = rmfield(metsToAdd,'metFormulas'); - end - if isfield(metsToAdd,'metCharges') - metsToAdd.metCharges(existingMets) = []; - end - metsToAdd.inchis(existingMets) = []; - if all(cellfun(@isempty,metsToAdd.inchis)) - metsToAdd = rmfield(metsToAdd,'inchis'); - end - metsToAdd.metNotes(existingMets) = []; - if all(cellfun(@isempty,metsToAdd.metNotes)) - metsToAdd = rmfield(metsToAdd,'metNotes'); + metsToAdd.metNames = metsToAdd.metNames(newMets); + metsToAdd.compartments = metsToAdd.compartments(newMets); + for f = optFields + metsToAdd.(f{1}) = metsToAdd.(f{1})(newMets); + if iscell(metsToAdd.(f{1})) && all(cellfun(@isempty,metsToAdd.(f{1}))) + metsToAdd = rmfield(metsToAdd,f{1}); + end end % Add metabolites newModel = addMets(newModel,metsToAdd,true,metPrefix); addedIdx = numel(newModel.mets)-numel(newMets)+1:numel(newModel.mets); - newModel = extractAndAddMiriam(newModel,raw(7:end),newMets,addedIdx,'met'); + newModel = extractAndAddMiriam(newModel,header(miriamCols),data(:,miriamCols),newMets,addedIdx,'met'); end end %% Genes if ~strcmp(genesInfo,'none') - % Gather all data, first about reaction stoichiometries... - fid = fopen(genesInfo); - raw = textscan(fid,['%q' repmat(' %q',1,10)],'Delimiter','\t'); - fclose(fid); - genesToAdd.genes = raw{1}(2:end); - genesToAdd.geneShortNames = raw{2}(2:end); + [header,data] = readCurationTable(genesInfo); + [col,miriamCols] = matchColumns(header,genesInfo,{'genes'},{ ... + 'genes', {'genes','gene'}; ... + 'geneShortNames', {'geneShortNames','geneShortName'}}); + genesToAdd.genes = data(:,col.genes); + if isfield(col,'geneShortNames') + genesToAdd.geneShortNames = data(:,col.geneShortNames); + end [~,notNewGene,existingGene] = intersect(genesToAdd.genes,newModel.genes); if ~isempty(notNewGene) warning(['The following genes are already present in the model, their annotation '... 'will be overwritten to match the genesInfo file: \n\t\t%s'],... strjoin(genesToAdd.genes(notNewGene),'\n\t\t')); - if isfield(newModel,'geneShortNames') - newModel.geneShortNames(existingGene) = genesToAdd.geneShortNames(notNewGene); + if isfield(genesToAdd,'geneShortNames') + newModel = writeExistingField(newModel,'geneShortNames',existingGene,genesToAdd.geneShortNames(notNewGene)); end - newModel = extractAndAddMiriam(newModel,raw(3:end),notNewGene,existingGene,'gene'); + newModel = extractAndAddMiriam(newModel,header(miriamCols),data(:,miriamCols),notNewGene,existingGene,'gene'); end - + % Continue with new genes toAdd = 1:numel(genesToAdd.genes); toAdd(notNewGene) = []; if ~isempty(toAdd) genesToAdd.genes = genesToAdd.genes(toAdd); - genesToAdd.geneShortNames = genesToAdd.geneShortNames(toAdd); - if all(cellfun(@isempty,genesToAdd.geneShortNames)) - genesToAdd = rmfield(genesToAdd,'geneShortNames'); + if isfield(genesToAdd,'geneShortNames') + genesToAdd.geneShortNames = genesToAdd.geneShortNames(toAdd); + if all(cellfun(@isempty,genesToAdd.geneShortNames)) + genesToAdd = rmfield(genesToAdd,'geneShortNames'); + end end % Add genes newModel = addGenesRaven(newModel,genesToAdd); addedIdx = numel(newModel.genes)-numel(toAdd)+1:numel(newModel.genes); - newModel = extractAndAddMiriam(newModel,raw(3:end),toAdd,addedIdx,'gene'); + newModel = extractAndAddMiriam(newModel,header(miriamCols),data(:,miriamCols),toAdd,addedIdx,'gene'); end end %% Reactions if ~any(strcmp({rxnsCoeffs,rxnsInfo},'none')) % Gather all data, first about reaction stoichiometries... - fid = fopen(rxnsCoeffs); - raw = textscan(fid,'%q %q %q %q %f','Delimiter','\t','HeaderLines',1); - fclose(fid); - rxnCheck.coeffsIdx = str2double(raw{1}); - rxns.rxnNames = raw{2}; - rxns.metNames = raw{3}; - rxns.comps = raw{4}; - rxns.coeffs = raw{5}; - rxnCheck.coeffs = strcat(raw{1},'***',rxns.rxnNames); + [header,data] = readCurationTable(rxnsCoeffs); + col = matchColumns(header,rxnsCoeffs,{'rxnIdx','rxnNames','metNames','comps','coefficient'},{ ... + 'rxnIdx', {'rxnIdx','index'}; ... + 'rxnNames', {'rxnNames','rxnName'}; ... + 'metNames', {'metNames','metName'}; ... + 'comps', {'comps','comp'}; ... + 'coefficient', {'coefficient','coefficients','coeffs','stoichCoeffs'}}); + rxnCheck.coeffsIdx = str2double(data(:,col.rxnIdx)); + rxns.rxnNames = data(:,col.rxnNames); + rxns.metNames = data(:,col.metNames); + rxns.comps = data(:,col.comps); + rxns.coeffs = str2double(data(:,col.coefficient)); + rxnCheck.coeffs = strcat(data(:,col.rxnIdx),'***',rxns.rxnNames); % ... and then additional reaction-specific data. - fid = fopen(rxnsInfo); - raw = textscan(fid,['%q' repmat(' %q',1,20)],'Delimiter','\t'); - fclose(fid); - rxnCheck.rxnsIdx = str2double(raw{1}(2:end)); - rxnCheck.rxns = strcat(raw{1}(2:end),'***',raw{2}(2:end)); - + [header,data] = readCurationTable(rxnsInfo); + [col,miriamCols] = matchColumns(header,rxnsInfo,{'rxnIdx','rxnNames'},{ ... + 'rxnIdx', {'rxnIdx','index'}; ... + 'rxnNames', {'rxnNames','rxnName'}; ... + 'grRules', {'grRules'}; ... + 'lb', {'lb'}; ... + 'ub', {'ub'}; ... + 'rev', {'rev'}; ... + 'subSystems', {'subSystems'}; ... + 'eccodes', {'eccodes'}; ... + 'rxnNotes', {'rxnNotes'}; ... + 'rxnReferences', {'rxnReferences'}; ... + 'rxnConfidenceScores', {'rxnConfidenceScores'}}); + rxnCheck.rxnsIdx = str2double(data(:,col.rxnIdx)); + rxnCheck.rxns = strcat(data(:,col.rxnIdx),'***',data(:,col.rxnNames)); + notMatching = setxor(rxnCheck.coeffs,rxnCheck.rxns); if numel(notMatching)>1 error(['The following reactions ánd/or their indices are not matched '... 'between the rxnsInfo and rxnsCoeffs files:\n\t\t%s'],... strjoin(regexprep(notMatching,'^\d+***',''),'\n\t\t')) end - - rxnsToAdd.rxnNames = raw{2}(2:end); - rxnsToAdd.grRules = raw{3}(2:end); - rxnsToAdd.lb = cellfun(@str2num,raw{4}(2:end)); - rxnsToAdd.ub = cellfun(@str2num,raw{5}(2:end)); - rxnsToAdd.rev = cellfun(@str2num,raw{6}(2:end)); - rxnsToAdd.subSystems = raw{7}(2:end); - rxnsToAdd.eccodes = raw{8}(2:end); - rxnsToAdd.rxnNotes = raw{9}(2:end); - rxnsToAdd.rxnReferences = raw{10}(2:end); - rxnsToAdd.rxnConfidenceScores = cellfun(@str2num,raw{11}(2:end),'UniformOutput',false); - emptyEntry = cellfun(@isempty,rxnsToAdd.rxnConfidenceScores); - if all(emptyEntry) - rxnsToAdd = rmfield(rxnsToAdd,'rxnConfidenceScores'); - elseif any(emptyEntry) % If some rxnConfidenceScores are given, assume 0 for those without rxnConfidenceScores specified - rxnsToAdd.rxnConfidenceScores(emptyEntry) = {0}; - rxnsToAdd.rxnConfidenceScores = cell2mat(rxnsToAdd.rxnConfidenceScores); - else - rxnsToAdd.rxnConfidenceScores = cell2mat(rxnsToAdd.rxnConfidenceScores); + + rxnsToAdd.rxnNames = data(:,col.rxnNames); + for f = intersect({'grRules','subSystems','eccodes','rxnNotes','rxnReferences'},fieldnames(col)','stable') + rxnsToAdd.(f{1}) = data(:,col.(f{1})); end - + for f = intersect({'lb','ub','rev'},fieldnames(col)','stable') + rxnsToAdd.(f{1}) = cellfun(@str2num,data(:,col.(f{1}))); + end + if isfield(col,'rxnConfidenceScores') + rxnsToAdd.rxnConfidenceScores = cellfun(@str2num,data(:,col.rxnConfidenceScores),'UniformOutput',false); + emptyEntry = cellfun(@isempty,rxnsToAdd.rxnConfidenceScores); + if all(emptyEntry) + rxnsToAdd = rmfield(rxnsToAdd,'rxnConfidenceScores'); + else % If some rxnConfidenceScores are given, assume 0 for those without rxnConfidenceScores specified + rxnsToAdd.rxnConfidenceScores(emptyEntry) = {0}; + rxnsToAdd.rxnConfidenceScores = cell2mat(rxnsToAdd.rxnConfidenceScores); + end + end + optFields = intersect({'subSystems','eccodes','rxnNotes','rxnReferences','rxnConfidenceScores'},fieldnames(rxnsToAdd)','stable'); + existingRxn=[]; notNewRxn=[]; for i=1:numel(rxnCheck.rxnsIdx) @@ -272,7 +291,7 @@ notNewRxn = [notNewRxn,i]; end end - + % If some reactions already existed, then replace its information if ~isempty(existingRxn) warning(['The following reactions are the same as existing model '... @@ -280,84 +299,126 @@ 'rxnCoeffs file: \n\t\t(Existing reaction ID): New reaction name\n\t\t%s'],... strjoin(strcat('(',newModel.rxns(existingRxn),{'): '},rxnsToAdd.rxnNames(notNewRxn)),'\n\t\t')); newModel.rxnNames(existingRxn) = rxnsToAdd.rxnNames(notNewRxn); - newModel.lb(existingRxn) = rxnsToAdd.lb(notNewRxn); - newModel.ub(existingRxn) = rxnsToAdd.ub(notNewRxn); - newModel.rev(existingRxn) = rxnsToAdd.rev(notNewRxn); - if isfield(newModel,'subSystems') - newModel.subSystems(existingRxn) = rxnsToAdd.subSystems(notNewRxn); + for f = intersect({'lb','ub','rev'},fieldnames(rxnsToAdd)','stable') + newModel.(f{1})(existingRxn) = rxnsToAdd.(f{1})(notNewRxn); end - if isfield(newModel,'eccodes') - newModel.eccodes(existingRxn) = rxnsToAdd.eccodes(notNewRxn); - end - if isfield(newModel,'rxnNotes') - newModel.rxnNotes(existingRxn) = rxnsToAdd.rxnNotes(notNewRxn); + for f = optFields + values = rxnsToAdd.(f{1})(notNewRxn); + if strcmp(f{1},'subSystems') && isfield(newModel,'subSystems') ... + && any(cellfun(@iscell,newModel.subSystems)) + values = cellfun(@(x) {x}, values, 'UniformOutput', false); + end + newModel = writeExistingField(newModel,f{1},existingRxn,values); end - if isfield(newModel,'rxnConfidenceScores') && isfield(rxnsToAdd,'rxnConfidenceScores') - newModel.rxnConfidenceScores(existingRxn) = rxnsToAdd.rxnConfidenceScores(notNewRxn); + if isfield(rxnsToAdd,'grRules') + emptyEntries = cellfun(@isempty,rxnsToAdd.grRules(notNewRxn)); + if ~all(emptyEntries) + newModel = changeGrRules(newModel,newModel.rxns(existingRxn(~emptyEntries)),rxnsToAdd.grRules(notNewRxn(~emptyEntries)),true); + end end - emptyEntries = cellfun(@isempty,rxnsToAdd.grRules); - if ~all(emptyEntries) - newModel = changeGrRules(newModel,newModel.rxns(existingRxn(~emptyEntries)),rxnsToAdd.grRules(notNewRxn(~emptyEntries)),true); - end - newModel = extractAndAddMiriam(newModel,raw(11:end),notNewRxn,existingRxn,'rxn'); + newModel = extractAndAddMiriam(newModel,header(miriamCols),data(:,miriamCols),notNewRxn,existingRxn,'rxn'); end - + % Continue with new reactions toAdd = 1:numel(rxnsToAdd.rxnNames); toAdd(notNewRxn) = []; if ~isempty(toAdd) rxnsToAdd.rxnNames = rxnsToAdd.rxnNames(toAdd); - rxnsToAdd.lb = rxnsToAdd.lb(toAdd); - rxnsToAdd.ub = rxnsToAdd.ub(toAdd); - rxnsToAdd.rev = rxnsToAdd.rev(toAdd); - rxnsToAdd.subSystems = rxnsToAdd.subSystems(toAdd); - if all(cellfun(@isempty,rxnsToAdd.subSystems)) - rxnsToAdd = rmfield(rxnsToAdd,'subSystems'); - end - rxnsToAdd.eccodes = rxnsToAdd.eccodes(toAdd); - if all(cellfun(@isempty,rxnsToAdd.eccodes)) - rxnsToAdd = rmfield(rxnsToAdd,'eccodes'); - end - rxnsToAdd.rxnNotes = rxnsToAdd.rxnNotes(toAdd); - if all(cellfun(@isempty,rxnsToAdd.rxnNotes)) - rxnsToAdd = rmfield(rxnsToAdd,'rxnNotes'); - end - if isfield(rxnsToAdd,'rxnConfidenceScores') - rxnsToAdd.rxnConfidenceScores = rxnsToAdd.rxnConfidenceScores(toAdd); + for f = intersect({'lb','ub','rev','grRules'},fieldnames(rxnsToAdd)','stable') + rxnsToAdd.(f{1}) = rxnsToAdd.(f{1})(toAdd); end - rxnsToAdd.rxnReferences = rxnsToAdd.rxnReferences(toAdd); - if all(cellfun(@isempty,rxnsToAdd.rxnReferences)) - rxnsToAdd = rmfield(rxnsToAdd,'rxnReferences'); + for f = optFields + rxnsToAdd.(f{1}) = rxnsToAdd.(f{1})(toAdd); + if iscell(rxnsToAdd.(f{1})) && all(cellfun(@isempty,rxnsToAdd.(f{1}))) + rxnsToAdd = rmfield(rxnsToAdd,f{1}); + end end - rxnsToAdd.grRules = rxnsToAdd.grRules(toAdd); rxnsToAdd.mets = rxnsToAdd.mets(toAdd); rxnsToAdd.stoichCoeffs = rxnsToAdd.stoichCoeffs(toAdd); rxnsToAdd.rxns = generateNewIds(model,'rxns',rxnPrefix,numel(rxnsToAdd.rxnNames)); - + newModel = addRxns(newModel,rxnsToAdd,1,[],false,false); rxnsModelIdx = numel(newModel.rxns)-numel(toAdd)+1:numel(newModel.rxns); - newModel = extractAndAddMiriam(newModel,raw(11:end),toAdd,rxnsModelIdx,'rxn'); + newModel = extractAndAddMiriam(newModel,header(miriamCols),data(:,miriamCols),toAdd,rxnsModelIdx,'rxn'); end end end -function newModel = extractAndAddMiriam(model,raw,inputIndex,modelIndex,type) -newModel=model; -miriamName = cell(numel(inputIndex),1); -miriamValues = cell(1,numel(inputIndex)); -for i=1:numel(raw) - miriamName{i} = raw{i}{1}; - miriamValues(1:numel(inputIndex),i) = raw{i}(inputIndex+1); +function [header,data] = readCurationTable(fileName) +% Read a tab-delimited curation table into its header row and a cell +% array of text with one column per header field. +fid = fopen(fileName); +if fid==-1 + error('Could not open the file "%s"',fileName); +end +headerLine = fgetl(fid); +if ~ischar(headerLine) + fclose(fid); + error('The file "%s" is empty',fileName); +end +headerLine = strrep(headerLine,char(65279),''); % byte order mark +header = strtrim(strsplit(headerLine,'\t')); +header = regexprep(header,'^"(.*)"$','$1'); +nCols = numel(header); +raw = textscan(fid,repmat('%q',1,nCols),'Delimiter','\t'); +fclose(fid); +nRows = max(cellfun(@numel,raw)); +data = repmat({''},nRows,nCols); +for i=1:nCols + data(1:numel(raw{i}),i) = raw{i}; +end end -emptyMiriam = all(cellfun(@isempty,miriamValues),1); + +function [col,miriamCols] = matchColumns(header,fileName,required,spec) +% Locate the core columns of a curation table by header name (case +% insensitive). spec has one row per core field: {field, accepted headers}. +% Every other column is a MIRIAM namespace. +col = struct(); +isCore = false(1,numel(header)); +for i=1:size(spec,1) + match = find(ismember(lower(header),lower(spec{i,2}))); + if ~isempty(match) + col.(spec{i,1}) = match(1); + isCore(match) = true; + end +end +missing = setdiff(required,fieldnames(col)','stable'); +if ~isempty(missing) + error('The file "%s" has no column with header: %s',fileName,strjoin(missing,', ')); +end +miriamCols = find(~isCore & ~cellfun(@isempty,header)); +end + +function model = writeExistingField(model,field,modelIdx,values) +% Write the values of an optional field onto existing entities. A field +% the model does not have is created from its ravenModelFields default, +% unless all supplied values are empty. +if ~isfield(model,field) + if iscell(values) && all(cellfun(@isempty,values)) + return + end + reg = ravenModelFields(); + entry = reg(strcmp({reg.name},field)); + n = numel(model.([entry.type 's'])); + if iscell(values) + model.(field) = repmat({entry.default},n,1); + else + model.(field) = repmat(entry.default,n,1); + end +end +model.(field)(modelIdx) = values; +end + +function newModel = extractAndAddMiriam(model,miriamName,miriamData,inputIndex,modelIndex,type) +newModel=model; +miriamValues = miriamData(inputIndex,:); +emptyMiriam = all(cellfun(@isempty,miriamValues),1); miriamName(emptyMiriam) = []; miriamValues(:,emptyMiriam) = []; -if ~isfield(newModel,[type 'Miriams']); +if ~isfield(newModel,[type 'Miriams']) newModel.([type 'Miriams'])=cell(numel(newModel.([type 's'])),1); end -if ~isempty(miriamName) - for i=1:numel(miriamName) - newModel = editMiriam(newModel,type,modelIndex,miriamName{i},miriamValues(:,i),'replace'); - end +for i=1:numel(miriamName) + newModel = editMiriam(newModel,type,modelIndex,miriamName{i},miriamValues(:,i),'replace'); end end diff --git a/tests/function_tests/tCuration.m b/tests/function_tests/tCuration.m index ae6327af..845dc6c5 100644 --- a/tests/function_tests/tCuration.m +++ b/tests/function_tests/tCuration.m @@ -39,6 +39,80 @@ function curateModelFromTablesOverwritesExactDuplicate(testCase) testCase.verifyEqual(out.rxnNames{1}, 'newRxn'); end + function curateModelFromTablesReadsMetColumnsByHeader(testCase) + % Columns are found by header, so a table without an inchi + % column and in another order must not write its notes over + % the InChIs; notes and charges reach existing metabolites, + % creating the fields the model lacks. + model = tCuration.toyModel(); + model.inchis = {'InChI=A';'InChI=B';'InChI=C';'InChI=H'}; + metsF = tCuration.writeTsv(testCase, { ... + 'comps\tmetNames\tcharge\tmetNotes\tformula', ... + 'c\tA\t-2\tnote A\tC1'}); + evalc('out = curateModelFromTables(model,metsF);'); + + testCase.verifyNumElements(out.mets, 4); + testCase.verifyEqual(out.inchis, model.inchis); + testCase.verifyEqual(out.metFormulas{1}, 'C1'); + testCase.verifyEqual(out.metNotes, {'note A';'';'';''}); + testCase.verifyEqual(out.metCharges(1), -2); + testCase.verifyTrue(all(isnan(out.metCharges(2:4)))); + end + + function curateModelFromTablesWritesExistingRxnFields(testCase) + % Notes, references, confidence scores and MIRIAM columns + % reach an existing reaction, and grRules are matched to the + % right reaction when existing and new reactions are mixed. + % The yeast-GEM tables name the index column "index". + model = tCuration.toyModel(); + model.genes = {'g1'}; + model.rxnGeneMat = sparse(1,1); + coeffsF = tCuration.writeTsv(testCase, { ... + 'index\trxnNames\tmetNames\tcomps\tcoefficient', ... + '1\tdupRxn\tA\tc\t-1', '1\tdupRxn\tB\tc\t-1', ... + '1\tdupRxn\tC\tc\t1', '1\tdupRxn\tH\tc\t1', ... + '2\tnewRxn\tA\tc\t-1', '2\tnewRxn\tB\tc\t-1', ... + '2\tnewRxn\tC\tc\t1'}); + infoF = tCuration.writeTsv(testCase, { ... + ['index\trxnNames\tgrRules\tlb\tub\trev\tsubSystems\t' ... + 'eccodes\trxnNotes\trxnReferences\trxnConfidenceScores\tkegg.reaction'], ... + '1\tdupRxn\t\t0\t1000\t0\t\t\ta note\tPMID:1\t3\tR00001', ... + '2\tnewRxn\tg1\t0\t1000\t0\t\t\t\t\t\t'}); + evalc(['out = curateModelFromTables(model,''none'',''rxnsCoeffs'',' ... + 'coeffsF,''rxnsInfo'',infoF);']); + + testCase.verifyNumElements(out.rxns, 2); + testCase.verifyEqual(out.rxnNames{1}, 'dupRxn'); + testCase.verifyEqual(out.rxnNotes{1}, 'a note'); + testCase.verifyEqual(out.rxnReferences{1}, 'PMID:1'); + testCase.verifyEqual(out.rxnConfidenceScores(1), 3); + testCase.verifyEqual(out.rxnMiriams{1}.name, {'kegg.reaction'}); + testCase.verifyEqual(out.rxnMiriams{1}.value, {'R00001'}); + testCase.verifyEqual(out.grRules{1}, ''); + testCase.verifyEqual(out.grRules{2}, 'g1'); + end + + function curateModelFromTablesWritesExistingGeneFields(testCase) + % yeast-GEM gene tables use the header geneShortName; it must + % be read as the short name, not as a MIRIAM namespace. + model = tCuration.toyModel(); + model.genes = {'g1'}; + model.rxnGeneMat = sparse(1,1); + genesF = tCuration.writeTsv(testCase, { ... + 'genes\tgeneShortName\tuniprot', 'g1\tGENE1\tP12345'}); + evalc('out = curateModelFromTables(model,''none'',''genesInfo'',genesF);'); + + testCase.verifyEqual(out.geneShortNames, {'GENE1'}); + testCase.verifyEqual(out.geneMiriams{1}.name, {'uniprot'}); + testCase.verifyEqual(out.geneMiriams{1}.value, {'P12345'}); + end + + function curateModelFromTablesRequiresKeyColumns(testCase) + model = tCuration.toyModel(); + metsF = tCuration.writeTsv(testCase, {'metNames\tformula', 'A\tC1'}); + testCase.verifyError(@() curateModelFromTables(model,metsF), ?MException); + end + end methods (Static, Access = private) @@ -92,5 +166,17 @@ function curateModelFromTablesOverwritesExactDuplicate(testCase) fclose(fid); end + function f = writeTsv(testCase, lines) + % Write tab-delimited lines (with \t escapes) to a temporary file. + d = tempname; mkdir(d); + testCase.addTeardown(@() rmdir(d,'s')); + f = fullfile(d,'table.tsv'); + fid = fopen(f,'w'); + for i=1:numel(lines) + fprintf(fid,[lines{i} '\n']); + end + fclose(fid); + end + end end From 53f4a7083f85ea888259947985d21e513592dd0b Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:52:22 +0200 Subject: [PATCH 10/14] docs: state how haveFlux treats a dead-end reaction with a forced flux --- src/analysis/haveFlux.m | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/analysis/haveFlux.m b/src/analysis/haveFlux.m index aa0ae4f0..fc04d630 100755 --- a/src/analysis/haveFlux.m +++ b/src/analysis/haveFlux.m @@ -39,6 +39,12 @@ % ----- % If a model has +/- Inf bounds then those are replaced with an arbitary % large value of +/- 10000 prior to solving. +% +% Dead-end reactions are removed (simplifyModel) before solving. A dead-end +% reaction with a forced flux (lb > 0 or ub < 0) makes the model infeasible, +% but it is removed like any other dead end, so the remaining reactions are +% tested as if it were absent. A forced flux elsewhere in an infeasible +% model gives false for all reactions. p=parseRAVENargs(varargin, {'cutOff',10^-6; 'rxns',[]; 'seed',[]}); if ~isempty(p.seed) From e584572f1e69cc119d7995378e37b8c264d59b2a Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:53:07 +0200 Subject: [PATCH 11/14] fix: keep the lower bound of unreversed reaction candidates in gapFillMILP Every draft reaction with lb >= 0 is a reversal candidate, and its lower bound was relaxed to -bigM with a coupling that only kept the flux non-negative while the reaction was not reversed. A forced flux such as a maintenance requirement (lb > 0) was therefore dropped, and the MILP could report a repair, or no repair, that does not reach the growth floor in the returned model. The coupling now keeps the original lower bound unless the reaction is reversed. --- src/gapfilling/gapFillMILP.m | 11 ++++++---- tests/function_tests/tGapfilling.m | 35 ++++++++++++++++++++++++++++++ 2 files changed, 42 insertions(+), 4 deletions(-) diff --git a/src/gapfilling/gapFillMILP.m b/src/gapfilling/gapFillMILP.m index e66d228c..ac3094f0 100644 --- a/src/gapfilling/gapFillMILP.m +++ b/src/gapfilling/gapFillMILP.m @@ -206,20 +206,23 @@ [ones(1, nUniv), -uniLb'], ... nUniv, nVar); -% (5) Reversal coupling: v_j + M * y_rev_j >= 0 [nRev rows, 'G'] -% (forces v_j >= 0 when y_rev_j = 0, i.e., preserves original irreversibility) +% (5) Reversal coupling: v_j + (M + lb_j) * y_rev_j >= lb_j [nRev rows, 'G'] +% (forces v_j >= lb_j when y_rev_j = 0, i.e., preserves the original +% lower bound, including a forced flux such as maintenance with lb_j > 0; +% allows v_j >= -M when y_rev_j = 1) yRevCols = nMergedRxns + (1:nRev); +revLb = mergedModel.lb(revCandIdx); A_rev = sparse(... [1:nRev, 1:nRev], ... [revCandIdx(:)', yRevCols], ... - [ones(1, nRev), M * ones(1, nRev)], ... + [ones(1, nRev), M + revLb(:)'], ... nRev, nVar); % --- Assemble --- prob.A = [A_SS; A_growth; A_uniUp; A_uniLo; A_rev]; prob.a = prob.A; prob.b = [zeros(nMergedMets, 1); minGrowth; ... - zeros(nUniv, 1); zeros(nUniv, 1); zeros(nRev, 1)]; + zeros(nUniv, 1); zeros(nUniv, 1); revLb(:)]; prob.csense = [repmat('E', 1, nMergedMets), 'G', ... repmat('L', 1, nUniv), ... repmat('G', 1, nUniv + nRev)]; diff --git a/tests/function_tests/tGapfilling.m b/tests/function_tests/tGapfilling.m index 7e6e8de2..8f5a2386 100644 --- a/tests/function_tests/tGapfilling.m +++ b/tests/function_tests/tGapfilling.m @@ -390,6 +390,41 @@ function gapFillMILPKeepsRenamedUniversalReaction(testCase) testCase.verifyGreaterThan(sol.f, 0); end + function gapFillMILPKeepsForcedFlux(testCase) + % A draft reaction with a forced flux (lb > 0, like a + % maintenance drain) is a reversal candidate, but unless it is + % reversed its lower bound must hold. Maint takes 5 of the 10 + % units of a, so growth 8 needs the extra uptake U1; dropping + % the forced flux would reach 8 with no repair at all. + evalc('setRavenSolver(''glpk'')'); + warnState = warning('off','RAVEN:glpkMILP'); + testCase.addTeardown(@() warning(warnState)); + gapModel = struct(); + gapModel.id='gapModel'; gapModel.rxns={'ExA';'R1';'Grow';'Maint'}; gapModel.rxnNames=gapModel.rxns; + gapModel.mets={'a';'b'}; gapModel.metNames=gapModel.mets; gapModel.metComps=[1;1]; + gapModel.comps={'c'}; gapModel.compNames={'c'}; + gapModel.S=sparse([1 -1 0 -1; 0 1 -1 0]); % ExA: =>a R1: a=>b Grow: b=> Maint: a=> + gapModel.lb=[0;0;0;5]; gapModel.ub=[10;1000;1000;1000]; gapModel.rev=[0;0;0;0]; + gapModel.c=[0;0;1;0]; gapModel.b=zeros(2,1); + gapModel.genes={}; gapModel.grRules={'';'';'';''}; gapModel.rxnGeneMat=sparse(4,0); + + modelDB = struct(); + modelDB.id='DB'; modelDB.rxns={'U1'}; modelDB.rxnNames=modelDB.rxns; + modelDB.mets={'a'}; modelDB.metNames=modelDB.mets; modelDB.metComps=1; + modelDB.comps={'c'}; modelDB.compNames={'c'}; + modelDB.S=sparse(1,1); modelDB.S(1,1)=1; % =>a + modelDB.lb=0; modelDB.ub=10; modelDB.rev=0; modelDB.c=0; modelDB.b=0; + modelDB.genes={}; modelDB.grRules={''}; modelDB.rxnGeneMat=sparse(1,0); + + evalc(['[addedRxns,reversedRxns,newModel,exitFlag] = gapFillMILP(gapModel, modelDB, ' ... + '''minGrowth'', 8, ''weights'', [10 2 3 4], ''verbose'', false);']); + testCase.verifyEqual(exitFlag, 1); + testCase.verifyEqual(addedRxns, {'U1'}); + testCase.verifyEmpty(reversedRxns); + sol = solveLP(newModel); + testCase.verifyGreaterThanOrEqual(sol.f, 8 - 1e-6); + end + function gapFillMILPReversesDirectionality(testCase) % gapFillMILP should reverse a reaction whose directionality is wrong. testCase.assumeMILPSolver(); From faf33c7bd00c510e145e7253427acfa04d376176 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 09:11:02 +0200 Subject: [PATCH 12/14] fix: create rxnNotes in getKEGGModelForOrganism when the reference model has none The kegg118 reference model has no reaction notes, so readYAMLmodel returns no rxnNotes field and adding the provenance note failed. --- src/reconstruction/kegg/getKEGGModelForOrganism.m | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/reconstruction/kegg/getKEGGModelForOrganism.m b/src/reconstruction/kegg/getKEGGModelForOrganism.m index de8ce653..18330504 100755 --- a/src/reconstruction/kegg/getKEGGModelForOrganism.m +++ b/src/reconstruction/kegg/getKEGGModelForOrganism.m @@ -340,7 +340,11 @@ model.geneMiriams{i,1}.name{1,1}='kegg.genes'; model.geneMiriams{i,1}.value{1,1}=strcat(lower(organismID),model.genes{i,1}); end - %Add the description to the reactions + %Add the description to the reactions. The reference model has no + %reaction notes, so readYAMLmodel leaves the field out + if ~isfield(model,'rxnNotes') + model.rxnNotes=repmat({''},numel(model.rxns),1); + end for i=1:numel(model.rxns) if ~isempty(model.rxnNotes{i}) model.rxnNotes(i)=strcat('Included by getKEGGModelForOrganism (without HMMs).',model.rxnNotes(i)); @@ -532,7 +536,11 @@ model.grRules = grRules; model.rxnGeneMat = rxnGeneMat; -%Add the description to the reactions +%Add the description to the reactions. The reference model has no reaction +%notes, so readYAMLmodel leaves the field out +if ~isfield(model,'rxnNotes') + model.rxnNotes=repmat({''},numel(model.rxns),1); +end for i=1:numel(model.rxns) if ~isempty(model.rxnNotes{i}) model.rxnNotes(i)=strcat('Included by getKEGGModelForOrganism (using HMMs).',model.rxnNotes(i)); From 1f336ed3b97c7dd9069b019f1e9c6ff1ad340c12 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 09:27:02 +0200 Subject: [PATCH 13/14] fix: correct a typo in the curateModelFromTables overwrite warning --- src/curation/curateModelFromTables.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/curation/curateModelFromTables.m b/src/curation/curateModelFromTables.m index 47c6d976..2b0e2d0d 100644 --- a/src/curation/curateModelFromTables.m +++ b/src/curation/curateModelFromTables.m @@ -136,7 +136,7 @@ warning(['The following metabolites are already present in the model, '... 'their annotation will be overwritten to match the metsInfo file. '... 'If you do not particularly want to curate their annotations, it '... - 'would be better to removes these metabolites from metsInfo:\n\t\t%s'],... + 'would be better to remove these metabolites from metsInfo:\n\t\t%s'],... strjoin(strcat(metsToAdd.metNames(existingMets),'[',metsToAdd.compartments(existingMets),']'),'\n\t\t')); end From d22836484e498060a7cff10874039880311cb24e Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 13:25:26 +0200 Subject: [PATCH 14/14] gapFillMILP adds database reactions only; gapReport stops on an infeasible model gapFillMILP no longer reverses draft reactions: the draft keeps its bounds and the MILP minimises the number of added reactions. The reversedRxns output and the weights parameter are removed. gapReport raises RAVEN:infeasibleModel when the model is infeasible under its own bounds or with its boundary closed, naming the forced-flux reactions. --- src/gapfilling/fillGaps.m | 4 +- src/gapfilling/gapFillMILP.m | 118 ++++++++--------------------- src/gapfilling/gapReport.m | 30 ++++++++ tests/function_tests/tGapfilling.m | 102 ++++++++++++++++++++----- 4 files changed, 147 insertions(+), 107 deletions(-) diff --git a/src/gapfilling/fillGaps.m b/src/gapfilling/fillGaps.m index ee9f13bb..0ebe1605 100755 --- a/src/gapfilling/fillGaps.m +++ b/src/gapfilling/fillGaps.m @@ -55,7 +55,7 @@ % - 'reference' : MILP-based connectivity maximisation % - 'fastLP' : L1-norm LP relaxation via FASTCORE (gapFillFastLP) % - 'swiftLP' : SWIFTCORE single-LP variant (gapFillFastLP) -% - 'gapfillMILP' : growth-floor MILP with directionality repair +% - 'gapfillMILP' : growth-floor MILP adding database reactions % (gapFillMILP). Uses only models{1} as the database. % - 'topological' : BFS metabolite producibility pre-screening % (gapFillTopological); no model modification. @@ -142,7 +142,7 @@ if numel(models) < 1 error('RAVEN:badInput', 'fillGaps: gapfillMILP requires at least one reference model.'); end - [addedRxns, ~, newModel, exitFlag] = gapFillMILP(model, models{1}, ... + [addedRxns, newModel, exitFlag] = gapFillMILP(model, models{1}, ... 'minGrowth', p.minGrowth, 'verbose', p.verbose, 'params', params); newConnected = {}; cannotConnect = {}; diff --git a/src/gapfilling/gapFillMILP.m b/src/gapfilling/gapFillMILP.m index ac3094f0..34da6e53 100644 --- a/src/gapfilling/gapFillMILP.m +++ b/src/gapfilling/gapFillMILP.m @@ -1,18 +1,18 @@ -function [addedRxns, reversedRxns, newModel, exitFlag] = gapFillMILP(model, universalModel, varargin) +function [addedRxns, newModel, exitFlag] = gapFillMILP(model, universalModel, varargin) % gapFillMILP Objective-based gap-filling via a global growth-floor MILP. % -% Finds the minimum-cost set of modifications to make the model's objective -% (biomass) >= minGrowth. Supports two repair mechanisms from Kumar et al. -% 2007 (Bioinformatics 23:1626-1635): -% 1. Directionality reversal of existing draft reactions (binary y_rev) -% 2. Addition of reactions from a universal database (binary y_db) +% Finds the smallest set of reactions from a universal database that makes +% the model's objective (biomass) >= minGrowth, following the database +% repair of Kumar et al. 2007 (BMC Bioinformatics 8:212, +% doi:10.1186/1471-2105-8-212). Reactions already in the draft model keep +% their bounds: a reaction with the wrong direction is a curation issue and +% is not changed here. % % The universal model should contain all candidate reactions, including % exchange and transport reactions if those repair classes are desired. % -% This is a SINGLE-level MILP (not bilevel). It is computationally more -% expensive than gapFillFastLP but supports the directionality-reversal -% repair mechanism and directly optimises for objective feasibility. +% This is a single-level MILP. It is computationally more expensive than +% gapFillFastLP but directly optimises for objective feasibility. % % Parameters % ---------- @@ -26,11 +26,6 @@ % 'minGrowth' (default []) % minimum objective value to achieve. If empty, runs FBA on the merged % model and uses 10% of the maximum as the floor. -% 'weights' (default [1 2 3 4]) -% cost weights [w_rev w_db w_exch w_trans]. w_rev is the cost per -% reversed draft reaction; w_db is the cost per added universal -% reaction. w_exch and w_trans are reserved for future use when -% exchange and transport reactions are tracked separately. % 'bigM' (default 1000) % Big-M constant for coupling constraints. Must be larger than any % expected flux magnitude. @@ -44,11 +39,8 @@ % ------- % addedRxns : cell % reaction IDs from universalModel that are selected by the MILP. -% reversedRxns : cell -% draft reaction IDs whose directionality is reversed by the MILP. % newModel : struct -% model with addedRxns incorporated and reversedRxns bounds updated -% (lb set to -bigM for reversed reactions). +% model with addedRxns incorporated. % exitFlag : double % 1 = optimal, -1 = infeasible, -2 = timeout or other failure. % @@ -65,25 +57,19 @@ % ---- Parse optional arguments ---- p = inputParser(); addParameter(p, 'minGrowth', [], @(x) isempty(x) || isnumeric(x)); -addParameter(p, 'weights', [1 2 3 4], @isnumeric); addParameter(p, 'bigM', 1000, @isnumeric); addParameter(p, 'params', [], @(x) isempty(x) || isstruct(x)); addParameter(p, 'verbose', true, @islogical); parse(p, varargin{:}); minGrowth = p.Results.minGrowth; -weights = p.Results.weights; M = p.Results.bigM; params = p.Results.params; verbose = p.Results.verbose; -w_rev = weights(1); -w_db = weights(2); - % ---- Initialise outputs ---- -addedRxns = {}; -reversedRxns = {}; -newModel = model; -exitFlag = -1; +addedRxns = {}; +newModel = model; +exitFlag = -1; % Check that the model has an objective if all(model.c == 0) @@ -105,16 +91,10 @@ univIdx = find(isUniv); % positions in merged model nUniv = numel(univIdx); -% Reversal candidates: draft reactions that are currently irreversible -% (lb >= 0). Reversing their directionality means allowing negative flux. -revCandMask = isDraft & mergedModel.lb >= 0; -revCandIdx = find(revCandMask); % positions in merged model -nRev = numel(revCandIdx); - if verbose fprintf('gapFillMILP: merged model has %d mets, %d rxns (%d draft, %d universal).\n', ... nMergedMets, nMergedRxns, sum(isDraft), nUniv); - fprintf('gapFillMILP: %d reversal candidates, %d database candidates.\n', nRev, nUniv); + fprintf('gapFillMILP: %d database candidates.\n', nUniv); end % ---- Determine minGrowth ---- @@ -153,35 +133,30 @@ % ---- Build MILP ---- % Variable ordering: % cols 1 : nMergedRxns → v (all fluxes, continuous) -% cols nMergedRxns+1 : +nRev → y_rev (binary, directionality reversal) -% cols nMergedRxns+nRev+1 : end → y_db (binary, universal reactions) -nVar = nMergedRxns + nRev + nUniv; +% cols nMergedRxns+1 : end → y_db (binary, universal reactions) +nVar = nMergedRxns + nUniv; % --- Variable bounds --- lb = mergedModel.lb; ub = mergedModel.ub; -% Relax reversal candidates to allow negative flux (Big-M on lower bound) -lb(revCandIdx) = -M; - % Universal reactions: allow them to be forced to 0 (extend lb toward 0) lb(univIdx) = min(mergedModel.lb(univIdx), 0); % Binary variable bounds -lb = [lb; zeros(nRev + nUniv, 1)]; -ub = [ub; ones(nRev + nUniv, 1)]; +lb = [lb; zeros(nUniv, 1)]; +ub = [ub; ones(nUniv, 1)]; -% --- Objective --- +% --- Objective: number of added universal reactions --- c = zeros(nVar, 1); -c(nMergedRxns + 1 : nMergedRxns + nRev) = w_rev; % y_rev cost -c(nMergedRxns + nRev + 1 : end) = w_db; % y_db cost +c(nMergedRxns + 1 : end) = 1; % --- Constraint blocks --- % (1) Steady state: S * v = 0 [nMergedMets rows, 'E'] -A_SS = [mergedModel.S, sparse(nMergedMets, nRev + nUniv)]; +A_SS = [mergedModel.S, sparse(nMergedMets, nUniv)]; % (2) Growth floor: c' * v >= minGrowth [1 row, 'G'] -A_growth = [sparse(mergedModel.c'), sparse(1, nRev + nUniv)]; +A_growth = [sparse(mergedModel.c'), sparse(1, nUniv)]; % (3) Universal upper coupling: v_k - ub_k * y_db_k <= 0 [nUniv rows, 'L'] % (forces v_k <= 0 when y_db_k = 0) @@ -189,7 +164,7 @@ % universal reaction (ub = Inf, RAVEN's own default when no upper-bound % annotation is set) would otherwise put an Inf coefficient into the % constraint matrix, and Inf*0 (y_db_k = 0) is NaN, not 0. -yDbCols = nMergedRxns + nRev + (1:nUniv); +yDbCols = nMergedRxns + (1:nUniv); uniUb = min(mergedModel.ub(univIdx), M); uniLb = max(mergedModel.lb(univIdx), -M); A_uniUp = sparse(... @@ -206,32 +181,19 @@ [ones(1, nUniv), -uniLb'], ... nUniv, nVar); -% (5) Reversal coupling: v_j + (M + lb_j) * y_rev_j >= lb_j [nRev rows, 'G'] -% (forces v_j >= lb_j when y_rev_j = 0, i.e., preserves the original -% lower bound, including a forced flux such as maintenance with lb_j > 0; -% allows v_j >= -M when y_rev_j = 1) -yRevCols = nMergedRxns + (1:nRev); -revLb = mergedModel.lb(revCandIdx); -A_rev = sparse(... - [1:nRev, 1:nRev], ... - [revCandIdx(:)', yRevCols], ... - [ones(1, nRev), M + revLb(:)'], ... - nRev, nVar); - % --- Assemble --- -prob.A = [A_SS; A_growth; A_uniUp; A_uniLo; A_rev]; +prob.A = [A_SS; A_growth; A_uniUp; A_uniLo]; prob.a = prob.A; -prob.b = [zeros(nMergedMets, 1); minGrowth; ... - zeros(nUniv, 1); zeros(nUniv, 1); revLb(:)]; +prob.b = [zeros(nMergedMets, 1); minGrowth; zeros(2 * nUniv, 1)]; prob.csense = [repmat('E', 1, nMergedMets), 'G', ... repmat('L', 1, nUniv), ... - repmat('G', 1, nUniv + nRev)]; + repmat('G', 1, nUniv)]; prob.c = c; -prob.osense = 1; % minimise weighted sum of binary variables +prob.osense = 1; % minimise the number of added reactions prob.lb = lb; prob.ub = ub; prob.vartype = repmat('C', 1, nVar); -prob.vartype(nMergedRxns + 1 : end) = 'B'; % y_rev and y_db are binary +prob.vartype(nMergedRxns + 1 : end) = 'B'; % y_db is binary if isempty(params) params.intTol = 1e-9; @@ -241,7 +203,7 @@ % ---- Solve ---- if verbose fprintf('gapFillMILP: solving MILP (%d vars, %d binary, %d constraints)...\n', ... - nVar, nRev + nUniv, size(prob.A, 1)); + nVar, nUniv, size(prob.A, 1)); end sol = optimizeProb(prob, params, verbose); @@ -255,11 +217,8 @@ % ---- Extract solution ---- binThreshold = 1e-3; -yRev = sol.full(nMergedRxns + 1 : nMergedRxns + nRev); -yDb = sol.full(nMergedRxns + nRev + 1 : end); - -revSelected = revCandIdx(yRev > binThreshold); -dbSelected = univIdx(yDb > binThreshold); +yDb = sol.full(nMergedRxns + 1 : end); +dbSelected = univIdx(yDb > binThreshold); % Map back to original IDs. dbSelected is already restricted to universal % reactions via univIdx/isDraft above (an id-membership test against @@ -267,26 +226,15 @@ % further filtering against universalModel.rxns is needed here -- doing so % would silently drop a universal reaction that mergeModels renamed % because its id collided with one already in the draft model. -reversedRxns = mergedModel.rxns(revSelected); -addedRxns = mergedModel.rxns(dbSelected); +addedRxns = mergedModel.rxns(dbSelected); if verbose - fprintf('gapFillMILP: reversed %d draft reaction(s), added %d universal reaction(s).\n', ... - numel(reversedRxns), numel(addedRxns)); + fprintf('gapFillMILP: added %d universal reaction(s).\n', numel(addedRxns)); end % ---- Build new model ---- newModel = model; -% Apply directionality reversals: allow negative flux -for k = 1:numel(reversedRxns) - rxnPos = find(getIndexes(newModel, reversedRxns(k), 'rxns', true)); - if ~isempty(rxnPos) - newModel.lb(rxnPos) = -M; - newModel.rev(rxnPos) = 1; - end -end - % Add universal reactions. Extracted from mergedModel by the dbSelected % indices directly, not by looking addedRxns' names back up in % universalModel: that lookup fails whenever mergeModels renamed one of diff --git a/src/gapfilling/gapReport.m b/src/gapfilling/gapReport.m index b1816068..b28bdc38 100755 --- a/src/gapfilling/gapReport.m +++ b/src/gapfilling/gapReport.m @@ -54,6 +54,14 @@ % the reactions that were added from the template models and which % model they were added from. % +% Notes +% ----- +% The model must be feasible under its own bounds, and, when the mass +% balancing section runs, also with its boundary closed. Otherwise +% gapReport stops with a RAVEN:infeasibleModel error that names the +% reactions with a forced flux; set the bounds of ATP maintenance and any +% other forced flux to zero first. +% % Examples % -------- % [noFluxRxns, noFluxRxnsRelaxed, subGraphs, notProducedMets, ... @@ -85,6 +93,11 @@ calculateINOUT=false; end +requireFeasible(model,'under its own bounds'); +if calculateINOUT + requireFeasible(closedModel,'once its boundary reactions are closed'); +end + model2=model; model2.b=[model2.b inf(numel(model2.mets),1)]; I=haveFlux(model); @@ -133,3 +146,20 @@ fprintf([num2str(numel(connectedFromTemplates)) ' unconnected reactions can be connected by including ' num2str(numel(addedFromTemplates)) ' reactions from\n' t '\n']); end end + +function requireFeasible(model,condition) +%Errors, naming the forced-flux reactions, if the model is infeasible +model.c(:)=0; +sol=solveLP(model); +if sol.stat~=-1 + return +end +forced=model.rxns(model.lb>0 | model.ub<0); +hint=''; +if ~isempty(forced) + hint=[' Reactions with a forced flux: ' strjoin(forced,', ') '.']; +end +error('RAVEN:infeasibleModel',['The model is infeasible ' condition ... + ', so gapReport cannot analyse it.' strrep(hint,'%','%%') ... + ' Set the bounds of ATP maintenance and any other forced flux to zero first.']); +end diff --git a/tests/function_tests/tGapfilling.m b/tests/function_tests/tGapfilling.m index 8f5a2386..d08bb14c 100644 --- a/tests/function_tests/tGapfilling.m +++ b/tests/function_tests/tGapfilling.m @@ -80,6 +80,26 @@ function gapReportRuns(testCase) testCase.verifyClass(noFluxRxns, 'cell'); end + function gapReportInfeasibleModelErrors(testCase) + % Maint needs 20 units of a but uptake allows 10. + m = testCase.tinyModel({'ExA';'Maint'}, sparse([1 -1]), [0;20], [10;1000]); + err = testCase.gapReportError(m); + testCase.verifySubstring(err.message, 'under its own bounds'); + testCase.verifySubstring(err.message, 'Maint'); + testCase.verifySubstring(err.message, 'ATP maintenance'); + end + + function gapReportInfeasibleOnceClosedErrors(testCase) + % Uptake of a comes from the unconstrained a_b: feasible with + % the boundary open, but Maint cannot run once it is closed. + m = testCase.tinyModel({'ExA';'Maint'}, sparse([1 -1; -1 0]), [0;1], [10;1000]); + m.mets = {'a';'a_b'}; m.metNames = m.mets; m.metComps = [1;1]; + m.b = zeros(2,1); m.unconstrained = [0;1]; + err = testCase.gapReportError(m); + testCase.verifySubstring(err.message, 'once its boundary reactions are closed'); + testCase.verifySubstring(err.message, 'Maint'); + end + function gapReportKeepsPercentInModelName(testCase) % A model.name containing "%" must survive intact in the % report header, not be truncated by fprintf misreading it as @@ -314,9 +334,8 @@ function gapFillMILPRepairsGrowth(testCase) modelDB = testCase.model; modelDB.id = 'DB'; gapModel = removeReactions(modelDB, modelDB.rxns(1:5)); gapModel.id = 'gapModel'; - evalc('[addedRxns, reversedRxns, newModel, exitFlag] = gapFillMILP(gapModel, modelDB, ''verbose'', false);'); + evalc('[addedRxns, newModel, exitFlag] = gapFillMILP(gapModel, modelDB, ''verbose'', false);'); testCase.verifyClass(addedRxns, 'cell'); - testCase.verifyClass(reversedRxns, 'cell'); testCase.verifyClass(newModel, 'struct'); testCase.verifyEqual(exitFlag, 1); % The repaired model should be able to produce objective flux. @@ -349,7 +368,7 @@ function gapFillMILPHandlesUnboundedUniversalReaction(testCase) modelDB.lb=0; modelDB.ub=Inf; modelDB.rev=0; modelDB.c=0; modelDB.b=0; modelDB.genes={}; modelDB.grRules={''}; modelDB.rxnGeneMat=sparse(1,0); - evalc(['[addedRxns,~,newModel,exitFlag] = gapFillMILP(gapModel, modelDB, ' ... + evalc(['[addedRxns,newModel,exitFlag] = gapFillMILP(gapModel, modelDB, ' ... '''verbose'', false);']); testCase.verifyEqual(exitFlag, 1); testCase.verifyEqual(addedRxns, {'ExA'}); @@ -382,7 +401,7 @@ function gapFillMILPKeepsRenamedUniversalReaction(testCase) modelDB.lb=0; modelDB.ub=1000; modelDB.rev=0; modelDB.c=0; modelDB.b=0; modelDB.genes={}; modelDB.grRules={''}; modelDB.rxnGeneMat=sparse(1,0); - evalc(['[addedRxns,~,newModel,exitFlag] = gapFillMILP(gapModel, modelDB, ' ... + evalc(['[addedRxns,newModel,exitFlag] = gapFillMILP(gapModel, modelDB, ' ... '''verbose'', false);']); testCase.verifyEqual(exitFlag, 1); testCase.verifyNotEmpty(addedRxns); @@ -392,10 +411,9 @@ function gapFillMILPKeepsRenamedUniversalReaction(testCase) function gapFillMILPKeepsForcedFlux(testCase) % A draft reaction with a forced flux (lb > 0, like a - % maintenance drain) is a reversal candidate, but unless it is - % reversed its lower bound must hold. Maint takes 5 of the 10 - % units of a, so growth 8 needs the extra uptake U1; dropping - % the forced flux would reach 8 with no repair at all. + % maintenance drain) keeps its lower bound. Maint takes 5 of + % the 10 units of a, so growth 8 needs the extra uptake U1; + % dropping the forced flux would reach 8 with no repair at all. evalc('setRavenSolver(''glpk'')'); warnState = warning('off','RAVEN:glpkMILP'); testCase.addTeardown(@() warning(warnState)); @@ -416,24 +434,44 @@ function gapFillMILPKeepsForcedFlux(testCase) modelDB.lb=0; modelDB.ub=10; modelDB.rev=0; modelDB.c=0; modelDB.b=0; modelDB.genes={}; modelDB.grRules={''}; modelDB.rxnGeneMat=sparse(1,0); - evalc(['[addedRxns,reversedRxns,newModel,exitFlag] = gapFillMILP(gapModel, modelDB, ' ... - '''minGrowth'', 8, ''weights'', [10 2 3 4], ''verbose'', false);']); + evalc(['[addedRxns,newModel,exitFlag] = gapFillMILP(gapModel, modelDB, ' ... + '''minGrowth'', 8, ''verbose'', false);']); testCase.verifyEqual(exitFlag, 1); testCase.verifyEqual(addedRxns, {'U1'}); - testCase.verifyEmpty(reversedRxns); + testCase.verifyEqual(newModel.lb(ismember(newModel.rxns, 'Maint')), 5); sol = solveLP(newModel); testCase.verifyGreaterThanOrEqual(sol.f, 8 - 1e-6); end - function gapFillMILPReversesDirectionality(testCase) - % gapFillMILP should reverse a reaction whose directionality is wrong. - testCase.assumeMILPSolver(); - model = testCase.model; - modelDB = model; modelDB.id = 'DB'; - gapModel = removeReactions(model, model.rxns(1:3)); - gapModel.id = 'gapModel'; - evalc('[~, reversedRxns, ~, ~] = gapFillMILP(gapModel, modelDB, ''verbose'', false);'); - testCase.verifyClass(reversedRxns, 'cell'); + function gapFillMILPKeepsDraftDirections(testCase) + % A draft reaction running in the wrong direction is not + % reversed: with only that repair available the MILP has no + % solution, and the draft bounds are returned unchanged. + evalc('setRavenSolver(''glpk'')'); + warnState = warning('off','RAVEN:glpkMILP'); + testCase.addTeardown(@() warning(warnState)); + gapModel = struct(); + gapModel.id='gapModel'; gapModel.rxns={'ExB';'R1';'Grow'}; gapModel.rxnNames=gapModel.rxns; + gapModel.mets={'a';'b'}; gapModel.metNames=gapModel.mets; gapModel.metComps=[1;1]; + gapModel.comps={'c'}; gapModel.compNames={'c'}; + gapModel.S=sparse([0 -1 -1; 1 1 0]); % ExB: =>b R1: a=>b Grow: a=> + gapModel.lb=[0;0;0]; gapModel.ub=[10;1000;1000]; gapModel.rev=[0;0;0]; + gapModel.c=[0;0;1]; gapModel.b=zeros(2,1); + gapModel.genes={}; gapModel.grRules={'';'';''}; gapModel.rxnGeneMat=sparse(3,0); + + modelDB = struct(); + modelDB.id='DB'; modelDB.rxns={'Other'}; modelDB.rxnNames=modelDB.rxns; + modelDB.mets={'b'}; modelDB.metNames=modelDB.mets; modelDB.metComps=1; + modelDB.comps={'c'}; modelDB.compNames={'c'}; + modelDB.S=sparse(1,1); modelDB.S(1,1)=-1; % b=> + modelDB.lb=0; modelDB.ub=1000; modelDB.rev=0; modelDB.c=0; modelDB.b=0; + modelDB.genes={}; modelDB.grRules={''}; modelDB.rxnGeneMat=sparse(1,0); + + evalc(['[addedRxns,newModel,exitFlag] = gapFillMILP(gapModel, modelDB, ' ... + '''minGrowth'', 1, ''verbose'', false);']); + testCase.verifyEqual(exitFlag, -1); + testCase.verifyEmpty(addedRxns); + testCase.verifyEqual(newModel.lb, gapModel.lb); end function fillGapsDispatchesFastLP(testCase) @@ -470,5 +508,29 @@ function fillGapsDispatchesGapfillMILP(testCase) m.b = [0;0;0]; m.genes = {}; m.grRules = {'';''}; m.rxnGeneMat = sparse(2,0); end + + function m = tinyModel(~, rxns, S, lb, ub) + % One-compartment model with metabolite a and no objective. + n = numel(rxns); + m = struct(); + m.id = 'tiny'; m.name = 'tiny'; + m.rxns = rxns; m.rxnNames = rxns; + m.mets = {'a'}; m.metNames = {'a'}; m.metComps = 1; + m.comps = {'c'}; m.compNames = {'c'}; + m.S = S; m.lb = lb; m.ub = ub; m.rev = zeros(n,1); m.c = zeros(n,1); + m.b = 0; + m.genes = {}; m.grRules = repmat({''},n,1); m.rxnGeneMat = sparse(n,0); + end + + function err = gapReportError(testCase, m) + % The error gapReport raises on m, checked to be RAVEN:infeasibleModel. + err = []; + try + evalc('gapReport(m);'); + catch err + end + testCase.assertNotEmpty(err); + testCase.verifyEqual(err.identifier, 'RAVEN:infeasibleModel'); + end end end