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