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)

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