Home > RAVEN > fitTasks.m

fitTasks

PURPOSE ^

fitTasks

SYNOPSIS ^

function [outModel addedRxns]=fitTasks(model,refModel,inputFile,printOutput,rxnScores,taskStructure,params)

DESCRIPTION ^

 fitTasks
   Fills gaps in a model by including reactions from a reference model,
   so that the resulting model can perform all the tasks in a task list

   model           model structure
   refModel        reference model from which to include reactions 
   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)
   rxnScores       scores for each of the reactions in the reference
                   model. Only negative scores are allowed. The solver will
                   try to maximize the sum of the scores for the included
                   reactions (opt, default is -1 for all reactions)
   taskStructure   structure with the tasks, as from parseTaskList. If
                   this is supplied then inputFile is ignored (opt)
   params          parameter structure as used by getMILPParams (opt)

   outModel        model structure with reactions added to perform the
                   tasks
   addedRxns       MxN matrix with the added reactions (M) from refModel 
                   for each task (N). An element is true if the corresponding
                   reaction is added in the corresponding task.
                   Failed tasks and SHOULD FAIL tasks are ignored

   This function fills gaps in a model by using a reference model, so
   that the resulting model can perform a list of metabolic tasks. The
   gap-filling is done in a task-by-task manner, rather than solving for
   all tasks at once. This means that the order of the tasks could influence
   the result.

   Usage: [outModel addedRxns]=fitTasks(model,refModel,inputFile,printOutput,...
           rxnScores,taskStructure,params)

   Rasmus Agren, 2013-04-22

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [outModel addedRxns]=fitTasks(model,refModel,inputFile,printOutput,rxnScores,taskStructure,params)
0002 % fitTasks
0003 %   Fills gaps in a model by including reactions from a reference model,
0004 %   so that the resulting model can perform all the tasks in a task list
0005 %
0006 %   model           model structure
0007 %   refModel        reference model from which to include reactions
0008 %   inputFile       a task list in Excel format. See the function
0009 %                   parseTaskList for details (opt if taskStructure is
0010 %                   supplied)
0011 %   printOutput     true if the results of the test should be displayed
0012 %                   (opt, default true)
0013 %   rxnScores       scores for each of the reactions in the reference
0014 %                   model. Only negative scores are allowed. The solver will
0015 %                   try to maximize the sum of the scores for the included
0016 %                   reactions (opt, default is -1 for all reactions)
0017 %   taskStructure   structure with the tasks, as from parseTaskList. If
0018 %                   this is supplied then inputFile is ignored (opt)
0019 %   params          parameter structure as used by getMILPParams (opt)
0020 %
0021 %   outModel        model structure with reactions added to perform the
0022 %                   tasks
0023 %   addedRxns       MxN matrix with the added reactions (M) from refModel
0024 %                   for each task (N). An element is true if the corresponding
0025 %                   reaction is added in the corresponding task.
0026 %                   Failed tasks and SHOULD FAIL tasks are ignored
0027 %
0028 %   This function fills gaps in a model by using a reference model, so
0029 %   that the resulting model can perform a list of metabolic tasks. The
0030 %   gap-filling is done in a task-by-task manner, rather than solving for
0031 %   all tasks at once. This means that the order of the tasks could influence
0032 %   the result.
0033 %
0034 %   Usage: [outModel addedRxns]=fitTasks(model,refModel,inputFile,printOutput,...
0035 %           rxnScores,taskStructure,params)
0036 %
0037 %   Rasmus Agren, 2013-04-22
0038 %
0039 
0040 if nargin<4
0041     printOutput=true;
0042 end
0043 if nargin<5
0044     rxnScores=ones(numel(refModel.rxns),1)*-1;
0045 end
0046 if isempty(rxnScores)
0047     rxnScores=ones(numel(refModel.rxns),1)*-1;
0048 end
0049 if nargin<6
0050     taskStructure=[];
0051 end
0052 if nargin<7
0053     params=[];
0054 end
0055 
0056 if strcmpi(model.id,refModel.id)
0057     fprintf('NOTE: The model and reference model have the same IDs. The ID for the reference model was set to "refModel" in order to keep track of the origin of reactions.\n'); 
0058     refModel.id='refModel';
0059 end
0060 
0061 if any(rxnScores>=0)
0062     throw(MException('','Only negative values are allowed in rxnScores'));
0063 end
0064 
0065 %Prepare the input models a little
0066 model.b=zeros(numel(model.mets),2);
0067 modelMets=upper(strcat(model.metNames,'[',model.comps(model.metComps),']'));
0068 %This is the mets in the reference model. Used if the tasks involve
0069 %metabolites that doesn't exist in the model
0070 largeModelMets=upper(strcat(refModel.metNames,'[',refModel.comps(refModel.metComps),']'));
0071 
0072 if ~isfield(model,'unconstrained')
0073     fprintf('WARNING: 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.\n'); 
0074 end
0075 
0076 if isempty(taskStructure)
0077    taskStructure=parseTaskList(inputFile); 
0078 end
0079 
0080 tModel=model;
0081 addedRxns=false(numel(refModel.rxns),numel(taskStructure));
0082 supressWarnings=false;
0083 nAdded=0;
0084 for i=1:numel(taskStructure)
0085     if ~taskStructure(i).shouldFail
0086         %Set the inputs
0087         if ~isempty(taskStructure(i).inputs)
0088             [I J]=ismember(upper(taskStructure(i).inputs),modelMets);
0089             K=ismember(upper(taskStructure(i).inputs),'ALLMETS');
0090             L=~cellfun('isempty',strfind(upper(taskStructure(i).inputs),'ALLMETSIN'));
0091             %Check that all metabolites are either real metabolites or
0092             %ALLMETS/ALLMETSIN
0093             goodMets=I|K|L;
0094             if ~all(goodMets)
0095                 %Not all of the inputs could be found in the small model. Check
0096                 %if they exist in the large model
0097                 [found metMatch]=ismember(upper(taskStructure(i).inputs(~goodMets)),largeModelMets);
0098                 if ~all(found)
0099                     throw(MException('',['Could not find all inputs in "[' taskStructure(i).id '] ' taskStructure(i).description '" in either model']));
0100                 else
0101                    %Otherwise add them to the model
0102                    met.metNames=refModel.metNames(metMatch);
0103                    met.compartments=refModel.comps(refModel.metComps(metMatch));
0104                    
0105                    %Add the metabolite both to the base model and the model
0106                    %used in the current task
0107                    model=addMets(model,met);
0108                    tModel=addMets(tModel,met);
0109                    modelMets=[modelMets;upper(taskStructure(i).inputs(~goodMets))];
0110                 end
0111                 
0112                 %By now the indexes might be getting a bit confusing, but
0113                 %this is to update the indexes of the "real" metabolites to
0114                 %point to the newly added ones
0115                 I(~goodMets)=true; %All the bad ones are fixed at this stage
0116                 J(~goodMets)=numel(modelMets)-numel(metMatch)+1:numel(modelMets);
0117             end
0118             if numel(J(I))~=numel(unique(J(I)))
0119                 throw(MException('',['The constraints on some input(s) in "[' taskStructure(i).id '] ' taskStructure(i).description '" are defined more than one time']));  
0120             end
0121             %If all metabolites should be added
0122             if any(K)
0123                 %Check if ALLMETS is the first metabolite. Otherwise print a
0124                 %warning since it will write over any other constraints that
0125                 %are set
0126                 if K(1)==0
0127                     fprintf(['WARNING: 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\n']);
0128                 end
0129                 %Use the first match of ALLMETS. There should only be one, but
0130                 %still..
0131                 tModel.b(:,1)=taskStructure(i).UBin(find(K,1))*-1;
0132             end
0133             %If metabolites in a specific compartment should be used
0134             if any(L)
0135                 L=find(L);
0136                 for j=1:numel(L)
0137                     %The compartment defined
0138                     compartment=upper(taskStructure(i).inputs{L(j)}(11:end-1));
0139                     %Check if it exists in the model
0140                     C=find(ismember(upper(model.comps),compartment));
0141                     if any(C)
0142                         %Match to metabolites
0143                         tModel.b(model.metComps==C,1)=taskStructure(i).UBin(L(j))*-1;
0144                     else
0145                         throw(MException('',['The compartment defined for ALLMETSIN in "[' taskStructure(i).id '] ' taskStructure(i).description '" does not exist']));  
0146                     end
0147                 end
0148             end
0149             %Then add the normal constraints
0150             if any(J(I))
0151                 tModel.b(J(I),1)=taskStructure(i).UBin(I)*-1;
0152                 tModel.b(J(I),2)=taskStructure(i).LBin(I)*-1;
0153             end
0154         end
0155         %Set the outputs
0156         if ~isempty(taskStructure(i).outputs)
0157             [I J]=ismember(upper(taskStructure(i).outputs),modelMets);
0158             K=ismember(upper(taskStructure(i).outputs),'ALLMETS');
0159             L=~cellfun('isempty',strfind(upper(taskStructure(i).outputs),'ALLMETSIN'));
0160             %Check that all metabolites are either real metabolites or
0161             %ALLMETS/ALLMETSIN
0162             goodMets=I|K|L;
0163             if ~all(goodMets)
0164                 %Not all of the outputs could be found in the small model. Check
0165                 %if they exist in the large model
0166                 [found metMatch]=ismember(upper(taskStructure(i).outputs(~goodMets)),largeModelMets);
0167                 if ~all(found)
0168                     throw(MException('',['Could not find all outputs in "[' taskStructure(i).id '] ' taskStructure(i).description '" in either model']));
0169                 else
0170                    %Otherwise add them to the model
0171                    met.metNames=refModel.metNames(metMatch);
0172                    met.compartments=refModel.comps(refModel.metComps(metMatch));
0173                    
0174                    %Add the metabolite both to the base model and the model
0175                    %used in the current task
0176                    model=addMets(model,met);
0177                    tModel=addMets(tModel,met);
0178                    modelMets=[modelMets;upper(taskStructure(i).outputs(~goodMets))];
0179                 end
0180                 
0181                 %By now the indexes might be getting a bit confusing, but
0182                 %this is to update the indexes of the "real" metabolites to
0183                 %point to the newly added ones
0184                 I(~goodMets)=true; %All the bad ones are fixed at this stage
0185                 J(~goodMets)=numel(modelMets)-numel(metMatch)+1:numel(modelMets);
0186             end
0187             if numel(J(I))~=numel(unique(J(I)))
0188                 throw(MException('',['The constraints on some output(s) in "[' taskStructure(i).id '] ' taskStructure(i).description '" are defined more than one time']));  
0189             end
0190             %If all metabolites should be added
0191             if any(K)
0192                 %Check if ALLMETS is the first metabolite. Otherwise print a
0193                 %warning since it will write over any other constraints that
0194                 %are set
0195                 if K(1)==0
0196                     fprintf(['WARNING: 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\n']);
0197                 end
0198                 %Use the first match of ALLMETS. There should only be one, but
0199                 %still..
0200                 tModel.b(:,2)=taskStructure(i).UBout(find(K,1));
0201             end
0202             %If metabolites in a specific compartment should be used
0203             if any(L)
0204                 L=find(L);
0205                 for j=1:numel(L)
0206                     %The compartment defined
0207                     compartment=upper(taskStructure(i).outputs{L(j)}(11:end-1));
0208                     %Check if it exists in the model
0209                     C=find(ismember(upper(model.comps),compartment));
0210                     if any(C)
0211                         %Match to metabolites
0212                         tModel.b(model.metComps==C,2)=taskStructure(i).UBout(L(j));
0213                     else
0214                         throw(MException('',['The compartment defined for ALLMETSIN in "[' taskStructure(i).id '] ' taskStructure(i).description '" does not exist']));  
0215                     end
0216                 end
0217             end
0218             %Then add the normal constraints
0219             if any(J(I))
0220                 tModel.b(J(I),1)=taskStructure(i).LBout(I);
0221                 tModel.b(J(I),2)=taskStructure(i).UBout(I);
0222             end
0223         end 
0224         
0225         %Add new rxns
0226         if ~isempty(taskStructure(i).equations)
0227             rxn.equations=taskStructure(i).equations;
0228             rxn.lb=taskStructure(i).LBequ;
0229             rxn.ub=taskStructure(i).UBequ;
0230             rxn.rxns=strcat({'TEMPORARY_'},num2str((1:numel(taskStructure(i).equations))'));
0231             tModel=addRxns(tModel,rxn,3);
0232         end
0233         %Add changed bounds
0234         if ~isempty(taskStructure(i).changed)
0235            tModel=setParam(tModel,'lb',taskStructure(i).changed,taskStructure(i).LBrxn);
0236            tModel=setParam(tModel,'ub',taskStructure(i).changed,taskStructure(i).UBrxn);
0237         end
0238 
0239         %Solve and print. Display a warning if the problem is not solveable
0240         sol=solveLP(tModel);
0241         if isempty(sol.x)
0242             %Only do gap-filling if it cannot be solved
0243             failed=false;
0244             try
0245                 [crap crap newRxns newModel exitFlag]=fillGaps(tModel,refModel,false,true,supressWarnings,rxnScores,params);
0246                 if exitFlag==-2
0247                     fprintf(['WARNING: "[' taskStructure(i).id '] ' taskStructure(i).description '" was aborted before reaching optimality. Consider increasing params.maxTime\n\n']);
0248                 end
0249             catch
0250                 fprintf(['WARNING: "[' taskStructure(i).id '] ' taskStructure(i).description '" could not be performed for any set of reactions\n\n']);
0251                 failed=true;
0252             end
0253             if failed==false
0254                 if ~isempty(newRxns)
0255                     nAdded=nAdded+numel(newRxns);
0256                     
0257                     %Add the reactions to the base model. It is not correct to use newModel
0258                     %directly, as it may contain reactions/constraints that are specific to
0259                     %this task
0260                     model=mergeModels({model,removeRxns(newModel,setdiff(newModel.rxns,newRxns),true,true)},true);
0261                     
0262                     %Keep track of the added reactions
0263                     addedRxns(ismember(refModel.rxns,newRxns),i)=true;
0264                 end
0265                 if printOutput==true
0266                     fprintf(['[' taskStructure(i).id '] ' taskStructure(i).description ': Added ' num2str(numel(newRxns)) ' reaction(s), ' num2str(nAdded) ' reactions added in total\n']);
0267                 end
0268             end
0269         else
0270             if printOutput==true
0271                 fprintf(['[' taskStructure(i).id '] ' taskStructure(i).description ': Added 0 reaction(s), ' num2str(nAdded) ' reactions added in total\n']);
0272             end
0273         end
0274         supressWarnings=true;
0275         
0276         %Print the output if chosen
0277         if taskStructure(i).printFluxes && printOutput
0278             if ~isempty(sol.x)
0279                 sol=solveLP(tModel,1);
0280                 printFluxes(tModel,sol.x,false,10^-5,[],'%rxnID (%eqn):%flux\n');
0281                 fprintf('\n');
0282             else
0283                 %If the problem wasn't solveable then the gap-filled model
0284                 %should be used
0285                 if failed==false
0286                     sol=solveLP(newModel,1);
0287                     printFluxes(newModel,sol.x,false,10^-5,[],'%rxnID (%eqn):%flux\n');
0288                     fprintf('\n');
0289                 end
0290             end
0291         end
0292         
0293         tModel=model;
0294         %Since new mets are added by merging the new reactions and not only
0295         %from the task sheet
0296         modelMets=upper(strcat(model.metNames,'[',model.comps(model.metComps),']'));
0297     else
0298         fprintf(['WARNING: "[' taskStructure(i).id '] ' taskStructure(i).description '" is set as SHOULD FAIL. Such tasks cannot be modelled using this approach and the task is therefore ignored\n\n']);
0299     end
0300 end
0301 outModel=model;
0302 end

Generated on Tue 23-Apr-2013 15:18:37 by m2html © 2005