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