fitTasks

 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

   Input:
   model           model structure
   refModel        reference model from which to include reactions
   inputFile       a task list in Excel format. See the function
                   parseTaskList for details (optional if taskStructure is
                   supplied)
   printOutput     true if the results of the test should be displayed
                   (optional, 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 (optional, default is -1 for all reactions)
   taskStructure   structure with the tasks, as from parseTaskList. If
                   this is supplied then inputFile is ignored (optional)
   params          parameter structure as used by getMILPParams (optional)


   Output:
   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)

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