diff --git a/INIT/ftINITInternalAlg.m b/INIT/ftINITInternalAlg.m index 19c0ddc4..0b990920 100644 --- a/INIT/ftINITInternalAlg.m +++ b/INIT/ftINITInternalAlg.m @@ -392,10 +392,17 @@ metVarInd = (1:nMetabolMets) + (length(prob.vartype) - nMetVars); if allowExcretion - prob.csense = [repmat('L', 1, length(milpModel.mets)), ... + %S*v >= 0: a metabolite may be produced in excess of what is consumed, + %and the surplus is implicitly excreted. 'L' would say S*v <= 0, which + %lets a metabolite be consumed without ever being produced -- free + %uptake, the opposite of excretion, and enough for the MILP to switch a + %reaction on by conjuring its substrates out of nothing. + prob.csense = [repmat('G', 1, length(milpModel.mets)), ... repmat('E', 1, length(prob.b) - length(milpModel.mets))]; else - prob.csense = '='; + %One character per row: optimizeProb maps csense elementwise for glpk and + %cobra, so a scalar '=' only happens to work on gurobi. + prob.csense = repmat('E', 1, length(prob.b)); end params.intTol = 10^-7; %This value is very important. If set too low diff --git a/testing/function_tests/tINIT.m b/testing/function_tests/tINIT.m index 5c77724c..71c5e1ac 100644 --- a/testing/function_tests/tINIT.m +++ b/testing/function_tests/tINIT.m @@ -259,6 +259,32 @@ function ftINITFullVsThreeStepRuns(testCase) testCase.verifyTrue(all(contains(mres2.rxns, expResult))); end + function ftINITSeriesVariantsRun(testCase) + % Only '1+1' and 'full' were ever exercised, so the 2-step series + % from the paper -- and the allowExcretion constraint they lean on + % -- had no coverage at all. + testCase.assumeMILPSolver(); + testModel = getTstModel(); + testParams = struct(); + evalc('prepData = prepINITModel(testModel, {}, {}, false, {}, ''s'');'); + arrayData.genes = testModel.genes; + arrayData.tissues = {'a'}; + arrayData.levels = getExprForRxnScore(getTstModelRxnScores()); + arrayData.threshold = 1; + + steps = getINITSteps([], '2+1'); + evalc(['resModel = ftINIT(prepData,arrayData.tissues{1},[],[],' ... + 'arrayData,[],steps,true,true,testParams,false);']); + % Same answer as the '1+1' series in ftINITPipelineRuns. + testCase.verifyEqual(resModel.rxns, {'R1';'R4';'R6';'R8';'R9';'R10'}); + + steps = getINITSteps([], '2+0'); + evalc(['resModel = ftINIT(prepData,arrayData.tissues{1},[],[],' ... + 'arrayData,[],steps,true,true,testParams,false);']); + % '2+0' skips step 3, so the GPR-less transport R2 survives. + testCase.verifyEqual(resModel.rxns, {'R1';'R2';'R4';'R6';'R8';'R9';'R10'}); + end + end end