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)


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