diff --git a/src/gapfilling/fitTasks.m b/src/gapfilling/fitTasks.m index 78868d16..d311cbbd 100755 --- a/src/gapfilling/fitTasks.m +++ b/src/gapfilling/fitTasks.m @@ -262,11 +262,14 @@ end %Add changed bounds if ~isempty(taskStructure(i).changed) - tModel=setParam(tModel,'lb',taskStructure(i).changed,taskStructure(i).LBrxn); - tModel=setParam(tModel,'ub',taskStructure(i).changed,taskStructure(i).UBrxn); + %NaN bounds mean "leave unchanged" + lbSet=~isnan(taskStructure(i).LBrxn); + ubSet=~isnan(taskStructure(i).UBrxn); + tModel=setParam(tModel,'lb',taskStructure(i).changed(lbSet),taskStructure(i).LBrxn(lbSet)); + tModel=setParam(tModel,'ub',taskStructure(i).changed(ubSet),taskStructure(i).UBrxn(ubSet)); if preMerged - tRefModel=setParam(tRefModel,'lb',taskStructure(i).changed,taskStructure(i).LBrxn); - tRefModel=setParam(tRefModel,'ub',taskStructure(i).changed,taskStructure(i).UBrxn); + tRefModel=setParam(tRefModel,'lb',taskStructure(i).changed(lbSet),taskStructure(i).LBrxn(lbSet)); + tRefModel=setParam(tRefModel,'ub',taskStructure(i).changed(ubSet),taskStructure(i).UBrxn(ubSet)); end end diff --git a/src/io/writeYAMLmodel.m b/src/io/writeYAMLmodel.m index 28879551..b9afd250 100755 --- a/src/io/writeYAMLmodel.m +++ b/src/io/writeYAMLmodel.m @@ -382,7 +382,15 @@ function emitMetaField(fid, key, value) if x == floor(x) && abs(x) < 1e16 s = sprintf('%.1f', x); % e.g. 1000 -> "1000.0" else - s = sprintf('%.15g', x); + % Shortest text that reads back as exactly x (15 to 17 significant + % digits), as Python's repr does, so 1/3 is written 0.3333333333333333 + % and a value is never rounded on the way out. + for digits = 15:17 + s = sprintf(['%.' num2str(digits) 'g'], x); + if str2double(s) == x + break + end + end end end @@ -549,7 +557,9 @@ function emitScalarCore(fid, ~, value, ~, bareNumeric) % because the text itself is not valid as a plain block scalar (leading % indicator character, a mid-string ": "/" #", leading/trailing % whitespace, or embedded control characters). Implements the YAML 1.2 -% rules for resolving implicit scalar types and analyzing scalar style. +% rules for resolving implicit scalar types and analyzing scalar style, +% and also quotes the YAML 1.1 booleans y, n, yes, no, on and off in any +% case, which YAML 1.1 loaders such as PyYAML read as booleans. if isempty(s) tf = true; return @@ -557,6 +567,7 @@ function emitScalarCore(fid, ~, value, ~, bareNumeric) resolverPattern = [ ... '^(true|True|TRUE|false|False|FALSE)$|' ... + '^([Yy]|[Nn]|[Yy][Ee][Ss]|[Nn][Oo]|[Oo][Nn]|[Oo][Ff][Ff])$|' ... '^[-+]?[0-9][0-9_]*\.[0-9_]*([eE][-+]?[0-9]+)?$|' ... '^[-+]?[0-9][0-9_]*[eE][-+]?[0-9]+$|' ... '^[-+]?\.[0-9_]+([eE][-+][0-9]+)?$|' ... diff --git a/src/manipulation/removeLowScoreGenes.m b/src/manipulation/removeLowScoreGenes.m index df2378d8..fb9b8c3d 100644 --- a/src/manipulation/removeLowScoreGenes.m +++ b/src/manipulation/removeLowScoreGenes.m @@ -19,8 +19,8 @@ % score will try to be removed. % % If all genes associated with a reaction have a negative score, then the -% least-negative gene will remain; if there is a tie, one will be selected -% at random. +% least-negative gene will remain; if there is a tie, the gene listed first +% in the rule will be kept. % % If a negative-scoring gene is a subunit in a complex, it will not be % removed; however, the entire complex may be removed. See the following @@ -149,7 +149,7 @@ % get IDs and indices of genes involved in rule -ruleGenes = unique(regexp(rule,'[^&|\(\) ]+','match')); +ruleGenes = unique(regexp(rule,'[^&|\(\) ]+','match'),'stable'); [~,geneInd] = ismember(ruleGenes,genes); % rules with one or no genes remain unchanged @@ -164,8 +164,8 @@ scoreMethod = isozymeScoring; negInd = gScores(geneInd) < 0; % NaNs will return false here if all(negInd) - % deterministic least-negative: sort descending; stable sort breaks - % ties by original gene order (left-to-right in the rule string) + % deterministic least-negative: sort descending; the stable sort + % breaks ties by the order of the genes in the rule string [~,sortOrder] = sort(gScores(geneInd),'descend'); maxInd = sortOrder(1); updatedRule = ruleGenes{maxInd}; @@ -179,7 +179,7 @@ end % update ruleGenes and their indices - ruleGenes = unique(regexp(updatedRule,'[^&|\(\) ]+','match')); + ruleGenes = unique(regexp(updatedRule,'[^&|\(\) ]+','match'),'stable'); [~,geneInd] = ismember(ruleGenes,genes); elseif ~contains(rule,'|') % rule contains enzyme complex diff --git a/src/manipulation/replaceMets.m b/src/manipulation/replaceMets.m index bfc31a9f..bcba7d0f 100755 --- a/src/manipulation/replaceMets.m +++ b/src/manipulation/replaceMets.m @@ -79,6 +79,12 @@ error('The to-be-replaced metabolite cannot be found in the model.'); end +% Either of the two rows can be the one that survives the merge: with identifiers the +% replacement stays, by name the replaced metabolite is renamed and stays while the +% replacement's row is deleted below. Both ids are kept so that contractModel is +% restricted to the metabolite that is left. +repMetIds = [repMetIds; model.mets(metIdx)]; + rxnsWithMet = find(model.S(metIdx,:)); if verbose==true fprintf('\n\nThe following reactions contain the to-be-replaced metabolite as reactant:\n') @@ -173,10 +179,13 @@ % This could now have created duplicate reactions. Contract model. % repIdx may point past the end, or at the wrong row, once idxDelete has -% shifted things above; re-resolve the replacement metabolite by id. +% shifted things above; re-resolve the metabolite that is left by id. An empty +% list must not reach contractModel, which would then check every reaction. % distReverse is passed by name rather than positionally as [], which % parseRAVENargs would take literally instead of falling back to its % documented default of true. repIdx = find(ismember(model.mets, repMetIds)); -model=contractModel(model,'mets',repIdx); +if ~isempty(repIdx) + model=contractModel(model,'mets',repIdx); +end end diff --git a/src/tasks/checkTasks.m b/src/tasks/checkTasks.m index e9e309bc..778f4eff 100755 --- a/src/tasks/checkTasks.m +++ b/src/tasks/checkTasks.m @@ -112,6 +112,13 @@ if ~isempty(taskStructure(i).inputs) [I, J]=ismember(upper(taskStructure(i).inputs),modelMets); J=J(I); %Only keep the ones with matches + %Several model metabolites can share a name[comp]. Only the first + %one is constrained, so tell the user when that happens + dupMets=duplicateMetNames(taskStructure(i).inputs(I),modelMets); + if ~isempty(dupMets) + EM=['Several model metabolites match input(s) in "[' taskStructure(i).id '] ' taskStructure(i).description '": ' strjoin(dupMets,', ') '. Only the first match is constrained']; + warning('RAVEN:warning', '%s', EM); + end K=ismember(upper(taskStructure(i).inputs),'ALLMETS'); L=~cellfun('isempty',strfind(upper(taskStructure(i).inputs),'ALLMETSIN')); %Check that all metabolites are either real metabolites or @@ -165,6 +172,13 @@ if ~isempty(taskStructure(i).outputs) [I, J]=ismember(upper(taskStructure(i).outputs),modelMets); J=J(I); %Only keep the ones with matches + %Several model metabolites can share a name[comp]. Only the first + %one is constrained, so tell the user when that happens + dupMets=duplicateMetNames(taskStructure(i).outputs(I),modelMets); + if ~isempty(dupMets) + EM=['Several model metabolites match output(s) in "[' taskStructure(i).id '] ' taskStructure(i).description '": ' strjoin(dupMets,', ') '. Only the first match is constrained']; + warning('RAVEN:warning', '%s', EM); + end K=ismember(upper(taskStructure(i).outputs),'ALLMETS'); L=~cellfun('isempty',strfind(upper(taskStructure(i).outputs),'ALLMETSIN')); %Check that all metabolites are either real metabolites or @@ -244,8 +258,11 @@ end %Add changed bounds if ~isempty(taskStructure(i).changed) - tModel=setParam(tModel,'lb',taskStructure(i).changed,taskStructure(i).LBrxn); - tModel=setParam(tModel,'ub',taskStructure(i).changed,taskStructure(i).UBrxn); + %NaN bounds mean "leave unchanged" + lbSet=~isnan(taskStructure(i).LBrxn); + ubSet=~isnan(taskStructure(i).UBrxn); + tModel=setParam(tModel,'lb',taskStructure(i).changed(lbSet),taskStructure(i).LBrxn(lbSet)); + tModel=setParam(tModel,'ub',taskStructure(i).changed(ubSet),taskStructure(i).UBrxn(ubSet)); end %Solve and print @@ -304,3 +321,13 @@ taskReport.ok=ok; end + +function dupMets = duplicateMetNames(names,modelMets) +% Return the (upper-case) names that match more than one model metabolite. +dupMets = unique(upper(names)); +isDup = false(size(dupMets)); +for k = 1:numel(dupMets) + isDup(k) = sum(strcmp(modelMets,dupMets{k})) > 1; +end +dupMets = dupMets(isDup); +end diff --git a/src/tasks/parseTaskList.m b/src/tasks/parseTaskList.m index 7a177db9..0c946d51 100755 --- a/src/tasks/parseTaskList.m +++ b/src/tasks/parseTaskList.m @@ -41,14 +41,20 @@ % Several IDs can be delimited by ";". If so, then the same bounds % are used for all reactions. If that is not wanted, then use several % rows for the task -% - CHANGED LB : lower bound for the reaction -% - CHANGED UB : upper bound for the reaction +% - CHANGED LB : lower bound for the reaction. If empty, the lower bound +% of the reaction is left as it is +% - CHANGED UB : upper bound for the reaction. If empty, the upper bound +% of the reaction is left as it is % - SHOULD FAIL : true if the correct behavior of the model is to not -% have a feasible solution given the constraints (optional, default +% have a feasible solution given the constraints. Any non-empty value +% except 0, false and no (any case) counts as true (optional, default % false) % - PRINT FLUX : true if the function should print the corresponding -% flux distribution for a task. Can be useful for testing (optional, -% default false) +% flux distribution for a task. Can be useful for testing. Read in +% the same way as SHOULD FAIL (optional, default false) +% +% Headers are matched ignoring case and surrounding whitespace. In the +% ";"-delimited cells, items are trimmed and empty items are dropped. % % Returns % ------- @@ -71,8 +77,10 @@ % reversible and 0 for irreversible) % - UBequ : array with upper bounds on equations (default, 1000) % - changed : cell array with reactions to change bounds for -% - LBrxn : array with lower bounds on changed reactions -% - UBrxn : array with upper bounds on changed reactions +% - LBrxn : array with lower bounds on changed reactions (NaN if empty, +% meaning the bound is left unchanged) +% - UBrxn : array with upper bounds on changed reactions (NaN if empty, +% meaning the bound is left unchanged) % % Examples % -------- @@ -136,8 +144,12 @@ %Captions columns={'ID';'DESCRIPTION';'IN';'IN LB';'IN UB';'OUT';'OUT LB';'OUT UB';'EQU';'EQU LB';'EQU UB';'CHANGED RXN';'CHANGED LB';'CHANGED UB';'SHOULD FAIL';'PRINT FLUX';'COMMENTS'}; -%Match the columns -[I, colI]=ismember(columns,raw(1,:)); +%Match the columns. Headers are matched ignoring case and surrounding +%whitespace, so "id" and " IN LB " are accepted. +header=raw(1,:); +isText=cellfun(@ischar,header); +header(isText)=upper(strtrim(header(isText))); +[I, colI]=ismember(columns,header); %Only ID is required. Columns that are absent get index 0 from ismember, which %cannot be used to index raw, so point them at a spare empty column instead: @@ -204,10 +216,10 @@ task.id=raw{2,colI(1)}; end task.description=raw{2,colI(2)}; -if ~isempty(raw{2,colI(15)}) +if isTruthy(raw{2,colI(15)}) task.shouldFail=true; end -if ~isempty(raw{2,colI(16)}) +if isTruthy(raw{2,colI(16)}) task.printFluxes=true; end if ~isempty(raw{2,colI(17)}) @@ -217,14 +229,14 @@ for i=2:size(raw,1) %Set the inputs if ischar(raw{i,colI(3)}) - inputs=regexp(raw{i,colI(3)},';','split'); + inputs=splitList(raw{i,colI(3)}); task.inputs=[task.inputs;inputs(:)]; task.LBin=[task.LBin;ones(numel(inputs),1)*raw{i,colI(4)}]; task.UBin=[task.UBin;ones(numel(inputs),1)*raw{i,colI(5)}]; end %Set the outputs if ischar(raw{i,colI(6)}) - outputs=regexp(raw{i,colI(6)},';','split'); + outputs=splitList(raw{i,colI(6)}); task.outputs=[task.outputs;outputs(:)]; task.LBout=[task.LBout;ones(numel(outputs),1)*raw{i,colI(7)}]; task.UBout=[task.UBout;ones(numel(outputs),1)*raw{i,colI(8)}]; @@ -249,10 +261,11 @@ end %Add changed bounds if ischar(raw{i,colI(12)}) - changed=regexp(raw{i,colI(12)},';','split'); + changed=splitList(raw{i,colI(12)}); task.changed=[task.changed;changed(:)]; - task.LBrxn=[task.LBrxn;ones(numel(changed),1)*raw{i,colI(13)}]; - task.UBrxn=[task.UBrxn;ones(numel(changed),1)*raw{i,colI(14)}]; + %An empty bound means "leave this bound as it is", stored as NaN + task.LBrxn=[task.LBrxn;ones(numel(changed),1)*boundOrNaN(raw{i,colI(13)})]; + task.UBrxn=[task.UBrxn;ones(numel(changed),1)*boundOrNaN(raw{i,colI(14)})]; end %Check if it should add more constraints @@ -271,10 +284,10 @@ task.id=raw{i+1,colI(1)}; end task.description=raw{i+1,colI(2)}; - if ~isempty(raw{i+1,colI(15)}) + if isTruthy(raw{i+1,colI(15)}) task.shouldFail=true; end - if ~isempty(raw{i+1,colI(16)}) + if isTruthy(raw{i+1,colI(16)}) task.printFluxes=true; end if ~isempty(raw{i+1,colI(17)}) @@ -286,3 +299,33 @@ %Should add more checks, such as unique IDs and missing headers end + +function tf = isTruthy(x) +% A flag cell is true unless it is empty, 0, or one of the texts "0", "false" +% and "no" (any case). +if isempty(x) + tf = false; +elseif ischar(x) || isstring(x) + tf = ~ismember(lower(strtrim(char(x))),{'','0','false','no'}); +elseif isnumeric(x) || islogical(x) + tf = any(x(:) ~= 0); +else + tf = true; +end +end + +function items = splitList(str) +% Split a ";"-delimited cell on ";", trim each item and drop empty ones, so +% that "A; B;" gives {'A';'B'}. +items = strtrim(regexp(str,';','split')); +items = items(~cellfun(@isempty,items)); +end + +function b = boundOrNaN(x) +% Bound for a changed reaction: NaN (leave unchanged) when the cell is empty. +if isempty(x) + b = NaN; +else + b = x; +end +end diff --git a/tests/function_tests/tINIT.m b/tests/function_tests/tINIT.m index 56214098..c7bd276a 100644 --- a/tests/function_tests/tINIT.m +++ b/tests/function_tests/tINIT.m @@ -62,6 +62,29 @@ function removeLowScoreGenesKeepsFieldsAligned(testCase) testCase.verifyEqual(full(mm.rxnGeneMat), ones(1, numel(mm.genes))); end + function removeLowScoreGenesTieKeepsFirstGeneInRule(testCase) + % With every isozyme negative and tied, the gene listed first in + % the rule is kept (not the alphabetically first one), and the + % order of the surviving genes follows the rule. + m = struct(); + m.id = 'test'; + m.rxns = {'R1';'R2'}; m.rxnNames = {'R1';'R2'}; + m.mets = {'A';'B'}; m.metNames = {'A';'B'}; m.metComps = [1;1]; + m.comps = {'c'}; m.compNames = {'c'}; + m.S = sparse([-1 -1;1 1]); m.lb = [0;0]; m.ub = [1000;1000]; + m.rev = [0;0]; m.c = [0;0]; m.b = [0;0]; + m.genes = {'G1';'G2';'G3';'G4'}; + m.grRules = {'G2 or G1';'G4 or G3 or G2'}; + m.rxnGeneMat = sparse([1 1 0 0;0 1 1 1]); + evalc('[mm, removed] = removeLowScoreGenes(m, [-1;-1;1;-1]);'); + testCase.verifyEqual(mm.grRules{1}, 'G2'); + testCase.verifyEqual(mm.grRules{2}, 'G3'); + testCase.verifyEqual(sort(removed), {'G1';'G4'}); + + evalc('mm = removeLowScoreGenes(m, [1;1;1;1]);'); + testCase.verifyEqual(mm.grRules, m.grRules); + end + function reverseRxnsRuns(testCase) % R1 = '=> a[s]'; R3 = 'a[c] <=> b[c] + c[c]' testModel = getTstModel(); diff --git a/tests/function_tests/tIO.m b/tests/function_tests/tIO.m index 86a06ae6..30f0ff6a 100644 --- a/tests/function_tests/tIO.m +++ b/tests/function_tests/tIO.m @@ -207,6 +207,54 @@ function writeReadYAMLRoundTrip(testCase) testCase.verifyEqual(numel(m2.rxns), numel(testCase.model.rxns)); end + function writeYAMLmodelQuotesYAML11Booleans(testCase) + % A YAML 1.1 loader such as PyYAML reads a plain y, n, yes, no, + % on or off in any case as a boolean, so these values are + % double-quoted. Longer words that start with them stay plain. + m = testCase.model; + quoted = {'NO'; 'Y'; 'n'; 'on'; 'yes'; 'off'}; + plain = {'NOS1'; 'onion'; 'off-target'}; + m.metNames(1:6) = quoted; + m.metFormulas(1:6) = quoted; + m.metNames(7:9) = plain; + f = [tempname '.yml']; + testCase.addTeardown(@() delete(f)); + evalc('writeYAMLmodel(m, f);'); + lines = splitlines(fileread(f)); + for i = 1:numel(quoted) + testCase.verifyTrue(any(strcmp(lines, [' - name: "' quoted{i} '"'])), quoted{i}); + testCase.verifyTrue(any(strcmp(lines, [' - formula: "' quoted{i} '"'])), quoted{i}); + end + for i = 1:numel(plain) + testCase.verifyTrue(any(strcmp(lines, [' - name: ' plain{i}])), plain{i}); + end + evalc('m2 = readYAMLmodel(f);'); + testCase.verifyEqual(m2.metNames(1:9), m.metNames(1:9)); + testCase.verifyEqual(m2.metFormulas(1:6), quoted); + end + + function writeYAMLmodelWritesShortestRoundTripFloats(testCase) + % Non-integers are written with as many digits as needed to read + % back exactly (as Python does), not cut at 15 digits. + m = testCase.model; + m.lb(1) = 0.1 + 0.2; + m.ub(1) = 1/3; + m.lb(2) = 1.5e-7; + m.ub(2) = 2.5; + f = [tempname '.yml']; + testCase.addTeardown(@() delete(f)); + evalc('writeYAMLmodel(m, f);'); + txt = fileread(f); + testCase.verifySubstring(txt, 'upper_bound: 0.3333333333333333'); + testCase.verifySubstring(txt, 'lower_bound: 0.30000000000000004'); + testCase.verifySubstring(txt, 'lower_bound: 1.5e-07'); + testCase.verifySubstring(txt, 'upper_bound: 2.5'); + evalc('m2 = readYAMLmodel(f);'); + testCase.verifyEqual(m2.lb(1), m.lb(1)); + testCase.verifyEqual(m2.ub(1), m.ub(1)); + testCase.verifyEqual(m2.lb(2), m.lb(2)); + end + function readYAMLmodelJoinsContinuationStartingWithDash(testCase) % A writer that wraps long scalars can break a name right % before a dash, so the continuation line starts with "-" diff --git a/tests/function_tests/tManipulation.m b/tests/function_tests/tManipulation.m index 67ec0bcd..77e57dfe 100644 --- a/tests/function_tests/tManipulation.m +++ b/tests/function_tests/tManipulation.m @@ -608,6 +608,39 @@ function replaceMetsByNameKeepsBShapeAndRuns(testCase) testCase.verifyEqual(numel(m2.mets), 2); % oxygen and o2 merged end + function replaceMetsByNameLeavesUnrelatedDuplicatesAlone(testCase) + % Merging "oxygen" into "o2" makes r1 a duplicate of r2, which is + % contracted. unrelated1 and unrelated2 are duplicates of each + % other that do not involve either metabolite; contractModel is + % restricted to the reactions of the merged metabolite, so they + % must stay. (By name the replacement's row is the one deleted, + % so the restriction has to follow the metabolite that is left.) + m = tManipulation.duplicateModel(); + evalc('m2 = replaceMets(m, ''oxygen'', ''o2'');'); + testCase.verifyEqual(numel(m2.mets), 3); + testCase.verifyTrue(all(ismember({'unrelated1';'unrelated2'}, m2.rxns))); + testCase.verifyEqual(sum(ismember({'r1';'r2'}, m2.rxns)), 1); + testCase.verifyEqual(numel(m2.rxns), 3); + end + + function replaceMetsByIdLeavesUnrelatedDuplicatesAlone(testCase) + m = tManipulation.duplicateModel(); + evalc('m2 = replaceMets(m, ''a'', ''b'', ''identifiers'', true);'); + testCase.verifyEqual(numel(m2.mets), 3); + testCase.verifyTrue(all(ismember({'unrelated1';'unrelated2'}, m2.rxns))); + testCase.verifyEqual(sum(ismember({'r1';'r2'}, m2.rxns)), 1); + end + + function replaceMetsByNameStillContractsTheNewDuplicate(testCase) + % The restriction must not switch the contraction off: r1 and + % r2 become identical once the two metabolites are one. + m = tManipulation.duplicateModel(); + evalc('m2 = replaceMets(m, ''oxygen'', ''o2'');'); + testCase.verifyLessThan(numel(m2.rxns), numel(m.rxns)); + rowO2 = find(strcmp(m2.metNames, 'o2')); + testCase.verifyEqual(numel(rowO2), 1); + end + function setExchangeBoundsFindsAllMultiExchangeMets(testCase) % Every metabolite exchanged by more than one reaction must be % reported, not just whichever one happens to line up between @@ -897,6 +930,25 @@ function expandModelLeavesDnfAlone(testCase) methods (Static, Access = private) + function m = duplicateModel() + % a and b are named oxygen and o2; r1 (a => x) and r2 (b => x) become + % duplicates once they are one metabolite. unrelated1 and unrelated2 + % (both d =>) are duplicates of each other from the start. + m = struct(); + m.id = 'dup'; + m.rxns = {'r1';'r2';'unrelated1';'unrelated2'}; + m.rxnNames = m.rxns; + m.mets = {'a';'b';'x';'d'}; + m.metNames = {'oxygen';'o2';'x';'d'}; + m.metComps = [1;1;1;1]; + m.comps = {'c'}; + m.compNames = {'cytosol'}; + m.S = sparse([-1 0 0 0; 0 -1 0 0; 1 1 0 0; 0 0 -1 -1]); + m.lb = -1000*ones(4,1); m.ub = 1000*ones(4,1); m.rev = ones(4,1); + m.c = zeros(4,1); m.b = zeros(4,1); + m.genes = {}; m.grRules = repmat({''},4,1); m.rxnGeneMat = sparse(4,0); + end + function m = twoMetModel() % Single reaction a -> b in one compartment. m = struct(); diff --git a/tests/function_tests/tTasks.m b/tests/function_tests/tTasks.m index 3b858ecf..11a71eb9 100644 --- a/tests/function_tests/tTasks.m +++ b/tests/function_tests/tTasks.m @@ -26,6 +26,80 @@ function parseTaskListReturnsStruct(testCase) testCase.verifyNotEmpty(ts); end + function checkTasksWarnsWhenSeveralMetabolitesShareAName(testCase) + m = testCase.taskTestModel(); + m.metNames{3} = 'a'; % bc is now also called a[s], like as + m.metComps(3) = 1; + task = testCase.taskTestStruct(); + lastwarn(''); + [~, report] = evalc('checkTasks(m, [], true, false, false, task);'); + testCase.verifySubstring(lastwarn, 'Several model metabolites match'); + testCase.verifyTrue(report.ok == 1); + end + + function parseTaskListReadsFlagsListsAndHeadersLeniently(testCase) + % Flags: "false"/"0"/"no" are false, other text is true. Lists + % separated by ";" are trimmed and empty items dropped. Headers + % are matched ignoring case and whitespace. An empty CHANGED + % bound is NaN (leave the bound unchanged). + rows = { + {'' 'id' 'Description' 'Should Fail' 'IN' ' in lb ' 'IN UB' 'OUT' 'OUT LB' 'OUT UB' 'EQU' 'EQU LB' 'EQU UB' 'changed rxn' 'CHANGED LB' 'CHANGED UB' 'PRINT FLUX' 'COMMENTS'} + {'' 'T1' 'first' 'false' 'A[c]; B[c];' '' '5' '' '' '' '' '' '' 'R1; R2' '' '10' '0' ''} + {'' 'T2' 'second' 'yes' ';C[c]' '2' '' '' '' '' '' '' '' 'R3' '-5' '' 'No' ''} + {'' 'T3' 'third' '0' '' '' '' 'D[c]' '' '' '' '' '' '' '' '' 'TRUE' ''} + }; + nCols = 20; + lines = cell(numel(rows),1); + for k = 1:numel(rows) + r = rows{k}; + r(end+1:nCols) = {''}; + lines{k} = strjoin(r, sprintf('\t')); + end + f = [tempname '.txt']; + fid = fopen(f,'w'); + fprintf(fid,'%s\n',lines{:}); + fclose(fid); + cleanup = onCleanup(@() delete(f)); + + ts = parseTaskList(f); + testCase.verifyEqual(numel(ts), 3); + + testCase.verifyFalse(ts(1).shouldFail); + testCase.verifyFalse(ts(1).printFluxes); + testCase.verifyEqual(ts(1).inputs, {'A[c]';'B[c]'}); + testCase.verifyEqual(ts(1).LBin, [0;0]); + testCase.verifyEqual(ts(1).UBin, [5;5]); + testCase.verifyEqual(ts(1).changed, {'R1';'R2'}); + testCase.verifyTrue(all(isnan(ts(1).LBrxn))); + testCase.verifyEqual(ts(1).UBrxn, [10;10]); + + testCase.verifyTrue(ts(2).shouldFail); + testCase.verifyFalse(ts(2).printFluxes); + testCase.verifyEqual(ts(2).inputs, {'C[c]'}); + testCase.verifyEqual(ts(2).LBin, 2); + testCase.verifyEqual(ts(2).LBrxn, -5); + testCase.verifyTrue(isnan(ts(2).UBrxn)); + + testCase.verifyFalse(ts(3).shouldFail); + testCase.verifyTrue(ts(3).printFluxes); + end + + function checkTasksLeavesNaNChangedBoundsAlone(testCase) + m = testCase.taskTestModel(); + task = testCase.taskTestStruct(); + task.changed = {'R2'}; + task.LBrxn = NaN; + task.UBrxn = 1000; + [~, report] = evalc('checkTasks(m, [], true, false, false, task);'); + testCase.verifyTrue(report.ok == 1); + + % R2 can no longer carry flux from a[s] to a[c] with UB 0, and + % its lower bound (-1000) is kept + task.UBrxn = 0; + [~, report] = evalc('checkTasks(m, [], true, false, false, task);'); + testCase.verifyTrue(report.ok == 0); + end + function parseTaskListMatchesExpected(testCase) % The text and Excel task definitions must parse to identical % structs with exactly the expected fields and values.