Home > RAVEN > checkTasks.m

checkTasks

PURPOSE ^

checkTasks

SYNOPSIS ^

function [taskReport essentialRxns taskStructure]=checkTasks(model,inputFile,printOutput,printOnlyFailed,getEssential,taskStructure)

DESCRIPTION ^

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

   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)

   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

   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)

   Rasmus Agren, 2013-11-17

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Mon 06-Jan-2014 14:58:12 by m2html © 2005