0001 function newModel=addRxns(model,rxnsToAdd,eqnType,compartment,allowNewMets)
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
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077 if nargin<4
0078 compartment=[];
0079 end
0080 if nargin<5
0081 allowNewMets=false;
0082 end
0083
0084 newModel=model;
0085
0086
0087 if isempty(rxnsToAdd)
0088 return;
0089 end
0090
0091
0092 if ~isnumeric(eqnType)
0093 throw(MException('','eqnType must be numeric'));
0094 else
0095 if ~ismember(eqnType,[1 2 3])
0096 throw(MException('','eqnType must be 1, 2, or 3'));
0097 end
0098 end
0099
0100 if eqnType==2 || (eqnType==1 && allowNewMets==true)
0101 if ~ischar(compartment)
0102 throw(MException('','compartment must be a string'));
0103 end
0104 if ~ismember(compartment,model.comps)
0105 throw(MException('','compartment must match one of the compartments in model.comps'));
0106 end
0107 end
0108
0109 if ~isfield(rxnsToAdd,'rxns')
0110 throw(MException('','rxns is a required field in rxnsToAdd'));
0111 else
0112
0113 rxnsToAdd.rxns=rxnsToAdd.rxns(:);
0114 end
0115 if ~isfield(rxnsToAdd,'equations')
0116 throw(MException('','equations is a required field in rxnsToAdd'));
0117 end
0118
0119 if any(ismember(rxnsToAdd.rxns,model.rxns))
0120 throw(MException('','One or more reaction id was already present in the model. Use changeRxns or remove the existing reactions before adding new ones'));
0121 end
0122
0123 if ~iscellstr(rxnsToAdd.rxns) && ~ischar(rxnsToAdd.rxns)
0124
0125 throw(MException('','rxnsToAdd.rxns must be a cell array of strings'));
0126 else
0127 rxnsToAdd.rxns=cellstr(rxnsToAdd.rxns);
0128 end
0129 if ~iscellstr(rxnsToAdd.equations) && ~ischar(rxnsToAdd.equations)
0130
0131 throw(MException('','rxnsToAdd.equations must be a cell array of strings'));
0132 else
0133 rxnsToAdd.equations=cellstr(rxnsToAdd.equations);
0134 end
0135
0136
0137 illegalCells=regexp(rxnsToAdd.rxns,'[^a-z_A-Z0-9]', 'once');
0138 if ~isempty(cell2mat(illegalCells))
0139 errorText='Illegal character(s) in reaction IDs:\n';
0140 for i=1:length(illegalCells)
0141 if ~isempty(illegalCells{i})
0142 errorText=[errorText rxnsToAdd.rxns{i} '\n'];
0143 end
0144 end
0145 throw(MException('',errorText));
0146 end
0147 if isfield(rxnsToAdd,'rxnNames')
0148
0149 illegalCells=regexp(rxnsToAdd.rxnNames,'["%<>\\]', 'once');
0150 if ~isempty(cell2mat(illegalCells))
0151 errorText='Illegal character(s) in reaction names:\n';
0152 for i=1:length(illegalCells)
0153 if ~isempty(illegalCells{i})
0154 errorText=[errorText rxnsToAdd.rxnNames{i} '\n'];
0155 end
0156 end
0157 throw(MException('',errorText));
0158 end
0159 end
0160
0161 nRxns=numel(rxnsToAdd.rxns);
0162 nOldRxns=numel(model.rxns);
0163 filler=cell(nRxns,1);
0164 filler(:)={''};
0165 largeFiller=cell(nOldRxns,1);
0166 largeFiller(:)={''};
0167
0168
0169 if numel(rxnsToAdd.equations)~=nRxns
0170 throw(MException('','rxnsToAdd.equations must have the same number of elements as rxnsToAdd.rxns'));
0171 end
0172
0173
0174
0175 [S mets badRxns reversible]=constructS(rxnsToAdd.equations);
0176 if any(badRxns)
0177 fprintf('WARNING: The following equations have one or more metabolites both as substrate and product. Only the net equations will be added\n');
0178 I=find(badRxns);
0179 for i=1:numel(I)
0180 fprintf(['\t' rxnsToAdd.rxns{I(i)} '\n']);
0181 end
0182 end
0183 newModel.rev=[newModel.rev;reversible];
0184 newModel.rxns=[newModel.rxns;rxnsToAdd.rxns(:)];
0185
0186 if isfield(rxnsToAdd,'rxnNames')
0187 if numel(rxnsToAdd.rxnNames)~=nRxns
0188 throw(MException('','rxnsToAdd.rxnNames must have the same number of elements as rxnsToAdd.rxns'));
0189 end
0190
0191 if ~isfield(newModel,'rxnNames')
0192 newModel.rxnNames=largeFiller;
0193 end
0194 newModel.rxnNames=[newModel.rxnNames;rxnsToAdd.rxnNames(:)];
0195 else
0196
0197 if isfield(newModel,'rxnNames')
0198 newModel.rxnNames=[newModel.rxnNames;filler];
0199 end
0200 end
0201
0202 if isfield(rxnsToAdd,'lb')
0203 if numel(rxnsToAdd.lb)~=nRxns
0204 throw(MException('','rxnsToAdd.lb must have the same number of elements as rxnsToAdd.rxns'));
0205 end
0206
0207 if ~isfield(newModel,'lb')
0208 newModel.lb=zeros(nOldRxns,1);
0209 newModel.lb(newModel.rev~=0)=-inf;
0210 end
0211 newModel.lb=[newModel.lb;rxnsToAdd.lb(:)];
0212 else
0213
0214 if isfield(newModel,'lb')
0215 I=zeros(nRxns,1);
0216 I(reversible~=0)=-inf;
0217 newModel.lb=[newModel.lb;I];
0218 end
0219 end
0220
0221 if isfield(rxnsToAdd,'ub')
0222 if numel(rxnsToAdd.ub)~=nRxns
0223 throw(MException('','rxnsToAdd.ub must have the same number of elements as rxnsToAdd.rxns'));
0224 end
0225
0226 if ~isfield(newModel,'ub')
0227 newModel.ub=inf(nOldRxns,1);
0228 end
0229 newModel.ub=[newModel.ub;rxnsToAdd.ub(:)];
0230 else
0231
0232 if isfield(newModel,'ub')
0233 newModel.ub=[newModel.ub;inf(nRxns,1)];
0234 end
0235 end
0236
0237 if isfield(rxnsToAdd,'c')
0238 if numel(rxnsToAdd.c)~=nRxns
0239 throw(MException('','rxnsToAdd.c must have the same number of elements as rxnsToAdd.rxns'));
0240 end
0241
0242 if ~isfield(newModel,'c')
0243 newModel.c=zeros(nOldRxns,1);
0244 end
0245 newModel.c=[newModel.c;rxnsToAdd.c(:)];
0246 else
0247
0248 if isfield(newModel,'c')
0249 newModel.c=[newModel.c;zeros(nRxns,1)];
0250 end
0251 end
0252
0253 if isfield(rxnsToAdd,'eccodes')
0254 if numel(rxnsToAdd.eccodes)~=nRxns
0255 throw(MException('','rxnsToAdd.eccodes must have the same number of elements as rxnsToAdd.rxns'));
0256 end
0257
0258 if ~isfield(newModel,'eccodes')
0259 newModel.eccodes=largeFiller;
0260 end
0261 newModel.eccodes=[newModel.eccodes;rxnsToAdd.eccodes(:)];
0262 else
0263
0264 if isfield(newModel,'eccodes')
0265 newModel.eccodes=[newModel.eccodes;filler];
0266 end
0267 end
0268
0269 if isfield(rxnsToAdd,'subSystems')
0270 if numel(rxnsToAdd.subSystems)~=nRxns
0271 throw(MException('','rxnsToAdd.subSystems must have the same number of elements as rxnsToAdd.rxns'));
0272 end
0273
0274 if ~isfield(newModel,'subSystems')
0275 newModel.subSystems=largeFiller;
0276 end
0277 newModel.subSystems=[newModel.subSystems;rxnsToAdd.subSystems(:)];
0278 else
0279
0280 if isfield(newModel,'subSystems')
0281 newModel.subSystems=[newModel.subSystems;filler];
0282 end
0283 end
0284 if isfield(rxnsToAdd,'rxnMiriams')
0285 if numel(rxnsToAdd.rxnMiriams)~=nRxns
0286 throw(MException('','rxnsToAdd.rxnMiriams must have the same number of elements as rxnsToAdd.rxns'));
0287 end
0288
0289 if ~isfield(newModel,'rxnMiriams')
0290 newModel.rxnMiriams=cell(nOldRxns,1);
0291 end
0292 newModel.rxnMiriams=[newModel.rxnMiriams;rxnsToAdd.rxnMiriams(:)];
0293 else
0294
0295 if isfield(newModel,'rxnMiriams')
0296 newModel.rxnMiriams=[newModel.rxnMiriams;cell(nRxns,1)];
0297 end
0298 end
0299 if isfield(rxnsToAdd,'grRules')
0300 if numel(rxnsToAdd.grRules)~=nRxns
0301 throw(MException('','rxnsToAdd.grRules must have the same number of elements as rxnsToAdd.rxns'));
0302 end
0303
0304 if ~isfield(newModel,'grRules')
0305 newModel.grRules=largeFiller;
0306 end
0307 newModel.grRules=[newModel.grRules;rxnsToAdd.grRules(:)];
0308 else
0309
0310 if isfield(newModel,'grRules')
0311 newModel.grRules=[newModel.grRules;filler];
0312 end
0313 end
0314
0315 if isfield(rxnsToAdd,'rxnFrom')
0316 if numel(rxnsToAdd.rxnFrom)~=nRxns
0317 throw(MException('','rxnsToAdd.rxnFrom must have the same number of elements as rxnsToAdd.rxns'));
0318 end
0319
0320 if ~isfield(newModel,'rxnFrom')
0321 newModel.rxnFrom=largeFiller;
0322 end
0323 newModel.rxnFrom=[newModel.rxnFrom;rxnsToAdd.rxnFrom(:)];
0324 else
0325
0326 if isfield(newModel,'rxnFrom')
0327 newModel.rxnFrom=[newModel.rxnFrom;filler];
0328 end
0329 end
0330
0331
0332
0333 if eqnType==1
0334 illegalCells=regexp(mets,'[^a-z_A-Z0-9]', 'once');
0335 if ~isempty(cell2mat(illegalCells))
0336 errorText='Illegal character(s) in metabolite IDs:\n';
0337 for i=1:length(illegalCells)
0338 if ~isempty(illegalCells{i})
0339 errorText=[errorText mets{i} '\n'];
0340 end
0341 end
0342 throw(MException('',errorText));
0343 end
0344 else
0345
0346 illegalCells=regexp(mets,'["%<>\\]', 'once');
0347 if ~isempty(cell2mat(illegalCells))
0348 errorText='Illegal character(s) in metabolite names:\n';
0349 for i=1:length(illegalCells)
0350 if ~isempty(illegalCells{i})
0351 errorText=[errorText mets{i} '\n'];
0352 end
0353 end
0354 throw(MException('',errorText));
0355 end
0356 end
0357
0358
0359
0360 if eqnType==1
0361 [I J]=ismember(mets,model.mets);
0362 if ~all(I)
0363 if allowNewMets==true
0364
0365 metsToAdd.mets=mets(~I);
0366 metsToAdd.metNames=metsToAdd.mets;
0367 metsToAdd.compartments=compartment;
0368 newModel=addMets(newModel,metsToAdd);
0369 else
0370 throw(MException('','One or more equations contain metabolites that are not in model.mets. Set allowNewMets to true to allow this function to add metabolites or use addMets to add them before calling this function'));
0371 end
0372 end
0373
0374 metIndexes=J;
0375 metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0376 end
0377
0378
0379 if eqnType==2 || eqnType==3
0380
0381 t2=strcat(model.metNames,'¤¤¤',model.comps(model.metComps));
0382 end
0383
0384
0385 if eqnType==2
0386
0387
0388 t1=strcat(mets,'¤¤¤',compartment);
0389 [I J]=ismember(t1,t2);
0390
0391 if ~all(I)
0392 if allowNewMets==true
0393
0394 metsToAdd.metNames=mets(~I);
0395 metsToAdd.compartments=compartment;
0396 newModel=addMets(newModel,metsToAdd);
0397 else
0398 throw(MException('','One or more equations contain metabolites that are not in model.mets. Set allowNewMets to true to allow this function to add metabolites or use addMets to add them before calling this function'));
0399 end
0400 end
0401
0402
0403 metIndexes=J;
0404 metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0405 end
0406
0407
0408 if eqnType==3
0409
0410 metNames=cell(numel(mets),1);
0411 compartments=metNames;
0412 for i=1:numel(mets)
0413 starts=max(strfind(mets{i},'['));
0414 ends=max(strfind(mets{i},']'));
0415
0416
0417 if isempty(starts) || isempty(ends) || ends<numel(mets{i})
0418 throw(MException('',['The metabolite ' mets{i} ' is not correctly formatted for eqnType=3']));
0419 end
0420
0421
0422 compartments{i}=mets{i}(starts+1:ends-1);
0423 I=ismember(compartments(i),newModel.comps);
0424 if ~I
0425 throw(MException('',['The metabolite ' mets{i} ' has a compartment that is not in model.comps']));
0426 end
0427 metNames{i}=mets{i}(1:starts-1);
0428 end
0429
0430
0431 t1=strcat(metNames,'¤¤¤',compartments);
0432 [I J]=ismember(t1,t2);
0433
0434 if ~all(I)
0435 if allowNewMets==true
0436
0437 metsToAdd.metNames=metNames(~I);
0438 metsToAdd.compartments=compartments(~I);
0439 newModel=addMets(newModel,metsToAdd);
0440 else
0441 throw(MException('','One or more equations contain metabolites that are not in model.metNames. Set allowNewMets to true to allow this function to add metabolites or use addMets to add them before calling this function'));
0442 end
0443 end
0444
0445
0446 metIndexes=J;
0447 metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0448 end
0449
0450
0451
0452 newModel.S=[newModel.S sparse(size(newModel.S,1),nRxns)];
0453 for i=1:nRxns
0454 newModel.S(metIndexes,nOldRxns+i)=S(:,i);
0455
0456
0457 if isfield(newModel,'grRules')
0458 rule=newModel.grRules{nOldRxns+i};
0459 rule=strrep(rule,'(','');
0460 rule=strrep(rule,')','');
0461 rule=strrep(rule,' or ',' ');
0462 rule=strrep(rule,' and ',' ');
0463 genes=regexp(rule,' ','split');
0464 [I J]=ismember(genes,newModel.genes);
0465 if ~all(I) && any(rule)
0466 throw(MException('',['Not all genes for reaction ' rxnsToAdd.rxns{i} ' were found in model.genes. If needed, add genes with addGenes before calling this function']));
0467 end
0468 if any(rule)
0469 newModel.rxnGeneMat(nOldRxns+i,J)=1;
0470 else
0471
0472
0473 newModel.rxnGeneMat=[newModel.rxnGeneMat;sparse(1,numel(newModel.genes))];
0474 end
0475 end
0476 end
0477 end