tinitTests

run this test case with the command
results = runtests('tinitTests.m')

0001 %run this test case with the command
0002 %results = runtests('tinitTests.m')
0003 function tests = tinitTests
0004 tests = functiontests(localfunctions);
0005 solverExist = [exist('gurobi','file'), exist('scip.mexw64','file')] ==3;
0006 if solverExist(1)
0007     try
0008         gurobi_read('solverTests.m'); % Random function call
0009     catch ME
0010         if ~startsWith(ME.message,'Gurobi error 10012') % Expected error code, others may indicate problems with license
0011             solverExist(1) = false;
0012         end
0013     end
0014 end
0015 if all(~solverExist)
0016     disp('No suitable solver (gurobi or scip) was found, ftINIT tests skipped.')
0017     skipTests = contains({tests.Name},'ftinit','IgnoreCase',true);
0018     tests(skipTests) = [];
0019 end
0020 end
0021 
0022 function testparsexpTask1List(testCase)
0023 sourceDir = fileparts(which(mfilename));
0024 taskStruct = parseTaskList(strcat(sourceDir, '/test_data/test_tasks.txt'));
0025 taskStructExcel = parseTaskList(strcat(sourceDir, '/test_data/test_tasks.xls'));
0026 %check that all fields in the first line are what we expect
0027 
0028 expTask1.id='ER';
0029 expTask1.description='Aerobic rephosphorylation of ATP from glucose';
0030 expTask1.shouldFail=true;
0031 expTask1.printFluxes=true;
0032 expTask1.comments='Messed up reaction';
0033 expTask1.inputs={'O2[s]';'glucose[s]'};
0034 expTask1.LBin=[23.6;23.6];
0035 expTask1.UBin=[23.8;23.8];
0036 expTask1.outputs={'H2O[s]';'CO2[s]'};
0037 expTask1.LBout=[26.1;26.1];
0038 expTask1.UBout=[26.2;26.2];
0039 expTask1.equations={'ATP[c] + H2O[c] => ADP[c] + Pi[c] + H+[c]'};
0040 expTask1.LBequ=30.2;
0041 expTask1.UBequ=30.6;
0042 expTask1.changed={'ATP[a] + H2O[a] => ADP[a] + Pi[a] + H+[a]'};
0043 expTask1.LBrxn=56.2;
0044 %check that the check works as expected
0045 verifyNotEqual(testCase,taskStruct(1),expTask1)
0046 expTask1.UBrxn=60;%now add the last
0047 
0048 
0049 verifyEqual(testCase,taskStruct(1),expTask1)
0050 verifyEqual(testCase,taskStructExcel(1),expTask1)
0051 
0052 %check that we have 2 tasks in total
0053 verifyEqual(testCase,length(taskStruct),2)
0054 verifyEqual(testCase,length(taskStructExcel),2)
0055 
0056 expTask2.id='BS';
0057 expTask2.description='ATP de novo synthesis';
0058 expTask2.shouldFail=false;
0059 expTask2.printFluxes=false;
0060 expTask2.comments='';
0061 expTask2.inputs={'O2[s]';'glucose[s]';'NH3[s]';'Pi[s]'};
0062 expTask2.LBin=[0;0;0;0];
0063 expTask2.UBin=[1000;1000;1000;1000];
0064 expTask2.outputs={'H2O[s]';'CO2[s]';'ATP[c]'};
0065 expTask2.LBout=[0;0;1];
0066 expTask2.UBout=[1.1;1.1;1.3];
0067 expTask2.equations={'ATP[c] + H2O[c] <=> ADP[c] + Pi[c] + H+[c]'};
0068 expTask2.LBequ=-1000;
0069 expTask2.UBequ=1000;
0070 expTask2.changed={};
0071 expTask2.LBrxn=[];
0072 expTask2.UBrxn=[];
0073 
0074 verifyEqual(testCase,taskStruct(2),expTask2)
0075 verifyEqual(testCase,taskStructExcel(2),expTask2)
0076 
0077 %and, check that the two formats produce the same
0078 verifyEqual(testCase,taskStruct,taskStructExcel)
0079 end
0080 
0081 
0082 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0083 % testModel
0084 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0085 
0086 function testModel = getTstModel()
0087 testModel = struct();
0088 testModel.id = 'testModel';
0089 testModel.rxns = {};
0090 testModel.S=[];
0091 testModel.rev=[];
0092 testModel.mets = {'as';'ac';'bc';'cc';'dc';'ec';'es';'fc'};
0093 testModel.metNames = {'a';'a';'b';'c';'d';'e';'e';'f'};
0094 testModel.comps = {'s';'c'};
0095 testModel.compNames = testModel.comps;
0096 testModel.metComps = [1;2;2;2;2;2;1;2];
0097 testModel.genes = {'G1';'G2';'G3';'G4';'G5';'G6';'G7';'G8';'G9';'G10'};
0098 testModel.grRules = {};
0099 testModel.rxnGeneMat = [];
0100 
0101 rxnsToAdd = struct();
0102 rxnsToAdd.rxns = {'R1';'R2';'R3';'R4';'R5';'R6';'R7';'R8';'R9';'R10'};
0103 rxnsToAdd.equations = {'=> a[s]';...
0104     'a[s] <=> a[c]';...
0105     'a[c] <=> b[c] + c[c]';...
0106     'a[c] <=> 2 d[c]';...
0107     'b[c] + c[c] => e[c]';...
0108     '2 d[c] => e[c]';...
0109     'e[c] => e[s]';...
0110     'e[s] =>';...
0111     'a[c] <=> f[c]';...
0112     'f[c] <=> e[c]'};
0113 rxnsToAdd.grRules = {'';'';'G3';'G4';'G5';'G6';'G7';'';'G9';'G10'};
0114 testModel = addRxns(testModel,rxnsToAdd, 3);
0115 testModel.c = [0;0;0;0;0;0;0;1;0;0];%optimize for output flux, if this is used, not sure
0116 testModel.ub = repmat(1000,10,1);
0117 testModel.lb = [0;-1000;-1000;-1000;0;0;0;0;-1000;-1000];
0118 testModel.rxnNames = testModel.rxns;
0119 testModel.b = repmat(0,8,1);
0120 end
0121 
0122 function testModelRxnScores = getTstModelRxnScores()
0123 testModelRxnScores = [-2;-2;-1;7;0.5;0.5;-1;-2;-3;3.5];
0124 end
0125 
0126 function testModelTasks = getTstModelTasks()
0127 testModelTasks = struct();
0128 testModelTasks.id = 'Gen e[s] from a[s]';
0129 testModelTasks.description = 'Gen e[s] from a[s]';
0130 testModelTasks.shouldFail = false;
0131 testModelTasks.printFluxes = false;
0132 testModelTasks.comments = '';
0133 testModelTasks.inputs = {'a[s]'};
0134 testModelTasks.LBin = 0;
0135 testModelTasks.UBin = inf;
0136 testModelTasks.outputs = {'e[s]'};
0137 testModelTasks.LBout = 1;
0138 testModelTasks.UBout = 1;
0139 testModelTasks.equations = {};
0140 testModelTasks.LBequ = [];
0141 testModelTasks.UBequ = [];
0142 testModelTasks.changed = {};
0143 testModelTasks.LBrxn = {};
0144 testModelTasks.UBrxn = {};
0145 end
0146 
0147 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0148 % testModel2 - not used directly for now, but indirectly
0149 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0150 
0151 function testModel2 = getTstModel2()
0152 testModel2 = struct();
0153 testModel2.id = 'testModel2';
0154 testModel2.rxns = {};
0155 testModel2.S=[];
0156 testModel2.rev=[];
0157 testModel2.mets = {'a';'b'};
0158 testModel2.metNames = {'a';'b'};
0159 testModel2.comps = {'s'};
0160 testModel2.compNames = testModel2.comps;
0161 testModel2.metComps = [1;1];
0162 testModel2.genes = {'G1';'G2';'G3';'G4'};
0163 testModel2.grRules = {};
0164 testModel2.rxnGeneMat = [];
0165 
0166 rxnsToAdd = struct();
0167 rxnsToAdd.rxns = {'R1';'R2';'R3';'R4'};
0168 rxnsToAdd.equations = {'a[s] <=>';...
0169     'a[s] => b[s]';...
0170     'a[s] <=> b[s]';...
0171     'b[s] =>'};
0172 rxnsToAdd.grRules = testModel2.genes;
0173 testModel2 = addRxns(testModel2,rxnsToAdd,3,'',true,true);
0174 testModel2.c = [0;0;0;1];%optimize for output flux, if this is used, not sure
0175 testModel2.ub = repmat(1000,4,1);
0176 testModel2.lb = [-1000;0;-1000;0];
0177 testModel2.rxnNames = testModel2.rxns;
0178 testModel2.b = zeros(2,1);
0179 end
0180 
0181 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0182 % testModel4
0183 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0184 
0185 function testModel4 = getTstModel4()
0186 testModel2 = getTstModel2();
0187 testModel4 = testModel2;
0188 
0189 rxnsToAdd = struct();
0190 rxnsToAdd.rxns = {'R5';'R6';'R7';'R8';'R9';'R10';'R11'};
0191 rxnsToAdd.equations = {'5 a[s] <=> 5 d[s]';...
0192     'e[s] <=> d[s]';
0193     'f[s] + g[s] <=> e[s]';
0194     'b[s] <=> f[s]';
0195     'h[s] <=> g[s]';
0196     'h[s] =>';
0197     'e[s] => g[s]'};
0198 rxnsToAdd.grRules = {'G5';'G6';'G7';'G8';'G9';'G10';'G11'};
0199 [~,testModel4] = evalc("addRxns(testModel4,rxnsToAdd, 3, [], true, true);");
0200 end
0201 
0202 function testModel4RxnScores = getTstModel4RxnScores()
0203 testModel4RxnScores = [-1;-1;2;-1;0.5;-2;1;1.3;-0.5;-0.4;8];
0204 end
0205 
0206 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0207 % testModel5
0208 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0209 
0210 function testModel5 = getTstModel5()
0211 testModel = getTstModel();
0212 
0213 rxnsToAdd = struct();
0214 rxnsToAdd.rxns = {'R11';'R12';'R13';'R14'};
0215 rxnsToAdd.equations = {'a[c] <=> g[c]';...
0216     'a[c] <=> g[c]';...
0217     'g[c] <=> e[c]';...
0218     'g[c] <=> e[c]'};
0219 rxnsToAdd.grRules = {'G11';'G12';'G13';'G14'};
0220 [~,testModel5] = evalc("addRxns(testModel,rxnsToAdd, 3, [], true, true);");
0221 end
0222 
0223 function testModel5RxnScores = getTstModel5RxnScores()
0224 testModel5RxnScores = [getTstModelRxnScores();-1;-1.5;-1;-1.5];
0225 end
0226 
0227 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0228 %T0001: testModel without tasks
0229 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0230 function testftINIT_T0001(testCase)
0231 solverExist = [exist('gurobi','file'), exist('scip.mexw64','file')] == 3;
0232 currSolver = getpref('RAVEN','solver');
0233 if solverExist(1)
0234     setRavenSolver('gurobi');
0235 elseif solverExist(2)
0236     setRavenSolver('scip');
0237 else
0238     error('No compatible solvers found');
0239 end
0240 %    detectedMets = {};
0241 testParams = struct();
0242 
0243 %    params.TimeLimit = 10;
0244 
0245 testModel = getTstModel();
0246 [~, prepDataTest1] = evalc('prepINITModel(testModel, {}, {}, false, {}, ''s'');');
0247 %check some things in the prepData
0248 %1. We expect 3 rxns in origRxnsToZero:
0249 verifyTrue(testCase, all(strcmp(prepDataTest1.refModel.rxns(prepDataTest1.toIgnoreExch) , {'R1';'R8'})))
0250 %note that R7 should not be there, since it has a GPR.
0251 
0252 arrayData1.genes = testModel.genes;
0253 arrayData1.tissues = {'a'};
0254 arrayData1.levels = getExprForRxnScore(getTstModelRxnScores());
0255 arrayData1.threshold = 1;
0256 
0257 [~,tst1ResModel1] = evalc('ftINIT(prepDataTest1,arrayData1.tissues{1},[],[],arrayData1,[],getINITSteps(),true,true,testParams,false);');
0258 
0259 %We expect R1 and R8 to be added since they have no GPRs and are exch rxns. The transport rxn R2 without GPR will however be removed,
0260 %since we in the standard setting run the third step with [1;0;0;0;1;0;0], meaning that such reactions will be removed
0261 %R7 however will not be added since it has a GPR
0262 verifyTrue(testCase, all(strcmp(tst1ResModel1.rxns,{'R1';'R4';'R6';'R8';'R9';'R10'})))
0263 
0264 %also test spontaneous
0265 [~, prepDataTest1] = evalc('prepINITModel(testModel, {}, {''R7'';''R10''}, false, {}, ''s'');');
0266 verifyTrue(testCase, all(strcmp(prepDataTest1.refModel.rxns(prepDataTest1.toIgnoreExch | prepDataTest1.toIgnoreSpont), {'R1';'R7';'R8';'R10'})))
0267 arrayData1.genes = testModel.genes;
0268 arrayData1.tissues = {'a'};
0269 arrayData1.levels = getExprForRxnScore(getTstModelRxnScores());
0270 arrayData1.threshold = 1;
0271 [~,tst1ResModel1] = evalc('ftINIT(prepDataTest1,arrayData1.tissues{1},[],[],arrayData1,[],getINITSteps(),true,true,testParams,false);');
0272 %the model should now change to include the "correct" path (including 'R2') and
0273 %skip R9/R10:
0274 verifyTrue(testCase, all(strcmp(tst1ResModel1.rxns,{'R1';'R2';'R4';'R6';'R7';'R8'})), 1)
0275 setRavenSolver(currSolver);
0276 end
0277 
0278 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0279 %T0002: Create a task that wants to generate e[s] from a[s] for testModel
0280 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0281 function testftINIT_T0002(testCase)
0282 solverExist = [exist('gurobi','file'), exist('scip.mexw64','file')] ==3;
0283 currSolver = getpref('RAVEN','solver');
0284 if solverExist(1)
0285     setRavenSolver('gurobi');
0286 elseif solverExist(2)
0287     setRavenSolver('scip');
0288 else
0289     error('No compatible solvers found');
0290 end
0291 testModel = getTstModel();
0292 testModelTasks = getTstModelTasks();
0293 testRxnScores = getTstModelRxnScores();
0294 testParams = struct();
0295 [~, prepDataTest1] = evalc('prepINITModel(testModel, testModelTasks, {}, false, {}, ''s'');');
0296 %We now expect to R2 and R7 to be essential. Note that R1 and R8 are not essential,
0297 %the exchange rxns are not used when checking tasks.
0298 %This is a bit complicated to check, because the essential rxns are expressed
0299 %as rxn ids of the linearly merged model. We expect those to be called R1 and R7:
0300 verifyTrue(testCase, all(strcmp(prepDataTest1.essentialRxns,{'R1';'R7'})))
0301 
0302 arrayData1.genes = testModel.genes;
0303 arrayData1.tissues = {'a'};
0304 arrayData1.levels = getExprForRxnScore(testRxnScores);
0305 arrayData1.threshold = 1;
0306 [~,tst1ResModel1] = evalc('ftINIT(prepDataTest1,arrayData1.tissues{1},[],[],arrayData1,[],getINITSteps(),true,true,testParams,false);');
0307 %Since both R2 and R7 are now essential, we expect all rxns to be on except R3 and
0308 %R5 (which have a negative total score and are not needed for the task)
0309 verifyTrue(testCase, all(strcmp(tst1ResModel1.rxns,{'R1';'R2';'R4';'R6';'R7';'R8';'R9';'R10'})))
0310 setRavenSolver(currSolver);
0311 end
0312 
0313 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0314 %T0003: The second step - gapfilling
0315 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0316 function testftINIT_T0003(testCase)
0317 %First generate a model with gaps. We can use testModel. First we add
0318 %boundary mets. Then we remove the exchange reactions
0319 %(which is required for filling gaps) and create a gap by removing the R7 reaction.
0320 testModel = getTstModel();
0321 testModelTasks = getTstModelTasks();
0322 testRxnScores = getTstModelRxnScores();
0323 
0324 mTempRef = closeModel(testModel);
0325 mTempRef = removeReactions(mTempRef, {'R1';'R8'});
0326 mTemp = removeReactions(mTempRef, {'R7'});
0327 mTemp.id = 'tmp';
0328 tmpRxnScores = testRxnScores([2;3;4;5;6;7;9;10]);
0329 %now check that R7 is added back
0330 [~,outModel,addedRxnMat] = evalc('ftINITFillGapsForAllTasks(mTemp,mTempRef,[],false,min(tmpRxnScores,-0.1),testModelTasks);');
0331 verifyTrue(testCase, all(strcmp(mTempRef.rxns(addedRxnMat),'R7')))%ok
0332 end
0333 
0334 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0335 %T0004: MergeLinear and groupRxnScores
0336 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0337 function testftINIT_T0004(testCase)
0338 %first testModel
0339 testModel = getTstModel();
0340 testRxnScores = getTstModelRxnScores();
0341 [reducedModel,origRxnIds,groupIds,reversedRxns]=mergeLinear(testModel, {});
0342 %We expect mergeLinear to merge {R1,R2}, {R3,R5}, {R4,R6}, {R7,R8}, {R9,R10}
0343 verifyTrue(testCase, all(groupIds == [1;1;2;3;2;3;4;4;5;5]))
0344 %we expect R1, R3, R4, R7 to be irreversible, R9 to be reversible
0345 verifyTrue(testCase, all(reducedModel.rev == [0;0;0;0;1]))
0346 verifyTrue(testCase, all(reducedModel.lb == [0;0;0;0;-1000]))
0347 
0348 newRxnScores=groupRxnScores(reducedModel, testRxnScores, origRxnIds, groupIds, ismember(origRxnIds, {'R1';'R2';'R8'}));
0349 verifyTrue(testCase, all(newRxnScores == [0;-0.5;7.5;-1;0.5]))
0350 
0351 %then testModel4
0352 testModel4 = getTstModel4();
0353 [reducedModel,origRxnIds,groupIds,reversedRxns]=mergeLinear(testModel4, {});
0354 %we expect {R5,R6},{R7,R8}, and{R9,R10} to be merged
0355 verifyTrue(testCase, all(groupIds == [0;0;0;0;1;1;2;2;3;3;0]))
0356 %check reversibility
0357 verifyTrue(testCase, all(reducedModel.rev == [1;0;1;0;1;1;0;0]))
0358 %check that some reactions have flipped direction when turned to irrev
0359 verifyTrue(testCase, strcmp(constructEquations(reducedModel, 'R9'),'g[s] => '))
0360 verifyTrue(testCase, all(find(reversedRxns) == [6;9]))
0361 end
0362 
0363 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0364 %T0006: reverseRxns
0365 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0366 
0367 function testftINIT_T0006(testCase)
0368 %R1 = '=> a[s]';...
0369 %R3 = 'a[c] <=> b[c] + c[c]';...
0370 testModel = getTstModel();
0371 
0372 tmpModel = reverseRxns(testModel, {'R1';'R3'});
0373 res = constructEquations(tmpModel, {'R1';'R3'});
0374 expRes = {'a[s] => ';'b[c] + c[c] <=> a[c]'};
0375 verifyTrue(testCase, all(strcmp(res,expRes)))
0376 end
0377 
0378 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0379 %T0007: rescaleModelForINIT
0380 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0381 function testftINIT_T0007(testCase)
0382 miniModel = struct();
0383 miniModel.S = [1,1000;-1,-40];
0384 miniModel.rxns = {'1';'2'};
0385 miniModel.mets = {'1';'2'};
0386 res = rescaleModelForINIT(miniModel,10);
0387 verifyTrue(testCase, abs(res.S(1,2) - res.S(2,2)*-10) < 10^-6)
0388 verifyTrue(testCase, abs((abs(res.S(1,2)) + abs(res.S(2,2)))/2) - 1 < 10^-6)
0389 end
0390 
0391 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0392 %T0008: testModel with metabolomics
0393 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0394 function testftINIT_T0008(testCase)
0395 solverExist = [exist('gurobi','file'), exist('scip.mexw64','file')] ==3;
0396 currSolver = getpref('RAVEN','solver');
0397 if solverExist(1)
0398     setRavenSolver('gurobi');
0399 elseif solverExist(2)
0400     setRavenSolver('scip');
0401 else
0402     error('No compatible solvers found');
0403 end
0404 
0405 testParams = struct();
0406 
0407 testModel = getTstModel();
0408 [~, prepDataTest1] = evalc('prepINITModel(testModel, {}, {}, false, {}, ''s'');');
0409 
0410 arrayData1.genes = testModel.genes;
0411 arrayData1.tissues = {'a'};
0412 arrayData1.levels = getExprForRxnScore(getTstModelRxnScores());
0413 arrayData1.threshold = 1;
0414 
0415 [~,tst1ResModel1] = evalc('ftINIT(prepDataTest1,arrayData1.tissues{1},[],[],arrayData1,[],getINITSteps(),true,true,testParams,false);');
0416 
0417 %First the same as in T0001 (we keep them here to make the test case understandable):
0418 %We expect R1 and R8 to be added since they have no GPRs and are exch rxns. The transport rxn R2 without GPR will however be removed,
0419 %since we in the standard setting run the third step with [1;0;0;0;1;0;0], meaning that such reactions will be removed
0420 %R7 however will not be added since it has a GPR
0421 verifyTrue(testCase, all(strcmp(tst1ResModel1.rxns,{'R1';'R4';'R6';'R8';'R9';'R10'})))
0422 
0423 %make R7 and R10 spontaneous (also same as in T0001)
0424 [~, prepDataTest1] = evalc('prepINITModel(testModel, {}, {''R7'';''R10''}, false, {}, ''s'');');
0425 verifyTrue(testCase, all(strcmp(prepDataTest1.refModel.rxns(prepDataTest1.toIgnoreExch | prepDataTest1.toIgnoreSpont), {'R1';'R7';'R8';'R10'})))
0426 [~,tst1ResModel1] = evalc('ftINIT(prepDataTest1,arrayData1.tissues{1},[],[],arrayData1,[],getINITSteps(),true,true,testParams,false);');
0427 %the model should now change to include the "correct" path (including 'R2') and
0428 %skip R9/R10:
0429 verifyTrue(testCase, all(strcmp(tst1ResModel1.rxns,{'R1';'R2';'R4';'R6';'R7';'R8'})))
0430 %now, test to add the metabolite f - this should turn back the favor to R9/R10:
0431 [~,tst1ResModel1] = evalc('ftINIT(prepDataTest1,arrayData1.tissues{1},[],[],arrayData1,{''f''},getINITSteps(),true,true,testParams,false);');
0432 %R9 should now be included, the presence of R2 is random
0433 if length(tst1ResModel1.rxns) == 7
0434     verifyTrue(testCase, all(strcmp(tst1ResModel1.rxns,{'R1';'R4';'R6';'R7';'R8';'R9';'R10'})))
0435 else
0436     verifyTrue(testCase, all(strcmp(tst1ResModel1.rxns,{'R1';'R2';'R4';'R6';'R7';'R8';'R9';'R10'})))
0437 end
0438 %now, test to add the metabolite a, e, f - this should give the same result:
0439 [~,tst1ResModel1] = evalc('ftINIT(prepDataTest1,arrayData1.tissues{1},[],[],arrayData1,{''f'';''a'';''e''},getINITSteps(),true,true,testParams,false);');
0440 %Should be the same as above (R2 is random)
0441 if length(tst1ResModel1.rxns) == 7
0442     verifyTrue(testCase, all(strcmp(tst1ResModel1.rxns,{'R1';'R4';'R6';'R7';'R8';'R9';'R10'})))
0443 else
0444     verifyTrue(testCase, all(strcmp(tst1ResModel1.rxns,{'R1';'R2';'R4';'R6';'R7';'R8';'R9';'R10'})))
0445 end
0446 
0447 %now, test to add the metabolite b - this should turn on R2 and R3/R5 and turn off R9:
0448 [~,tst1ResModel1] = evalc('ftINIT(prepDataTest1,arrayData1.tissues{1},[],[],arrayData1,{''b''},getINITSteps(),true,true,testParams,false);');
0449 verifyTrue(testCase, all(strcmp(tst1ResModel1.rxns,{'R1';'R2';'R3';'R4';'R5';'R6';'R7';'R8'})))
0450 
0451 %now on model 5 to test reactions that are not linearly merged (testModel has only merged rxns)
0452 testModel5 = getTstModel5();
0453 arrayData1.genes = testModel5.genes;
0454 arrayData1.tissues = {'a'};
0455 arrayData1.levels = getExprForRxnScore(getTstModel5RxnScores());
0456 arrayData1.threshold = 1;
0457 
0458 [~, prepDataTest5] = evalc('prepINITModel(testModel5, {}, {''R7'';''R10''}, false, {}, ''s'');');
0459 [~,tst1ResModel1] = evalc('ftINIT(prepDataTest5,arrayData1.tissues{1},[],[],arrayData1,{},getINITSteps(),true,true,testParams,false);');
0460 %We expect the 'true' path, i.e. through R2, not R9/R10 or R11-R14
0461 verifyTrue(testCase, all(strcmp(tst1ResModel1.rxns,{'R1';'R2';'R4';'R6';'R7';'R8'})), 1)
0462 
0463 %now add metabolite g
0464 %modify the scores a bit
0465 [~, prepDataTest5] = evalc('prepINITModel(testModel5, {}, {''R10''}, false, {}, ''s'');');
0466 arrayData1.levels(7) = getExprForRxnScore(-1.1); %modify to avoid randomness
0467 [~,tst1ResModel1] = evalc('ftINIT(prepDataTest5,arrayData1.tissues{1},[],[],arrayData1,{''g''},getINITSteps(),true,true,testParams,false);');
0468 %We expect R2 to be replaced with R11 and R13
0469 verifyTrue(testCase, all(strcmp(tst1ResModel1.rxns,{'R1';'R4';'R6';'R8';'R11';'R13'})), 1)
0470 setRavenSolver(currSolver);
0471 end
0472 
0473 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0474 %T0009: getExprFromRxnScore
0475 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0476 function testftINIT_T0009(testCase)
0477 testModel = getTstModel();
0478 [~, prepDataTest1] = evalc('prepINITModel(testModel, {}, {}, false, {}, ''s'');');
0479 
0480 arrayData1.genes = testModel.genes;
0481 arrayData1.tissues = {'a'};
0482 arrayData1.levels = getExprForRxnScore(getTstModelRxnScores());
0483 arrayData1.threshold = 1;
0484 
0485 rxnScores = scoreComplexModel(prepDataTest1.refModel,[],arrayData1,arrayData1.tissues{1},[]);
0486 expRes = getTstModelRxnScores();
0487 verifyTrue(testCase, all(abs(rxnScores - expRes) < 10^-10)) %ok
0488 end
0489 
0490 function testModelL = getTstModelL()
0491 testModelL = struct();
0492 testModelL.id = 'testModel';
0493 testModelL.rxns = {};
0494 testModelL.S=[];
0495 testModelL.rev=[];
0496 testModelL.metNames = {'e1';'e2';'e3';'e4';'e5';'e6';'e7';'e8';'e9';'e1';'e2';'e3';'e4';'e5';'e6';'e7';'e8';'e9';'x1';'x2';'x3';'x4';'x5';'x6';'x7';'x8';'x9';'x10';'x11'};
0497 testModelL.comps = {'s';'c'};
0498 testModelL.compNames = testModelL.comps;
0499 testModelL.metComps = [1;1;1;1;1;1;1;1;1;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2];
0500 testModelL.mets = strcat(testModelL.metNames, testModelL.comps(testModelL.metComps));
0501 
0502 testModelL.grRules = {};
0503 testModelL.rxnGeneMat = [];
0504 
0505 testModelL.genes = {'Ge1';'Ge2';'Ge4';'Ge5';'Ge7';'Ge9'; 'Gr1';'Gr2';'Gr3';'Gr5';'Gr6';'Gr7';'Gr8';'Gr9';'Gr10';'Gr11';'Gr12';'Gr14';'Gr15'};
0506 
0507 testModelL.ub = [];
0508 testModelL.lb = [];
0509 
0510 rxnsToAdd = struct();
0511 rxnsToAdd.rxns = {  'S1';'S2';'S3';'S4';'S5';'S6';'S7';'S8';'S9';'E1';'E2';'E2b';'E3';'E4';'E5';'E6';'E7';'E8';'E9';'R1';'R2';'R3';'R4';'R5';'R6';'R7';'R8';'R9';'R10';'R11';'R12';'R13';'R14';'R15'};
0512 rxnsToAdd.grRules = {'';  '';  '';  '';  '';  '';  '';  '';  ''; 'Ge1';'Ge2';'';'';'Ge4';'Ge5';'';'Ge7';'';'Ge9'; 'Gr1';'Gr2';'Gr3';'';'Gr5';'Gr6';'Gr7';'Gr8';'Gr9';'Gr10';'Gr11';'Gr12';'';'Gr14';'Gr15'};
0513 rxnsToAdd.equations = {'e1[s] <=>';...
0514     'e2[s] <=>';...
0515     'e3[s] <=>';...
0516     'e4[s] <=>';...
0517     'e5[s] <=>';...
0518     'e6[s] <=>';...
0519     'e7[s] <=>';...
0520     'e8[s] <=>';...
0521     'e9[s] <=>';...
0522     'e1[s] <=> e1[c]';...
0523     'e2[s] <=> e2[c]';...
0524     'e2[s] <=> e2[c]';... %b variant
0525     'e3[s] <=> e3[c]';...
0526     'e4[s] <=> e4[c]';...
0527     'e5[s] <=> e5[c]';...
0528     'e6[s] <=> e6[c]';...
0529     'e7[s] <=> e7[c]';...
0530     'e8[s] <=> e8[c]';...
0531     'e9[s] <=> e9[c]';...
0532     'e1[c] + e2[c] <=> x1[c]';... %R1
0533     'e1[c] + e3[c] => x2[c] + x3[c]';... %R2
0534     'e4[c] + x3[c] => x4[c] + x5[c]';... %R3
0535     'e5[c] + e6[c] + x4[c] => 2 x2[c] + x6[c]';... %R4
0536     'x1[c] + x2[c] <=> x7[c] + 2 x8[c]';... %R5
0537     'x2[c] + x8[c] => x3[c] + x9[c]';... %R6
0538     'x4[c] <=> x9[c]';... %R7
0539     'x5[c] <=> x9[c]';... %R8
0540     'x6[c] <=> x10[c]';... %R9
0541     'x6[c] <=> x11[c]';... %R10
0542     'x10[c] + 2 x11[c] => e7[c]';... %R11
0543     'x9[c] + x10[c] <=> e8[c]';... %R12
0544     'x7[c] + x8[c] + x9[c] => e9[c]';... %R13
0545     'x6[c] => x9[c]';... %R14
0546     'x3[c] => x9[c]'... %R15
0547     };
0548 testModelL = addRxns(testModelL,rxnsToAdd, 3);
0549 testModelL.c = [0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];%optimize for output flux, if this is used, not sure
0550 testModelL.rxnNames = testModelL.rxns;
0551 testModelL.b = repmat(0,length(testModelL.mets),1);
0552 end
0553 
0554 function testModelLGeneScores = getTstModelLGeneScores()
0555 %testModelL.genes = {'Ge1';'Ge2';'Ge4';'Ge5';'Ge7';'Ge9'; 'Gr1';'Gr2';'Gr3';'Gr5';'Gr6';'Gr7';'Gr8';'Gr9';'Gr10';'Gr11';'Gr12';'Gr14';'Gr15'};
0556 testModelLGeneScores = [3; -1;   8;    6;    -5;    5;     4;    5;    2;    3;    6;    1;    3;    1;    -3;    1;     3;      1;    2];
0557 end
0558 
0559 
0560 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0561 %T0050: Test ftINIT on a slightly larger and more complex model
0562 %       Specifically tests if the three-step variant and the "full"
0563 %       variant gives similar results
0564 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0565 function testftINIT_T0050(testCase)
0566 solverExist = [exist('gurobi','file'), exist('scip.mexw64','file')] ==3;
0567 currSolver = getpref('RAVEN','solver');
0568 if solverExist(1)
0569     setRavenSolver('gurobi');
0570 elseif solverExist(2)
0571     setRavenSolver('scip');
0572 else
0573     error('No compatible solvers found');
0574 end
0575 
0576 testModelL = getTstModelL();
0577 testModelLGeneScores = getTstModelLGeneScores();
0578 testParams = struct();
0579 
0580 arrayDataL = struct();
0581 arrayDataL.genes = testModelL.genes;
0582 arrayDataL.tissues = {'t1'};
0583 arrayDataL.levels = getExprForRxnScore(testModelLGeneScores,1);
0584 arrayDataL.threshold = 1;
0585 
0586 %Run prep data
0587 [~, prepDataL] = evalc('prepINITModel(testModelL, [], {}, false, {}, ''s'');');
0588 
0589 [~,mres] = evalc('ftINIT(prepDataL,arrayDataL.tissues{1},[],[],arrayDataL,[],getINITSteps(),true,true,testParams,false);');
0590 [~,mres2] = evalc('ftINIT(prepDataL,arrayDataL.tissues{1},[],[],arrayDataL,[],getINITSteps([], ''full''),true,true,testParams,false);');
0591 
0592 expResult = {  'S1';'S2';'S3';'S4';'S5';'S6';'S7';'S8';'S9';'E1';'E2';'E3';'E4';'E5';'E6';'E8';'E9';'R1';'R2';'R3';'R4';'R5';'R6';'R7';'R8';'R9';'R12';'R13';'R14';'R15'};
0593 
0594 verifyTrue(testCase, all(contains(mres.rxns,expResult)))
0595 verifyTrue(testCase, all(contains(mres2.rxns,expResult)))
0596 
0597 %run the old tINIT version (in Human-GEM)
0598 %this is just to show that they become different, not really part of the test case
0599 %paramsL2 = struct();
0600 %paramsL2.TimeLimit = 1000;
0601 %testModelL2 = closeModel(testModelL);
0602 %init_modelOrig = getINITModel2(testModelL2,arrayDataL.tissues{1},[],[],arrayDataL,[],true,[],true,true,[],paramsL2);
0603 
0604 %in this call, I have modified the code - the possibility to turn off met secretion + don't allow flux in both directions is not possible.
0605 %The following line, around line 337, is changed
0606 %from:
0607 %[~, deletedRxnsInINIT, metProduction] = runINIT(simplifyModel(cModel),rxnScores,metabolomicsData,essentialRxnsForTasks,0,true,false,params);
0608 %to:
0609 %[~, deletedRxnsInINIT, metProduction] = runINIT(simplifyModel(cModel),rxnScores,metabolomicsData,essentialRxnsForTasks,0,false,true,params);
0610 %init_modelOrigNoSecrOneDirOnly = getINITModel2(testModelL2,arrayDataL.tissues{1},[],[],arrayDataL,[],true,[],true,true,[],paramsL2);
0611 
0612 %The models init_modelOrigNoSecrOneDirOnly and mres2 are very similar, (only one exch rxn differ, which is expected)
0613 %init_modelOrig is quite different, with a lot of gaps, and worse. So, the conclusion is that the new version does a pretty good job.
0614 setRavenSolver(currSolver);
0615 end