0001 function [outModel addedRxns]=fitTasks(model,refModel,inputFile,printOutput,rxnScores,taskStructure,params)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
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
0066 model.b=zeros(numel(model.mets),2);
0067 modelMets=upper(strcat(model.metNames,'[',model.comps(model.metComps),']'));
0068
0069
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
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
0092
0093 goodMets=I|K|L;
0094 if ~all(goodMets)
0095
0096
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
0102 met.metNames=refModel.metNames(metMatch);
0103 met.compartments=refModel.comps(refModel.metComps(metMatch));
0104
0105
0106
0107 model=addMets(model,met);
0108 tModel=addMets(tModel,met);
0109 modelMets=[modelMets;upper(taskStructure(i).inputs(~goodMets))];
0110 end
0111
0112
0113
0114
0115 I(~goodMets)=true;
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
0122 if any(K)
0123
0124
0125
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
0130
0131 tModel.b(:,1)=taskStructure(i).UBin(find(K,1))*-1;
0132 end
0133
0134 if any(L)
0135 L=find(L);
0136 for j=1:numel(L)
0137
0138 compartment=upper(taskStructure(i).inputs{L(j)}(11:end-1));
0139
0140 C=find(ismember(upper(model.comps),compartment));
0141 if any(C)
0142
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
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
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
0161
0162 goodMets=I|K|L;
0163 if ~all(goodMets)
0164
0165
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
0171 met.metNames=refModel.metNames(metMatch);
0172 met.compartments=refModel.comps(refModel.metComps(metMatch));
0173
0174
0175
0176 model=addMets(model,met);
0177 tModel=addMets(tModel,met);
0178 modelMets=[modelMets;upper(taskStructure(i).outputs(~goodMets))];
0179 end
0180
0181
0182
0183
0184 I(~goodMets)=true;
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
0191 if any(K)
0192
0193
0194
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
0199
0200 tModel.b(:,2)=taskStructure(i).UBout(find(K,1));
0201 end
0202
0203 if any(L)
0204 L=find(L);
0205 for j=1:numel(L)
0206
0207 compartment=upper(taskStructure(i).outputs{L(j)}(11:end-1));
0208
0209 C=find(ismember(upper(model.comps),compartment));
0210 if any(C)
0211
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
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
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
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
0240 sol=solveLP(tModel);
0241 if isempty(sol.x)
0242
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
0258
0259
0260 model=mergeModels({model,removeRxns(newModel,setdiff(newModel.rxns,newRxns),true,true)},true);
0261
0262
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
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
0284
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
0295
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