checkTasks

 checkTasks
   Performs a set of simulations as defined in a task file.

   Input:
   model           a model structure
   inputFile       a task list in Excel format. See the function
                   parseTaskList for details (opt if taskStructure is
                   supplied)
   printOutput     true if the results of the test should be displayed
                   (opt, default true)
   printOnlyFailed true if only tasks that failed should be displayed
                   (opt, default false)
   getEssential    true if the essential reactions should be calculated for
                   all the tasks. This option is used with runINIT (opt,
                   default false)
   taskStructure   structure with the tasks, as from parseTaskList. If
                   this is supplied then inputFile is ignored (opt)


   Output:
   taskReport          structure with the results
       id              cell array with the id of the task
       description     cell array with the description of the task
       ok              boolean array with true if the task was successful
   essentialRxns       MxN matrix with the essential reactions (M) for each
                       task (N). An element is true if the corresponding
                       reaction is essential in the corresponding task.
                       Failed tasks and SHOULD FAIL tasks are ignored.
                       This is used by the INIT algorithm (if tasks
                       are supplied). If getEssential=false then
                       essentialRxns=false(nRxns,nTasks)
   taskStructure       structure with the tasks, as from parseTaskList
   essentialFluxes     The fluxes of the essential rxns - same structure as essentialRxns
   

   This function is used for defining a set of tasks for a model to
   perform. The tasks are defined by defining constraints on the model,
   and if the problem is feasible, then the task is considered successful.
   In general, each row can contain one constraint on uptakes, one
   constraint on outputs, one new equation, and one change of reaction
   bounds. If more bounds are needed to define the task, then several rows
   can be used for each task.

   Usage: [taskReport, essentialReactions, taskStructure]=checkTasks(model,inputFile,...
           printOutput,printOnlyFailed,getEssential,taskStructure)

0001 function [taskReport, essentialRxns, taskStructure, essentialFluxes]=checkTasks(model,inputFile,printOutput,printOnlyFailed,getEssential,taskStructure)
0002 % checkTasks
0003 %   Performs a set of simulations as defined in a task file.
0004 %
0005 %   Input:
0006 %   model           a model structure
0007 %   inputFile       a task list in Excel format. See the function
0008 %                   parseTaskList for details (opt if taskStructure is
0009 %                   supplied)
0010 %   printOutput     true if the results of the test should be displayed
0011 %                   (opt, default true)
0012 %   printOnlyFailed true if only tasks that failed should be displayed
0013 %                   (opt, default false)
0014 %   getEssential    true if the essential reactions should be calculated for
0015 %                   all the tasks. This option is used with runINIT (opt,
0016 %                   default false)
0017 %   taskStructure   structure with the tasks, as from parseTaskList. If
0018 %                   this is supplied then inputFile is ignored (opt)
0019 %
0020 %
0021 %   Output:
0022 %   taskReport          structure with the results
0023 %       id              cell array with the id of the task
0024 %       description     cell array with the description of the task
0025 %       ok              boolean array with true if the task was successful
0026 %   essentialRxns       MxN matrix with the essential reactions (M) for each
0027 %                       task (N). An element is true if the corresponding
0028 %                       reaction is essential in the corresponding task.
0029 %                       Failed tasks and SHOULD FAIL tasks are ignored.
0030 %                       This is used by the INIT algorithm (if tasks
0031 %                       are supplied). If getEssential=false then
0032 %                       essentialRxns=false(nRxns,nTasks)
0033 %   taskStructure       structure with the tasks, as from parseTaskList
0034 %   essentialFluxes     The fluxes of the essential rxns - same structure as essentialRxns
0035 %
0036 %
0037 %   This function is used for defining a set of tasks for a model to
0038 %   perform. The tasks are defined by defining constraints on the model,
0039 %   and if the problem is feasible, then the task is considered successful.
0040 %   In general, each row can contain one constraint on uptakes, one
0041 %   constraint on outputs, one new equation, and one change of reaction
0042 %   bounds. If more bounds are needed to define the task, then several rows
0043 %   can be used for each task.
0044 %
0045 %   Usage: [taskReport, essentialReactions, taskStructure]=checkTasks(model,inputFile,...
0046 %           printOutput,printOnlyFailed,getEssential,taskStructure)
0047 
0048 if nargin<3 || isempty(printOutput)
0049     printOutput=true;
0050 end
0051 if nargin<4 || isempty(printOnlyFailed)
0052     printOnlyFailed=false;
0053 end
0054 if nargin<5 || isempty(getEssential)
0055     getEssential=false;
0056 end
0057 
0058 %Prepare the input model a little
0059 model.b=zeros(numel(model.mets),2);
0060 
0061 modelMets=upper(strcat(model.metNames,'[',model.comps(model.metComps),']'));
0062 if ~isfield(model,'unconstrained')
0063     EM='Exchange metabolites should normally not be removed from the model when using checkTasks. Inputs and outputs are defined in the task file instead. Use importModel(file,false) to import a model with exchange metabolites remaining';
0064     dispEM(EM,false);
0065 end
0066 
0067 %Parse the task file
0068 if nargin<6
0069     taskStructure=parseTaskList(inputFile);
0070 end
0071 
0072 essentialRxns=false(numel(model.rxns),numel(taskStructure));
0073 essentialFluxes = NaN(numel(model.rxns),numel(taskStructure));
0074 
0075 tModel=model;
0076 taskReport=[];
0077 for i=1:numel(taskStructure)
0078     taskReport.id{i,1}=taskStructure(i).id;
0079     taskReport.description{i,1}=taskStructure(i).description;
0080     %Set the inputs
0081     if ~isempty(taskStructure(i).inputs)
0082         [I, J]=ismember(upper(taskStructure(i).inputs),modelMets);
0083         J=J(I); %Only keep the ones with matches
0084         K=ismember(upper(taskStructure(i).inputs),'ALLMETS');
0085         L=~cellfun('isempty',strfind(upper(taskStructure(i).inputs),'ALLMETSIN'));
0086         %Check that all metabolites are either real metabolites or
0087         %ALLMETS/ALLMETSIN
0088         if ~all(I|K|L)
0089             fprintf(['ERROR: Could not find all inputs in "[' taskStructure(i).id '] ' taskStructure(i).description '"\n']);
0090             taskReport.ok(i,1)=false;
0091             tModel=model;
0092             continue;
0093         end
0094         if numel(J)~=numel(unique(J))
0095             EM=['The constraints on some input(s) in "[' taskStructure(i).id '] ' taskStructure(i).description '" are defined more than one time'];
0096             dispEM(EM);
0097         end
0098         %If all metabolites should be added
0099         if any(K)
0100             %Check if ALLMETS is the first metabolite. Otherwise print a
0101             %warning since it will write over any other constraints that
0102             %are set
0103             if K(1)==0
0104                 EM=['ALLMETS is used as an input in "[' taskStructure(i).id '] ' taskStructure(i).description '" but it it not the first metabolite in the list. Constraints defined for the metabolites before it will be over-written'];
0105                 dispEM(EM,false);
0106             end
0107             %Use the first match of ALLMETS. There should only be one, but
0108             %still..
0109             tModel.b(:,1)=taskStructure(i).UBin(find(K,1))*-1;
0110         end
0111         %If metabolites in a specific compartment should be used
0112         if any(L)
0113             L=find(L);
0114             for j=1:numel(L)
0115                 %The compartment defined
0116                 compartment=upper(taskStructure(i).inputs{L(j)}(11:end-1));
0117                 %Check if it exists in the model
0118                 C=find(ismember(upper(model.comps),compartment));
0119                 if any(C)
0120                     %Match to metabolites
0121                     tModel.b(model.metComps==C,1)=taskStructure(i).UBin(L(j))*-1;
0122                 else
0123                     EM=['The compartment defined for ALLMETSIN in "[' taskStructure(i).id '] ' taskStructure(i).description '" does not exist'];
0124                     dispEM(EM);
0125                 end
0126             end
0127         end
0128         %Then add the normal constraints
0129         if any(J)
0130             tModel.b(J,1)=taskStructure(i).UBin(I)*-1;
0131             tModel.b(J,2)=taskStructure(i).LBin(I)*-1;
0132         end
0133     end
0134     %Set the outputs
0135     if ~isempty(taskStructure(i).outputs)
0136         [I, J]=ismember(upper(taskStructure(i).outputs),modelMets);
0137         J=J(I); %Only keep the ones with matches
0138         K=ismember(upper(taskStructure(i).outputs),'ALLMETS');
0139         L=~cellfun('isempty',strfind(upper(taskStructure(i).outputs),'ALLMETSIN'));
0140         %Check that all metabolites are either real metabolites or
0141         %ALLMETS/ALLMETSIN
0142         if ~all(I|K|L)
0143             fprintf(['ERROR: Could not find all outputs in "[' taskStructure(i).id '] ' taskStructure(i).description '"\n']);
0144             taskReport.ok(i,1)=false;
0145             tModel=model;
0146             continue;
0147         end
0148         if numel(J)~=numel(unique(J))
0149             EM=['The constraints on some output(s) in "[' taskStructure(i).id '] ' taskStructure(i).description '" are defined more than one time'];
0150             dispEM(EM);
0151         end
0152         %If all metabolites should be added
0153         if any(K)
0154             %Check if ALLMETS is the first metabolite. Otherwise print a
0155             %warning since it will write over any other constraints that
0156             %are set
0157             if K(1)==0
0158                 EM=['ALLMETS is used as an output in "[' taskStructure(i).id '] ' taskStructure(i).description '" but it it not the first metabolite in the list. Constraints defined for the metabolites before it will be over-written'];
0159                 dispEM(EM,false);
0160             end
0161             %Use the first match of ALLMETS. There should only be one, but
0162             %still..
0163             tModel.b(:,2)=taskStructure(i).UBout(find(K,1));
0164         end
0165         %If metabolites in a specific compartment should be used
0166         if any(L)
0167             L=find(L);
0168             for j=1:numel(L)
0169                 %The compartment defined
0170                 compartment=upper(taskStructure(i).outputs{L(j)}(11:end-1));
0171                 %Check if it exists in the model
0172                 C=find(ismember(upper(model.comps),compartment));
0173                 if any(C)
0174                     %Match to metabolites
0175                     tModel.b(model.metComps==C,2)=taskStructure(i).UBout(L(j));
0176                 else
0177                     EM=['The compartment defined for ALLMETSIN in "[' taskStructure(i).id '] ' taskStructure(i).description '" does not exist'];
0178                     dispEM(EM);
0179                 end
0180             end
0181         end
0182         %Then add the normal constraints
0183         if any(J)
0184             %Verify that IN and OUT bounds are consistent. Cannot require
0185             %that a metabolite is simultaneously input AND output at some
0186             %nonzero flux.
0187             I = find(I);  % otherwise indexing becomes confusing
0188             nonzero_LBin = tModel.b(J,2) < 0;
0189             nonzero_LBout = taskStructure(i).LBout(I) > 0;
0190             if any(nonzero_LBin & nonzero_LBout)
0191                 EM=['The IN LB and OUT LB in "[' taskStructure(i).id '] ' taskStructure(i).description '" cannot be nonzero for the same metabolite'];
0192                 dispEM(EM);
0193             end
0194             tModel.b(J(nonzero_LBout),1)=taskStructure(i).LBout(I(nonzero_LBout));
0195             tModel.b(J,2)=taskStructure(i).UBout(I);
0196         end
0197     end
0198     %Add new rxns
0199     if ~isempty(taskStructure(i).equations)
0200         rxn.equations=taskStructure(i).equations;
0201         rxn.lb=taskStructure(i).LBequ;
0202         rxn.ub=taskStructure(i).UBequ;
0203         rxn.rxns=strcat({'TEMPORARY_'},num2str((1:numel(taskStructure(i).equations))'));
0204         %Allow for new metabolites to be added. This is because it should
0205         %be possible to add, say, a whole new pathway
0206         tModel=addRxns(tModel,rxn,3,[],true);
0207     end
0208     %Add changed bounds
0209     if ~isempty(taskStructure(i).changed)
0210         tModel=setParam(tModel,'lb',taskStructure(i).changed,taskStructure(i).LBrxn);
0211         tModel=setParam(tModel,'ub',taskStructure(i).changed,taskStructure(i).UBrxn);
0212     end
0213     
0214     %Solve and print
0215     sol=solveLP(tModel);
0216     if ~isempty(sol.x)
0217         %assign the fluxes
0218         essentialFluxes(:,i) = sol.x(1:numel(model.rxns));
0219         
0220         if ~taskStructure(i).shouldFail
0221             taskReport.ok(i,1)=true;
0222             if printOnlyFailed==false && printOutput==true
0223                 fprintf(['PASS: [' taskStructure(i).id '] ' taskStructure(i).description '\n']);
0224             end
0225             %Calculate the essential reactions
0226             if getEssential==true
0227                 [~, taskEssential]=getEssentialRxns(tModel);
0228                 %This is because there could be more reactions in tModel
0229                 %than in model
0230                 essentialRxns(taskEssential(taskEssential<=numel(model.rxns)),i)=true;
0231             end
0232         else
0233             taskReport.ok(i,1)=false;
0234             if printOutput==true
0235                 fprintf(['PASS (should fail): [' taskStructure(i).id '] ' taskStructure(i).description '\n']);
0236             end
0237         end
0238     else
0239         if ~taskStructure(i).shouldFail
0240             taskReport.ok(i,1)=false;
0241             if printOutput==true
0242                 fprintf(['FAIL: [' taskStructure(i).id '] ' taskStructure(i).description '\n']);
0243             end
0244         else
0245             taskReport.ok(i,1)=true;
0246             if printOnlyFailed==false && printOutput==true
0247                 fprintf(['FAIL (should fail): [' taskStructure(i).id '] ' taskStructure(i).description '\n']);
0248             end
0249         end
0250     end
0251     if taskStructure(i).printFluxes && ~isempty(sol.x)
0252         sol=solveLP(tModel,1);
0253         if ~isempty(sol.x)
0254             printFluxes(tModel,sol.x,false,10^-6,[],'%rxnID (%eqn):%flux\n');
0255             fprintf('\n');
0256         end
0257     end
0258     tModel=model;
0259 end
0260 end