ftINITFillGapsForAllTasks

 ftINITFillGapsForAllTasks
   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
   This is similar to the old fitTasks function but optimized for ftINIT.

   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
   params          parameter structure as used by getMILPParams
   verbose         if true, the MILP progression will be shown. 


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