diff --git a/gapfilling/gapFillTopological.m b/gapfilling/gapFillTopological.m index ad5d23ee..721b6bc8 100644 --- a/gapfilling/gapFillTopological.m +++ b/gapfilling/gapFillTopological.m @@ -31,7 +31,8 @@ % -------------------------------- % 'seeds' (default []) % Metabolite IDs (cell array) available from the medium. Default: -% metabolites involved in uptake exchange reactions (lb < 0). +% 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. @@ -67,9 +68,22 @@ % ---- Identify seed metabolites ---- if isempty(seeds) - % Default: metabolites in exchange reactions that allow uptake (lb < 0) - [~, exchIdx] = getExchangeRxns(model); - uptakeIdx = exchIdx(model.lb(exchIdx) < 0); + % Default: metabolites that an exchange reaction can supply. Ask + % getExchangeRxns rather than testing lb < 0 here: RAVEN writes an uptake + % exchange either as "=> met" (addExchangeRxns 'in', so lb = 0 and a + % positive flux supplies the metabolite) or as "met =>" with lb < 0. Only + % the second has lb < 0, so a hand-rolled lb < 0 test silently misses + % every 'in' exchange. + % + % 'uptake' alone is not enough either: it means the bounds allow *only* + % uptake, so a reversible exchange ("met <=>") is classified 'reverse' + % even though it supplies the metabolite just as well. The union of the + % two is exactly "forward flux produces it, or reverse flux does", which + % is what analyse_topology's lower_bound < 0 test means under cobra's + % convention, where every exchange is written "met <=>". + [~, onlyUptakeIdx] = getExchangeRxns(model, 'uptake'); + [~, reversibleIdx] = getExchangeRxns(model, 'reverse'); + uptakeIdx = union(onlyUptakeIdx, reversibleIdx); if isempty(uptakeIdx) warning(['gapFillTopological: no uptake exchange reactions found. ' ... 'Provide seeds manually via the ''seeds'' option.']); @@ -150,6 +164,20 @@ qHead = 1; reachable(seedMetIdx) = true; +% Reactions with no substrates have nothing to wait for and fire straight +% away. They are in no subOf list, so the countdown below never reaches them +% and their products would otherwise stay unreachable forever. This is where +% "=> met" uptake exchanges enter the scope. +for j = find(subCountFwd == 0)' + firedFwd(j) = true; + for p = rxnProds{j} + if ~reachable(p) + reachable(p) = true; + queue(end+1) = p; %#ok + end + end +end + while qHead <= numel(queue) m = queue(qHead); qHead = qHead + 1; diff --git a/testing/function_tests/tGapfilling.m b/testing/function_tests/tGapfilling.m index 0514c072..b0da89b0 100644 --- a/testing/function_tests/tGapfilling.m +++ b/testing/function_tests/tGapfilling.m @@ -140,6 +140,37 @@ function gapFillTopologicalIdentifiesGaps(testCase) testCase.verifyNumElements(result.reachableMets, numel(gapModel.mets)); end + function gapFillTopologicalReachesThroughInExchange(testCase) + % A -> B -> C fed by an 'in' exchange on A. addExchangeRxns writes + % 'in' as "=> A", which has lb = 0, so a seed test of lb < 0 finds + % no uptake at all and every metabolite comes back blocked. + m = testCase.chainModel(); + m = addExchangeRxns(m, 'in', {'A'}); + evalc(['result = gapFillTopological(m, m, ''targets'', ' ... + '{''A'',''B'',''C''}, ''verbose'', false);']); + testCase.verifyTrue(all(result.reachableMets)); + testCase.verifyEmpty(result.blockedMets); + end + + function gapFillTopologicalReachesThroughReversibleExchange(testCase) + % "A <=>" supplies A too, but is 'reverse' rather than 'uptake'. + m = testCase.chainModel(); + m = addExchangeRxns(m, 'both', {'A'}); + evalc(['result = gapFillTopological(m, m, ''targets'', ' ... + '{''A'',''B'',''C''}, ''verbose'', false);']); + testCase.verifyTrue(all(result.reachableMets)); + end + + function gapFillTopologicalBlocksWithoutExchange(testCase) + % Without any exchange nothing is producible: the positive control + % for the two tests above. + m = testCase.chainModel(); + evalc(['result = gapFillTopological(m, m, ''seeds'', {}, ''targets'', ' ... + '{''A'',''B'',''C''}, ''verbose'', false);']); + testCase.verifyFalse(any(result.reachableMets)); + testCase.verifyNumElements(result.blockedMets, 3); + end + function gapFillMILPRepairsGrowth(testCase) testCase.assumeMILPSolver(); modelDB = testCase.model; modelDB.id = 'DB'; @@ -190,4 +221,18 @@ function fillGapsDispatchesGapfillMILP(testCase) end end + + methods (Access = private) + function m = chainModel(~) + % A -> B -> C, no exchanges; the caller adds the ones it needs. + m = struct(); + m.rxns = {'R1';'R2'}; m.rxnNames = {'R1';'R2'}; + m.mets = {'A';'B';'C'}; m.metNames = {'A';'B';'C'}; m.metComps = [1;1;1]; + m.comps = {'c'}; m.compNames = {'c'}; + m.S = sparse([-1 0; 1 -1; 0 1]); + m.lb = [0;0]; m.ub = [1000;1000]; m.rev = [0;0]; m.c = [0;1]; + m.b = [0;0;0]; + m.genes = {}; m.grRules = {'';''}; m.rxnGeneMat = sparse(2,0); + end + end end