diff --git a/CHANGELOG.md b/CHANGELOG.md index 41cdf996..f197f7a1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,62 @@ # Changelog +## Unreleased + +### Manipulation + +* `mergeCompartments`: a merged metabolite is the copy already in compartment `s`, if any, + otherwise the copy whose id sorts first, with every field of that copy (was: the first copy + in model order), so `ACCOA_c` and `ACCOA_m` become `ACCOA_c` whatever their order. +* `simplifyModel` with `constrainReversible`: on a model without a reversible reaction this step + changes nothing (was: an error "logical indices ... outside of the array bounds", as + `getAllowedBounds` took the empty reaction list for all reactions). `prepINITModel` therefore + runs on a model whose reactions are all irreversible after the first simplification, with the + result of raven-toolbox's `prep_init_model`. +* `removeGenes`: the help states that, when at least one of the requested genes is in the model, + genes no remaining reaction uses are removed as well; when none is, no gene, reaction or bound + changes. + +### Gap-filling + +* `gapFillTopological`: without `targets` and without an objective reaction it raises + `RAVEN:badInput` (was: an error without identifier). Help states the default seeds and + targets. + +### Context-specific models (ftINIT / tINIT) + +* `ftINIT` with `resolveTies`: each main-extraction step takes its fluxes from the + least-total-flux solution for the tie-broken on/off state, so the direction and magnitude + at which the next step forces a reaction no longer depend on the solver's vertex. + That phase stops at a relative or absolute gap of 1e-6. +* `ftINIT` with `resolveTies`: among the least-total-flux solutions, each main-extraction + step takes its fluxes from the one with the least id-rank-weighted backward flux, so the + direction of a flux that carries the same total either way no longer depends on the + solver. The id-rank phases rank a merged reaction by the smallest original reaction id in + its group, so they no longer depend on which member names it. +* `ftINITInternalAlg` with `resolveTies`: when the primary solve stops at the time limit, the + tie-break phases are skipped and the primary solution is kept, as in raven-toolbox (was: the + phases ran on the unproven incumbent). The count and id-rank phases get half of + `params.TimeLimit` each, rounded down and at least 1 s (was: the full limit). +* `ftINIT` with `resolveTies`: the last step skips the minimal-flux and direction phases, as + no step reads its fluxes (new `ftINITInternalAlg` argument `pinFluxes`). The reaction set + does not change; the fluxes and `metProduction` of the last step are those of the id-rank + phase. +* `ftINITInternalAlg`: new optional positional arguments `rankKeys` (14th, the key each + reaction is sorted by in the id-rank phases, default `model.rxns`) and `pinFluxes` (15th, + default `true`). The "tie resolution did not converge" warning of `resolveTies` also names + `phase2c-minflux` and `phase2d-direction` when the minimal-flux or direction phase does not + prove optimal. +* `mergeLinear`: groups are numbered 1, 2, ... in the order of the first reaction of each + group in `origRxnIds` (was: the order in which the reactions were merged). +* `mergeLinear`: a reaction flipped after the first merge pass is the one marked in + `reversedRxns` (was: the reaction at the same index of the original model). + +### Utilities + +* `grRuleToDNF`: help states the grammar it reads, that of `parseGrRule` (`and`/`or` in any + case, `&`, `&&`, `|`, `||`, operators as whole words, so `G1&G2` is one gene id), and that a + malformed rule raises `RAVEN:badGrRule`. + ## 3.0.0b1 — 2026-09-07 First beta toward 3.0. Changes below are relative to RAVEN 2.9.3, the last 2.x release: diff --git a/src/INIT/ftINIT.m b/src/INIT/ftINIT.m index 1c66ba91..61de0623 100644 --- a/src/INIT/ftINIT.m +++ b/src/INIT/ftINIT.m @@ -71,10 +71,21 @@ % if true, pin each step's degenerate MILP to a canonical answer instead % of relying on the solver Seed alone: hold the score objective at its % optimum, then minimise the count of removable ("on") reactions, then, -% among the sparsest, minimise their summed reaction-id rank. Applies to -% every main-extraction step (ftINITInternalAlg) and to the final -% gap-filling step (fitTasks' resolveTies) -- see raven-gecko-parity#104 -% and raven-toolbox's resolve_ties (default false). +% among the sparsest, minimise their summed reaction-id rank. This +% applies to every main-extraction step and to the final gap-filling +% step (fitTasks' resolveTies) -- see raven-gecko-parity#104 and +% raven-toolbox's resolve_ties. In each main-extraction step but the +% last (ftINITInternalAlg) the fluxes, which set the direction and +% magnitude a reaction is forced at in the next step, then come from the +% least-total-flux solution for that on/off state and, among those, from +% the one with the least id-rank-weighted backward flux (a merged +% reaction ranked by the smallest original id in its group). The last +% step skips these flux phases, as no step reads its fluxes. A step +% whose primary solve stops at the time limit keeps its primary +% solution without any of these phases; the count and id-rank phases +% get half of the step's TimeLimit each, and a warning +% names each phase that does not prove optimal or finds no solution +% (default false). % proveAbsGap : double % if set, prove each step to this fixed *absolute* MILP gap in a single % solve, instead of the relative-gap escalation below (stp.MILPParams). @@ -191,6 +202,18 @@ %essential rxn. fluxes = ones(length(prepData.minModel.rxns), 1).*0.1; +%resolveTies ranks a merged reaction by its group's smallest original id, +%which does not depend on which member names the merged reaction. +rankKeys = prepData.minModel.rxns; +[inRef, refInd] = ismember(rankKeys, prepData.refModel.rxns); +for i = find(inRef)' + g = prepData.groupIds(refInd(i)); + if g ~= 0 + members = sort(prepData.refModel.rxns(prepData.groupIds == g)); + rankKeys{i} = members{1}; + end +end + for initStep = 1:length(INITSteps) disp(['ftINIT: Running step ' num2str(initStep)]) stp = INITSteps{initStep}; @@ -199,6 +222,9 @@ error('RAVEN:badInput', '%s', 'RxnsToIgnoreMask may not cover rxns not covered in previous steps, but the other way around is fine.'); end rxnsToIgnoreLastStep = stp.RxnsToIgnoreMask; + %The fluxes of a step set the forced directions of the next one, so the + %last step skips the flux phases of resolveTies. + pinFluxes = initStep < numel(INITSteps); mm = prepData.minModel; @@ -256,7 +282,7 @@ params.MIPGap = 0; params.MIPGapAbs = proveAbsGap; try - [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, [], fluxes, verbose, 0.1, resolveTies); + [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, [], fluxes, verbose, 0.1, resolveTies, rankKeys, pinFluxes); fluxes1(toRev) = -fluxes1(toRev); success = true; if fullMipRes.stat ~= 1 @@ -310,7 +336,7 @@ if ~isempty(fullMipRes) startVals = fullMipRes.full; end - [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, startVals, fluxes, verbose, 0.1, resolveTies); + [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, startVals, fluxes, verbose, 0.1, resolveTies, rankKeys, pinFluxes); %This is a bit tricky - since we reversed some reactions, those fluxes also need to be reversed fluxes1(toRev) = -fluxes1(toRev); diff --git a/src/INIT/ftINITInternalAlg.m b/src/INIT/ftINITInternalAlg.m index 9e8ebc84..af29c682 100644 --- a/src/INIT/ftINITInternalAlg.m +++ b/src/INIT/ftINITInternalAlg.m @@ -1,4 +1,4 @@ -function [deletedRxns,metProduction,res,turnedOnRxns,fluxes]=ftINITInternalAlg(model,rxnScores,metData,essentialRxns,prodWeight,allowExcretion,remPosRev,params,startVals,fluxes,verbose,forceOn,resolveTies) +function [deletedRxns,metProduction,res,turnedOnRxns,fluxes]=ftINITInternalAlg(model,rxnScores,metData,essentialRxns,prodWeight,allowExcretion,remPosRev,params,startVals,fluxes,verbose,forceOn,resolveTies,rankKeys,pinFluxes) % ftINITInternalAlg % This function runs the MILP for a step in ftINIT. % @@ -44,9 +44,30 @@ % rank. Positive-score reactions carry a continuous, not a % free binary, on/off variable (pinned near 1 by the score % objective), so only the negative-score ones are a genuine -% tie to resolve -- matches raven-toolbox's resolve_ties / -% _resolve_ties (raven-gecko-parity#104) (optional, default -% false) +% tie to resolve. The fluxes are then taken from the +% least-total-flux solution for that on/off state and, +% among those, from the one with the least id-rank-weighted +% backward flux, so the directions and forced magnitudes +% the next ftINIT step reads from them are pinned too -- +% matches raven-toolbox's resolve_ties +% (raven-gecko-parity#104). The phases are skipped when +% the primary solve stops at params.TimeLimit, and the +% primary solution is kept. Phases 2a and 2b get half of +% params.TimeLimit each (rounded down, at least 1 s), the +% flux phases the full limit. A warning names each phase +% that does not prove optimal or finds no solution; when a +% flux phase finds none, the fluxes of the phase before it +% are kept (optional, default false) +% rankKeys cell array with, for each reaction in model.rxns, the key +% that the id-rank phases of resolveTies sort it by +% (optional, default model.rxns). ftINIT passes the smallest +% original reaction id of each merge group, so the order +% does not depend on which member names the merged reaction +% pinFluxes if false, resolveTies skips phases 2c and 2d and the +% fluxes and metProduction are those of phase 2b +% (optional, default true). +% ftINIT passes false in its last step, whose fluxes no +% later step reads % % deletedRxns reactions which were deleted by the algorithm (only % rxns included in the problem) @@ -70,9 +91,15 @@ if nargin < 12 || isempty(forceOn) forceOn = 0.1; end +if nargin < 14 || isempty(rankKeys) + rankKeys = model.rxns; +end if nargin < 13 || isempty(resolveTies) resolveTies = false; end +if nargin < 15 || isempty(pinFluxes) + pinFluxes = true; +end if isempty(essentialRxns) essentialRxns={}; end @@ -453,7 +480,10 @@ error('RAVEN:badInput', '%s', EM); end -if resolveTies && (nNegIrrev + nNegRev) > 0 +if resolveTies && (nNegIrrev + nNegRev) > 0 && ~res.hitTimeLimit + %Skipped when the primary solve stopped at the time limit: the floor + %constraint would pin an unproven incumbent, so the phases cannot select + %among true optima. %This MILP is highly degenerate: many reaction subsets reach the same %score optimum, and the solver returns an arbitrary one -- reproducible %for a fixed solver build + Seed, but fragile to a solver version bump @@ -483,6 +513,10 @@ tieParams = params; tieParams.MIPGap = 0; tieParams.MIPGapAbs = 0.4; + %Phases 2a and 2b get half of the time limit each, as in raven-toolbox. + if isfield(tieParams,'TimeLimit') + tieParams.TimeLimit = max(1,floor(tieParams.TimeLimit*0.5)); + end unproven = {}; @@ -508,7 +542,7 @@ tieProb2b.b = [tieProb.b; kmin+0.5]; tieProb2b.csense = [tieProb.csense 'L']; - [~, sortOrder] = sort(model.rxns(negRxnIdx)); + [~, sortOrder] = sort(rankKeys(negRxnIdx)); ranks = zeros(numel(negRxnIdx),1); ranks(sortOrder) = (1:numel(negRxnIdx))'; idObj = zeros(size(prob.c)); @@ -522,6 +556,106 @@ if isFeasible2b %Adopt the tie-resolved solution for the on/off read-out below. res = res2b; + + %Phases 2c and 2d change the fluxes, not the on/off state; the next + %ftINIT step reads the fluxes, and pinFluxes is false when no step + %follows. + if pinFluxes + %Phase 2c: the least total absolute flux for this on/off state. + %The next step forces each turned-on reaction in the direction + %of its flux, at min(0.99*|flux|,0.1), so the flux read-out + %comes from a minimal-flux solution instead of an arbitrary + %vertex. A binary on/off variable is fixed at its value; a + %continuous one that is on keeps at least its value, one that + %is off at most its value. Directions of reversible scored + %reactions stay free, so this is still a MILP. + %It stops at a relative or absolute gap of 1e-6, as raven-toolbox's + %phase 2c. + nCols = numel(prob.c); + minProb = tieProb2b; + if isfield(minProb,'start') + minProb = rmfield(minProb,'start'); + end + x = res2b.full; + negVal = round(x(negCols)); + minProb.lb(negCols) = negVal; + minProb.ub(negCols) = negVal; + posCols = [onoffPosIrrev, onoffPosRev]; + posVal = min(max(x(posCols),0),1); + posOn = posVal >= 0.5; + minProb.lb(posCols(posOn)) = posVal(posOn); + minProb.ub(posCols(~posOn)) = posVal(~posOn); + %|v| is v for a forward-only flux, -v for a backward-only one, + %and an added t >= v, t >= -v otherwise. + fwdOnly = minProb.lb(1:nRxns) >= 0; + bwdOnly = ~fwdOnly & minProb.ub(1:nRxns) <= 0; + tRxns = find(~fwdOnly & ~bwdOnly); + nT = numel(tRxns); + selV = sparse(1:nT, tRxns, 1, nT, nCols); + minProb.A = [minProb.A sparse(size(minProb.A,1),nT); -selV speye(nT); selV speye(nT)]; + minProb.a = minProb.A; + minProb.b = [minProb.b; zeros(2*nT,1)]; + minProb.csense = [minProb.csense repmat('G',1,2*nT)]; + minProb.lb = [minProb.lb; zeros(nT,1)]; + minProb.ub = [minProb.ub; inf(nT,1)]; + minProb.vartype = [minProb.vartype repmat('C',1,nT)]; + minProb.c = [zeros(nCols,1); ones(nT,1)]; + minProb.c(fwdOnly) = 1; + minProb.c(bwdOnly) = -1; + minParams = params; + minParams.MIPGap = 1e-6; + minParams.MIPGapAbs = 1e-6; + res2c = optimizeProb(minProb,minParams,verbose); + [isFeasible2c, isOptimal2c] = checkSolution(res2c); + if ~isOptimal2c + unproven{end+1} = 'phase2c-minflux'; + end + if isFeasible2c + %Only the fluxes change; the status stays that of phase 2b. + res.full = res2c.full(1:nCols); + + %Phase 2d: the least total flux can leave a direction open + %when two directions carry the same total flux, for example + %a pair of reversible reactions closing a loop either way. + %Keep the on/off state of phase 2c, hold the total absolute + %flux at most at its phase-2c value plus + %max(1e-7*total, 1e-6), and minimise the backward flux + %max(0,-v) of the reactions whose flux can take either sign, + %weighted by id rank from the end: sorted by rankKeys, the + %first of n reactions gets weight n and the last 1, so a + %reaction with a lower id prefers its forward direction. + if nT > 0 + total = minProb.c' * res2c.full; + dirProb = minProb; + nCols2c = numel(minProb.c); + [~, ord] = sort(rankKeys(tRxns)); + w = zeros(nT,1); + w(ord) = (nT:-1:1)'; + dirProb.A = [minProb.A sparse(size(minProb.A,1),nT); ... + minProb.c' sparse(1,nT); ... + selV sparse(nT,nT) speye(nT)]; + dirProb.a = dirProb.A; + dirProb.b = [minProb.b; total + max(abs(total)*1e-7, 1e-6); zeros(nT,1)]; + dirProb.csense = [minProb.csense 'L' repmat('G',1,nT)]; + dirProb.lb = [minProb.lb; zeros(nT,1)]; + dirProb.ub = [minProb.ub; inf(nT,1)]; + dirProb.vartype = [minProb.vartype repmat('C',1,nT)]; + dirProb.c = [zeros(nCols2c,1); w]; + %The phase-2c solution is the MIP start: the tight cap on + %the total flux leaves the solver's heuristics little room + %to find a first solution. + dirProb.start = [res2c.full; max(0, -res2c.full(tRxns))]; + res2d = optimizeProb(dirProb,minParams,verbose); + [isFeasible2d, isOptimal2d] = checkSolution(res2d); + if ~isOptimal2d + unproven{end+1} = 'phase2d-direction'; + end + if isFeasible2d + res.full = res2d.full(1:nCols); + end + end + end + end end end if ~isempty(unproven) diff --git a/src/INIT/mergeLinear.m b/src/INIT/mergeLinear.m index 068fab65..66847e98 100644 --- a/src/INIT/mergeLinear.m +++ b/src/INIT/mergeLinear.m @@ -13,7 +13,8 @@ % reducedModel an updated model structure % origRxnIds Cell array: original rxn ids, used together with the groupIds variable % groupIds Vector of ids: merged rxns have the same id, same length as origRxnIds -% 0 means a reaction was not merged. +% 0 means a reaction was not merged. The groups are numbered 1, 2, ... +% in the order of their first reaction in origRxnIds. % reversedRxns Vector of booleans, saying for each rxn if it is reversed in the combined rxns. % % Usage: [reducedModel,origRxnIds,groupIds]=mergeLinear(model,noMergeRxns) @@ -196,7 +197,12 @@ notInUse=sum(reducedModel.S~=0,2)==0; reducedModel=removeMets(reducedModel,notInUse); end - + +%Number the groups in the order of their first reaction in origRxnIds, so +%the numbers do not depend on the order in which the reactions were merged +merged=groupIds~=0; +[~,~,groupIds(merged)]=unique(groupIds(merged),'stable'); + function [model1,reversedRxns1] = flipRxn(model1, rxnInd, groupIds1, reversedRxns1) model1.S(:,rxnInd) = model1.S(:,rxnInd)*-1; %swap the bounds @@ -208,10 +214,13 @@ %now take care of the reversedRxns - if this is a group, reverse all of the %reactions in the group in the reversedRxns index - they will all be reversed at the - %same time since they are the same rxn. - rxnIndices = rxnInd; - if groupIds1(rxnInd) > 0 - rxnIndices = find(groupIds1 == groupIds1(rxnInd)); + %same time since they are the same rxn. rxnInd indexes model1, which loses + %reactions after each pass, while groupIds1 and reversedRxns1 follow + %origRxnIds, so the reaction is looked up by its id. + origInd = origRxnIdMap(model1.rxns{rxnInd}); + rxnIndices = origInd; + if groupIds1(origInd) > 0 + rxnIndices = find(groupIds1 == groupIds1(origInd)); end reversedRxns1(rxnIndices) = ~reversedRxns1(rxnIndices); end diff --git a/src/gapfilling/gapFillTopological.m b/src/gapfilling/gapFillTopological.m index fe6e335e..bf832c02 100644 --- a/src/gapfilling/gapFillTopological.m +++ b/src/gapfilling/gapFillTopological.m @@ -34,8 +34,11 @@ % metabolites of every exchange reaction whose bounds let it supply the % metabolite, i.e. getExchangeRxns 'uptake' plus 'reverse'. % 'targets' (default []) -% Metabolite IDs (cell array) that should be produced. Default: -% substrates of the objective (biomass) reaction. +% Metabolite IDs (cell array) that should be produced. Default: every +% metabolite of the objective (biomass) reaction, substrates and +% products; the first one when there are several objective reactions. +% Without targets and without an objective reaction the function raises +% RAVEN:badInput. % 'verbose' (default true) % print a summary of the results. % @@ -99,10 +102,10 @@ % ---- Identify target metabolites ---- if isempty(targets) - % Default: substrates of the objective (biomass) reaction + % Default: every metabolite of the objective (biomass) reaction objRxns = find(model.c ~= 0); if isempty(objRxns) - error(['gapFillTopological: no objective function defined. ' ... + error('RAVEN:badInput', ['gapFillTopological: no objective function defined. ' ... 'Provide targets manually via the ''targets'' option.']); end biomassIdx = objRxns(1); diff --git a/src/manipulation/mergeCompartments.m b/src/manipulation/mergeCompartments.m index 80c7982e..ce4e80e9 100755 --- a/src/manipulation/mergeCompartments.m +++ b/src/manipulation/mergeCompartments.m @@ -42,8 +42,13 @@ % be useful for example to ensure that there are metabolic capabilities to % synthesize all metabolites. % -% If the metabolite IDs reflect the compartment that they are in the IDs may -% no longer be representative. +% Metabolites with the same name are merged into one. The merged metabolite +% keeps the id and the other fields of one of the copies: the copy that is +% already in compartment "s" if there is one, otherwise the copy whose id +% comes first in character code order (for example ACCOA_c before ACCOA_m). +% With keepUnconstrained only the copies that are not unconstrained are +% candidates. If the metabolite IDs reflect the compartment that they are in +% the IDs may no longer be representative. p=parseRAVENargs(varargin, {'keepUnconstrained',false; 'deleteRxnsWithOneMet',false; 'distReverse',true}); keepUnconstrained=p.keepUnconstrained; @@ -66,9 +71,8 @@ end end -%Loop through each metabolite, and if it is not unconstrained then change -%the S matrix to use the metabolite with the lowest index in model.comps -%instead +%Loop through each metabolite name, and if a copy is not unconstrained then +%change the S matrix to use the copy chosen as described in the Notes %Group metabolite indices by name once, rather than re-scanning all of %model.metNames (via ismember) for every unique name, which is quadratic %in a genome-scale model's (unique names x mets) count. @@ -79,19 +83,26 @@ I=false(numel(model.metNames),1); I(sameNameIdx{i})=true; - %Find the first of those that is not unconstrained. This is the one - %that the other "un-unconstrained" should be changed to. + %Choose the copy that the other "un-unconstrained" copies are merged + %into: one already in compartment "s" if there is one, otherwise the + %one whose id sorts first. if keepUnconstrained==true - mergeTo=find(I & model.unconstrained==false,1); + candidates=find(I & model.unconstrained==false); %This could happen if there is only one metabolite and it is %unconstrained - if isempty(mergeTo) + if isempty(candidates) continue; end else - mergeTo=find(I,1); + candidates=find(I); end + inMerged=candidates(strcmp(model.comps(model.metComps(candidates)),'s')); + if ~isempty(inMerged) + candidates=inMerged; + end + [~,first]=sort(model.mets(candidates)); + mergeTo=candidates(first(1)); I(mergeTo)=false; %Do not do anything for itself I=find(I); diff --git a/src/manipulation/removeGenes.m b/src/manipulation/removeGenes.m index e2a50f32..00f3e30f 100755 --- a/src/manipulation/removeGenes.m +++ b/src/manipulation/removeGenes.m @@ -23,7 +23,15 @@ % Returns % ------- % reducedModel : struct -% an updated model structure. +% an updated model structure. Genes in genesToRemove that are not in the +% model are ignored. When none of them is in the model, no gene, reaction +% or bound changes (with standardizeRules true the gene rules are still +% written in the standard format). Otherwise, besides genesToRemove, every +% gene that no remaining reaction used before the removal is deleted as +% well, also when it was not requested (for example a gene that no +% reaction uses). A gene whose only reaction is kept but blocked +% (removeBlockedRxns false), such as the other subunit of a complex that +% lost one, stays. % % Examples % -------- diff --git a/src/manipulation/simplifyModel.m b/src/manipulation/simplifyModel.m index e5dd214f..d579230d 100755 --- a/src/manipulation/simplifyModel.m +++ b/src/manipulation/simplifyModel.m @@ -185,30 +185,33 @@ if constrainReversible==true revs=find(reducedModel.rev); - [I,J]=getAllowedBounds(reducedModel,revs); + %getAllowedBounds treats an empty reaction list as all reactions + if ~isempty(revs) + [I,J]=getAllowedBounds(reducedModel,revs); - I=abs(I); - J=abs(J); + I=abs(I); + J=abs(J); - %Get the "small" values - K=I<10^-10; - L=J<10^-10; + %Get the "small" values + K=I<10^-10; + L=J<10^-10; - %Keep the values where only one direction is zero (small) - I(K==L)=[]; - J(K==L)=[]; - revs(K==L)=[]; + %Keep the values where only one direction is zero (small) + I(K==L)=[]; + J(K==L)=[]; + revs(K==L)=[]; - %Change the reversibility of the remaining reactions - reducedModel.rev(revs(J>10^-10))=0; %Ignore reverse direction - reducedModel.lb(revs(J>10^-10))=0; + %Change the reversibility of the remaining reactions + reducedModel.rev(revs(J>10^-10))=0; %Ignore reverse direction + reducedModel.lb(revs(J>10^-10))=0; - toSwitch=revs(I>10^-10); - reducedModel.rev(toSwitch)=0; %Change directionality - reducedModel.ub(toSwitch)=reducedModel.lb(toSwitch)*-1; - reducedModel.lb(toSwitch)=0; - reducedModel.S(:,toSwitch)=reducedModel.S(:,toSwitch).*-1; - reducedModel.c(toSwitch)=reducedModel.c(toSwitch)*-1; + toSwitch=revs(I>10^-10); + reducedModel.rev(toSwitch)=0; %Change directionality + reducedModel.ub(toSwitch)=reducedModel.lb(toSwitch)*-1; + reducedModel.lb(toSwitch)=0; + reducedModel.S(:,toSwitch)=reducedModel.S(:,toSwitch).*-1; + reducedModel.c(toSwitch)=reducedModel.c(toSwitch)*-1; + end end if groupLinear==true if suppressWarnings==false diff --git a/src/utils/grRuleToDNF.m b/src/utils/grRuleToDNF.m index 34da5a64..f2ea4465 100644 --- a/src/utils/grRuleToDNF.m +++ b/src/utils/grRuleToDNF.m @@ -8,7 +8,10 @@ % Parameters % ---------- % rule : char or string or struct -% a grRule, or a tree from parseGrRule. +% a grRule, or a tree from parseGrRule. A grRule is read by parseGrRule: +% genes joined by "and"/"or" in any case or by "&", "&&", "|", "||", +% each operator separated from the gene ids by whitespace or a bracket, +% with AND binding tighter than OR. % % Returns % ------- @@ -18,6 +21,10 @@ % % Notes % ----- +% A malformed rule (an operator without an operand, unbalanced or empty +% parentheses, two genes without an operator between them) raises +% RAVEN:badGrRule. +% % Gene order within a clause follows the order the genes appear in the rule. % Duplicates are not removed: "G1 and G1" stays a two-element clause, since % collapsing it would silently rewrite the user's rule. diff --git a/tests/function_tests/tGapfilling.m b/tests/function_tests/tGapfilling.m index acd68455..b4492fa6 100644 --- a/tests/function_tests/tGapfilling.m +++ b/tests/function_tests/tGapfilling.m @@ -350,6 +350,21 @@ function gapFillTopologicalBlocksWithoutExchange(testCase) testCase.verifyNumElements(result.blockedMets, 3); end + function gapFillTopologicalDefaultTargetsAreObjectiveMets(testCase) + % Objective R2 (B -> C): the default targets are its substrate B + % and its product C. Without an exchange both are blocked. + m = testCase.chainModel(); + evalc('result = gapFillTopological(m, m, ''verbose'', false);'); + testCase.verifyEqual(sort(result.blockedMets(:)), {'B';'C'}); + end + + function gapFillTopologicalDefaultTargetsNeedObjective(testCase) + m = testCase.chainModel(); + m.c = [0;0]; + testCase.verifyError(@() gapFillTopological(m, m, 'seeds', {'A'}, ... + 'verbose', false), 'RAVEN:badInput'); + end + function gapFillMILPRepairsGrowth(testCase) testCase.assumeMILPSolver(); modelDB = testCase.model; modelDB.id = 'DB'; diff --git a/tests/function_tests/tINIT.m b/tests/function_tests/tINIT.m index 4564e29d..dccb1630 100644 --- a/tests/function_tests/tINIT.m +++ b/tests/function_tests/tINIT.m @@ -179,6 +179,75 @@ function mergeLinearGroupsAndScores(testCase) testCase.verifyTrue(all(find(reversedRxns) == [6;9])); end + function mergeLinearNumbersGroupsByFirstRxn(testCase) + % Group numbers follow the first reaction of each group in + % model.rxns, not the merge order. The metabolites are visited + % in model.mets order, so c (R5, R6) merges first. + model = struct(); + model.id = 'numbering'; + model.rxns = {'R1';'R2';'R3';'R4';'R5';'R6'}; + model.rxnNames = model.rxns; + model.mets = {'z';'c';'w';'x';'a';'m';'b';'y'}; + model.metNames = model.mets; + model.comps = {'c'}; + model.compNames = {'c'}; + model.metComps = ones(8,1); + subs = {'x';'a';'m';'b';'z';'c'}; + prods = {'a';'m';'b';'y';'c';'w'}; + model.S = sparse(8,6); + for i = 1:6 + model.S(strcmp(model.mets,subs{i}),i) = -1; + model.S(strcmp(model.mets,prods{i}),i) = 1; + end + model.lb = zeros(6,1); + model.ub = 1000*ones(6,1); + model.rev = zeros(6,1); + model.c = zeros(6,1); + model.b = zeros(8,1); + model.genes = cell(0,1); + model.grRules = repmat({''},6,1); + model.rxnGeneMat = sparse(6,0); + evalc('[reducedModel,~,groupIds] = mergeLinear(model, {});'); + testCase.verifyEqual(groupIds(:), [1;1;1;1;2;2]); + testCase.verifyEqual(sort(reducedModel.rxns), {'R1';'R5'}); + end + + function mergeLinearFlipsRxnInLaterPass(testCase) + % R1 and R2 merge through A in the first pass and R2 is removed, + % so R4 sits at index 3 of the reduced model. The second pass + % merges R1 and R4 through Y and flips R4, which must be the + % reaction marked as reversed (not R3, the third original one). + model = struct(); + model.id = 'lateflip'; + model.rxns = {'R1';'R2';'R3';'R4';'R5';'R6'}; + model.rxnNames = model.rxns; + model.mets = {'A';'Y';'Q';'W'}; + model.metNames = model.mets; + model.comps = {'c'}; + model.compNames = {'c'}; + model.metComps = ones(4,1); + S = zeros(4,6); + S([1 2],1) = [1;1]; % R1: -> A + Y + S([1 2],2) = [-1;1]; % R2: A -> Y + S(3,3) = -1; % R3: Q -> + S([4 2],4) = [-1;1]; % R4: W <=> Y + S(4,5) = 1; % R5: -> W + S(4,6) = -1; % R6: W -> + model.S = sparse(S); + model.lb = [0;0;0;-1000;0;0]; + model.ub = 1000*ones(6,1); + model.rev = double(model.lb < 0); + model.c = zeros(6,1); + model.b = zeros(4,1); + model.genes = cell(0,1); + model.grRules = repmat({''},6,1); + model.rxnGeneMat = sparse(6,0); + evalc('[reducedModel,origRxnIds,groupIds,reversedRxns] = mergeLinear(model, {});'); + testCase.verifyEqual(groupIds(:), [1;1;0;1;0;0]); + testCase.verifyEqual(origRxnIds(reversedRxns), {'R4'}); + testCase.verifyEqual(sort(reducedModel.rxns), {'R1';'R3';'R5';'R6'}); + end + function scoreModelRuns(testCase) arrayData.genes = testCase.model.genes; arrayData.tissues = {'t1'}; @@ -384,6 +453,165 @@ function ftINITResolveTiesAndProveAbsGapStillRunCorrectly(testCase) testCase.verifyTrue(all(strcmp(resBoth.rxns, expected))); end + function ftINITResolveTiesRanksMergedRxnBySmallestOriginalId(testCase) + % a[c] reaches m[c] by path P (Rm, Rc) or path Q (Rn, Rb), each + % path one linear merge group with equal negative scores. Model + % order and group numbers put P first, and the merged reactions + % are named by their first members Rm and Rn, but the smallest + % original id of Q (Rb) sorts before that of P (Rc). resolveTies + % keeps Q; ranking by the merged id, the model order or the group + % number would keep P. + testCase.assumeMILPSolver(); + testModel = getTstModelTwoPaths(); + evalc('prepData = prepINITModel(testModel, getTstModelTasks(), {}, false, {}, ''s'');'); + g = @(r) prepData.groupIds(strcmp(prepData.refModel.rxns, r)); + testCase.verifyEqual(g('Rm'), g('Rc')); + testCase.verifyEqual(g('Rn'), g('Rb')); + testCase.verifyLessThan(g('Rm'), g('Rn')); + testCase.verifyTrue(all(ismember({'Rm';'Rn'}, prepData.minModel.rxns))); + + arrayData.genes = testModel.genes; + arrayData.tissues = {'a'}; + arrayData.levels = getExprForRxnScore([-1;-1;-1;-1;5;5]); + arrayData.threshold = 1; + evalc(['resModel = ftINIT(prepData,arrayData.tissues{1},[],[],' ... + 'arrayData,[],getINITSteps(),true,true,struct(),false,' ... + '''resolveTies'',true);']); + testCase.verifyTrue(all(ismember({'Rn';'Rb'}, resModel.rxns))); + testCase.verifyFalse(any(ismember({'Rm';'Rc'}, resModel.rxns))); + end + + function ftINITPinsFluxesInEveryStepButTheLast(testCase) + % ftINIT passes pinFluxes true to ftINITInternalAlg in each step + % whose fluxes a later step reads and false in the last step. A + % wrapper of ftINITInternalAlg earlier on the path records the + % argument. + testCase.assumeMILPSolver(); + global TINIT_INTERNALALG TINIT_PINFLUXES + TINIT_INTERNALALG = @ftINITInternalAlg; + TINIT_PINFLUXES = false(1,0); + d = tempname; mkdir(d); + fid = fopen(fullfile(d,'ftINITInternalAlg.m'),'w'); + fprintf(fid,'%s\n', ... + 'function varargout = ftINITInternalAlg(varargin)', ... + 'global TINIT_INTERNALALG TINIT_PINFLUXES', ... + 'TINIT_PINFLUXES(end+1) = varargin{15};', ... + '[varargout{1:nargout}] = TINIT_INTERNALALG(varargin{:});', ... + 'end'); + fclose(fid); + testCase.addTeardown(@() rmdir(d,'s')); + addpath(d,'-begin'); + testCase.addTeardown(@() rmpath(d)); + testCase.addTeardown(@() clear('global','TINIT_INTERNALALG','TINIT_PINFLUXES')); + + testModel = getTstModelTwoPaths(); + evalc('prepData = prepINITModel(testModel, getTstModelTasks(), {}, false, {}, ''s'');'); + arrayData.genes = testModel.genes; + arrayData.tissues = {'a'}; + arrayData.levels = getExprForRxnScore([-1;-1;-1;-1;5;5]); + arrayData.threshold = 1; + steps = getINITSteps(); + evalc(['resModel = ftINIT(prepData,arrayData.tissues{1},[],[],' ... + 'arrayData,[],steps,true,true,struct(),false,' ... + '''resolveTies'',true);']); + testCase.verifyEqual(TINIT_PINFLUXES, [true(1,numel(steps)-1) false]); + testCase.verifyTrue(all(ismember({'Rn';'Rb'}, resModel.rxns))); + end + + function ftINITInternalAlgResolveTiesGivesMinimalFluxes(testCase) + % With resolveTies the fluxes come from the least-total-flux + % solution for the tie-broken on/off state. The positive + % reversible R3 (m <=> q) can be on in either direction: backward + % (import q, feed the essential E) carries 0.4 in total, forward + % (a via the free R1 into R3 and E) 0.8. The free loop L1/L2 + % carries nothing and the negative-score dead end R2 stays off. + testCase.assumeMILPSolver(); + m.rxns = {'R1';'R2';'R3';'E';'L1';'L2';'EX_a';'EX_p';'EX_q'}; + m.mets = {'a';'m';'p';'q';'n';'z'}; + % R1 R2 R3 E L1 L2 EXa EXp EXq + m.S = sparse([-1 -1 0 0 0 0 -1 0 0; % a + 1 0 -1 -1 -1 1 0 0 0; % m + 0 0 0 1 0 0 0 -1 0; % p + 0 0 1 0 0 0 0 0 -1; % q + 0 0 0 0 1 -1 0 0 0; % n + 0 1 0 0 0 0 0 0 0]); % z + m.lb = [0;0;-1000;0;0;0;-1000;-1000;-1000]; + m.ub = 1000*ones(9,1); + m.rev = double(m.lb < 0); + scores = [0;-1;1;0;0;0;0;0;0]; + fluxIn = 0.1*ones(9,1); + [deleted,~,~,turnedOn,fluxes] = ftINITInternalAlg(m,scores,[],{'E'}, ... + [],false,false,struct(),[],fluxIn,false,0.1,true); + testCase.verifyTrue(turnedOn(3)); + testCase.verifyTrue(deleted(2)); + % E is forced at min(0.99*0.1,0.1) and needs the 0.1 of R3. + expected = [0;0;-0.1;0.1;0;0;0;0.1;-0.1]; + testCase.verifyEqual(fluxes, expected, 'AbsTol', 1e-6); + end + + function ftINITInternalAlgResolveTiesFixesEqualFluxDirections(testCase) + % The essential E consumes y, which the reversible R1 (y <=> x) + % or R2 (y <=> z) can supply when run backward from the imports + % IN_x or IN_z. Both routes carry the same total flux; with + % resolveTies the higher-ranked reaction carries the backward + % flux: R2 when ranked by id, R1 when rankKeys rank it last. N is + % a negative-score dead end. + testCase.assumeMILPSolver(); + [m, scores] = getTstModelTwoRoutes(); + fluxIn = ones(7,1); + [~,~,~,~,fluxes] = ftINITInternalAlg(m,scores,[],{'E'}, ... + [],false,false,struct(),[],fluxIn,false,0.1,true); + testCase.verifyEqual(fluxes(3:4), [0;-0.1], 'AbsTol', 1e-6); + keys = {'IN_x';'IN_z';'R9';'R2';'E';'OUT_p';'N'}; + [~,~,~,~,fluxes] = ftINITInternalAlg(m,scores,[],{'E'}, ... + [],false,false,struct(),[],fluxIn,false,0.1,true,keys); + testCase.verifyEqual(fluxes(3:4), [-0.1;0], 'AbsTol', 1e-6); + end + + function ftINITInternalAlgResolveTiesSkippedAfterTimeLimit(testCase) + % When the primary solve stops at the time limit, resolveTies + % keeps the primary solution and runs no tie-break phase. A + % wrapper of optimizeProb counts the solves and reports each as + % stopped at the time limit. + testCase.assumeMILPSolver(); + installOptimizeProbCounter(testCase, true); + global TINIT_OPTIMIZEPROB_CALLS + m = getTstModelTwoSources(); + [~,~,res,turnedOn] = ftINITInternalAlg(m,[-1;-1;0;0;0],[],{'E'}, ... + [],false,false,struct(),[],0.1*ones(5,1),false,0.1,true); + testCase.verifyEqual(TINIT_OPTIMIZEPROB_CALLS, 1); + testCase.verifyTrue(res.hitTimeLimit); + testCase.verifyEqual(nnz(turnedOn(1:2)), 1); + end + + function ftINITInternalAlgResolveTiesWithoutPinFluxesSkipsFluxPhases(testCase) + % With pinFluxes false, resolveTies runs the primary solve and + % phases 2a and 2b, not 2c and 2d, and still keeps the lower-id R1. + testCase.assumeMILPSolver(); + installOptimizeProbCounter(testCase, false); + global TINIT_OPTIMIZEPROB_CALLS + m = getTstModelTwoSources(); + [~,~,~,turnedOn] = ftINITInternalAlg(m,[-1;-1;0;0;0],[],{'E'}, ... + [],false,false,struct(),[],0.1*ones(5,1),false,0.1,true,[],false); + testCase.verifyEqual(TINIT_OPTIMIZEPROB_CALLS, 3); + testCase.verifyEqual(turnedOn(1:2), [true false]); + TINIT_OPTIMIZEPROB_CALLS = 0; + ftINITInternalAlg(m,[-1;-1;0;0;0],[],{'E'}, ... + [],false,false,struct(),[],0.1*ones(5,1),false,0.1,true); + testCase.verifyEqual(TINIT_OPTIMIZEPROB_CALLS, 5); + end + + function prepINITModelHandlesModelWithoutReversibleRxns(testCase) + % Every reaction is irreversible, so the second simplification + % has no reversible reaction to constrain. The merge gives the + % same reactions and groups as raven-toolbox's prep_init_model. + testModel = getTstModelIrreversible(); + evalc('prepData = prepINITModel(testModel, [], {}, false, {}, ''s'');'); + testCase.verifyEqual(prepData.minModel.rxns, {'EX_a';'R1';'R2';'T2'}); + testCase.verifyEqual(prepData.groupIds, [1;1;0;2;2;3;3]); + testCase.verifyEqual(prepData.refModel.rxns, testModel.rxns); + end + function ftINITWithTaskRuns(testCase) % A task requiring e[s] from a[s] forces R2 and R7 to be essential. testCase.assumeMILPSolver(); @@ -816,6 +1044,133 @@ function scoreModelRejectsUnknownDataPrecedence(testCase) testModel.b = repmat(0,8,1); end +function installOptimizeProbCounter(testCase, fakeTimeLimit) +%Put a wrapper of optimizeProb first on the path for the rest of the test. It +%counts the solves in the global TINIT_OPTIMIZEPROB_CALLS and, if +%fakeTimeLimit, reports each solve as stopped at the time limit. +global TINIT_OPTIMIZEPROB TINIT_OPTIMIZEPROB_CALLS TINIT_OPTIMIZEPROB_FAKETL +TINIT_OPTIMIZEPROB = @optimizeProb; +TINIT_OPTIMIZEPROB_CALLS = 0; +TINIT_OPTIMIZEPROB_FAKETL = fakeTimeLimit; +d = tempname; mkdir(d); +fid = fopen(fullfile(d,'optimizeProb.m'),'w'); +fprintf(fid,'%s\n', ... + 'function res = optimizeProb(varargin)', ... + 'global TINIT_OPTIMIZEPROB TINIT_OPTIMIZEPROB_CALLS TINIT_OPTIMIZEPROB_FAKETL', ... + 'res = TINIT_OPTIMIZEPROB(varargin{:});', ... + 'TINIT_OPTIMIZEPROB_CALLS = TINIT_OPTIMIZEPROB_CALLS + 1;', ... + 'if TINIT_OPTIMIZEPROB_FAKETL, res.hitTimeLimit = true; end', ... + 'end'); +fclose(fid); +testCase.addTeardown(@() rmdir(d,'s')); +addpath(d,'-begin'); +testCase.addTeardown(@() rmpath(d)); +testCase.addTeardown(@() clear('global','TINIT_OPTIMIZEPROB','TINIT_OPTIMIZEPROB_CALLS','TINIT_OPTIMIZEPROB_FAKETL')); +end + +function m = getTstModelTwoSources() +%R1 and R2 (a => m) can each feed the essential E (m => p): a tie for +%resolveTies. +m.rxns = {'R1';'R2';'E';'EX_a';'EX_p'}; +m.mets = {'a';'m';'p'}; +m.S = sparse([-1 -1 0 -1 0; + 1 1 -1 0 0; + 0 0 1 0 -1]); +m.lb = [0;0;0;-1000;-1000]; +m.ub = 1000*ones(5,1); +m.rev = double(m.lb < 0); +end + +function [m, scores] = getTstModelTwoRoutes() +%The essential E consumes y, which the reversible R1 (y <=> x) or R2 +%(y <=> z) can supply when run backward from the imports IN_x or IN_z. N is a +%negative-score dead end. +m.rxns = {'IN_x';'IN_z';'R1';'R2';'E';'OUT_p';'N'}; +m.mets = {'x';'y';'z';'p';'d'}; +% INx INz R1 R2 E OUTp N +m.S = sparse([1 0 1 0 0 0 -1; % x + 0 0 -1 -1 -1 0 0; % y + 0 1 0 1 0 0 0; % z + 0 0 0 0 1 -1 0; % p + 0 0 0 0 0 0 1]); % d +m.lb = [0;0;-1000;-1000;0;0;0]; +m.ub = 1000*ones(7,1); +m.rev = double(m.lb < 0); +scores = [0;0;0;0;0;0;-1]; +end + +function testModel = getTstModelIrreversible() +%a[s] is imported and reaches b[c] by R1 or by R2, R3; b[c] is exported. +%Every reaction is irreversible. +testModel = struct(); +testModel.id = 'irreversible'; +testModel.rxns = {}; +testModel.S = []; +testModel.rev = []; +testModel.mets = {'as';'ac';'bc';'cc';'bs'}; +testModel.metNames = {'a';'a';'b';'c';'b'}; +testModel.comps = {'s';'c'}; +testModel.compNames = testModel.comps; +testModel.metComps = [1;2;2;2;1]; +testModel.genes = {'G1';'G2';'G3'}; +testModel.grRules = {}; +testModel.rxnGeneMat = []; + +rxnsToAdd = struct(); +rxnsToAdd.rxns = {'EX_a';'T1';'R1';'R2';'R3';'T2';'EX_b'}; +rxnsToAdd.equations = {'=> a[s]';... + 'a[s] => a[c]';... + 'a[c] => b[c]';... + 'a[c] => c[c]';... + 'c[c] => b[c]';... + 'b[c] => b[s]';... + 'b[s] =>'}; +rxnsToAdd.grRules = {'';'';'G1';'G2';'G3';'';''}; +evalc('testModel = addRxns(testModel,rxnsToAdd, 3);'); +testModel.c = [0;0;0;0;0;0;1]; +testModel.ub = repmat(1000,7,1); +testModel.lb = zeros(7,1); +testModel.rxnNames = testModel.rxns; +testModel.b = zeros(5,1); +end + +function testModel = getTstModelTwoPaths() +%a[c] reaches m[c] by Rm, Rc or by Rn, Rb; E, R7 and R8 export e[s]. The +%paths are reversible. +testModel = struct(); +testModel.id = 'twoPaths'; +testModel.rxns = {}; +testModel.S = []; +testModel.rev = []; +testModel.mets = {'as';'ac';'xc';'yc';'mc';'ec';'es'}; +testModel.metNames = {'a';'a';'x';'y';'m';'e';'e'}; +testModel.comps = {'s';'c'}; +testModel.compNames = testModel.comps; +testModel.metComps = [1;2;2;2;2;2;1]; +testModel.genes = {'Gm';'Gc';'Gn';'Gb';'GE';'G7'}; +testModel.grRules = {}; +testModel.rxnGeneMat = []; + +rxnsToAdd = struct(); +rxnsToAdd.rxns = {'R1';'R2';'Rm';'Rc';'Rn';'Rb';'E';'R7';'R8'}; +rxnsToAdd.equations = {'=> a[s]';... + 'a[s] <=> a[c]';... + 'a[c] <=> x[c]';... + 'x[c] <=> m[c]';... + 'a[c] <=> y[c]';... + 'y[c] <=> m[c]';... + 'm[c] => e[c]';... + 'e[c] => e[s]';... + 'e[s] =>'}; +rxnsToAdd.grRules = {'';'';'Gm';'Gc';'Gn';'Gb';'GE';'G7';''}; +evalc('testModel = addRxns(testModel,rxnsToAdd, 3);'); +testModel.c = [0;0;0;0;0;0;0;0;1]; +testModel.ub = repmat(1000,9,1); +testModel.lb = [0;-1000;-1000;-1000;-1000;-1000;0;0;0]; +testModel.rxnNames = testModel.rxns; +testModel.b = zeros(7,1); +end + function testModelRxnScores = getTstModelRxnScores() testModelRxnScores = [-2;-2;-1;7;0.5;0.5;-1;-2;-3;3.5]; end diff --git a/tests/function_tests/tManipulation.m b/tests/function_tests/tManipulation.m index e8cfa3b7..aa9ea83f 100644 --- a/tests/function_tests/tManipulation.m +++ b/tests/function_tests/tManipulation.m @@ -510,6 +510,23 @@ function mergeCompartmentsDropsStaleCompMiriams(testCase) end end + function mergeCompartmentsKeepsFirstSortedId(testCase) + % The merged metabolite keeps the id that sorts first, not the + % id of the copy that comes first in the model. + m = tManipulation.mergeToyModel({'ACCOA_m';'ACCOA_c';'X_c'}, {'c';'m'}, [2;1;1]); + evalc('m2 = mergeCompartments(m);'); + testCase.verifyEqual(sort(m2.mets), {'ACCOA_c';'X_c'}); + testCase.verifyEqual(m2.comps, {'s'}); + end + + function mergeCompartmentsPrefersCopyInMergedComp(testCase) + % A copy that is already in compartment s keeps its id, even when + % another copy has an id that sorts before it. + m = tManipulation.mergeToyModel({'ACCOA_c';'ACCOA_z';'X_c'}, {'c';'s'}, [1;2;1]); + evalc('m2 = mergeCompartments(m);'); + testCase.verifyEqual(sort(m2.mets), {'ACCOA_z';'X_c'}); + end + function mergeModelsReturnsStruct(testCase) modelB = removeReactions(testCase.model, testCase.model.rxns(1:5), true, true, true); evalc('merged = mergeModels({testCase.model; modelB});'); @@ -986,6 +1003,23 @@ function expandModelLeavesDnfAlone(testCase) methods (Static, Access = private) + function m = mergeToyModel(mets, comps, metComps) + % Two copies of acetyl-CoA (mets 1 and 2) and X (met 3): + % R1 X => copy 1, R2 copy 2 => X. + m = struct(); + m.id = 'merge'; + m.rxns = {'R1';'R2'}; + m.rxnNames = m.rxns; + m.mets = mets; + m.metNames = {'acetyl-CoA';'acetyl-CoA';'X'}; + m.metComps = metComps; + m.comps = comps; + m.compNames = comps; + m.S = sparse([1 0; 0 -1; -1 1]); + m.lb = [0;0]; m.ub = [1000;1000]; m.rev = [0;0]; m.c = [0;0]; m.b = zeros(3,1); + m.genes = {}; m.grRules = {'';''}; m.rxnGeneMat = sparse(2,0); + end + 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 diff --git a/tests/function_tests/tUtils.m b/tests/function_tests/tUtils.m index 858ed0b2..b85024ca 100644 --- a/tests/function_tests/tUtils.m +++ b/tests/function_tests/tUtils.m @@ -181,6 +181,32 @@ function grRuleToDNFDistributes(testCase) testCase.verifyEqual(numel(grRuleToDNF('(g1 or g2) and (g3 or g4)')), 4); end + function grRuleToDNFTakesTheParseGrRuleGrammar(testCase) + % Doubled symbols and operators in any case are operators; an + % operator glued to a gene id is part of that id. + testCase.verifyEqual(grRuleToDNF('G1 && (G2 || G3)'), {{'G1','G2'},{'G1','G3'}}); + testCase.verifyEqual(grRuleToDNF('G1 Or G2 aNd G3'), {{'G1'},{'G2','G3'}}); + testCase.verifyEqual(grRuleToDNF('G1&G2 or G3'), {{'G1&G2'},{'G3'}}); + end + + function grRuleToDNFWhitespaceIsIsspace(testCase) + % isspace is true for char(28) and char(8195) and false for + % the no-break space char(160) and char(8239), which are then + % part of a word. + testCase.verifyEqual(grRuleToDNF(['G1' char(28) 'or' char(8195) 'G2']), {{'G1'},{'G2'}}); + for c = [133 160 8199 8239] + testCase.verifyError(@() grRuleToDNF(['G1' char(c) 'and G2']), 'RAVEN:badGrRule'); + testCase.verifyEqual(grRuleToDNF(['G1' char(c) 'G2']), {{['G1' char(c) 'G2']}}); + end + end + + function grRuleToDNFRejectsMalformed(testCase) + rules = {'G1 and', '(G1 or G2', 'G1 G2', '()', 'G1 or or G2', 'G1 xor G2'}; + for k = 1:numel(rules) + testCase.verifyError(@() grRuleToDNF(rules{k}), 'RAVEN:badGrRule', rules{k}); + end + end + function grRuleToDNFSimpleCases(testCase) testCase.verifyEqual(grRuleToDNF(''), {}); testCase.verifyEqual(grRuleToDNF('G1'), {{'G1'}});