0001 function newModel=addRxns(model,rxnsToAdd,eqnType,compartment,allowNewMets,allowNewGenes)
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
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104 if nargin<3
0105 eqnType=1;
0106 elseif ~isnumeric(eqnType)
0107 EM='eqnType must be numeric';
0108 dispEM(EM);
0109 elseif ~ismember(eqnType,[1 2 3])
0110 EM='eqnType must be 1, 2, or 3';
0111 dispEM(EM);
0112 end
0113
0114 if nargin<4
0115 compartment=[];
0116 else
0117 compartment=char(compartment);
0118 end
0119 if nargin<5
0120 allowNewMets=false;
0121 elseif ~islogical(allowNewMets)
0122 allowNewMets=char(allowNewMets);
0123 end
0124 if nargin<6
0125 allowNewGenes=false;
0126 end
0127
0128 if allowNewGenes & isfield(rxnsToAdd,'grRules')
0129 genesToAdd.genes = strjoin(convertCharArray(rxnsToAdd.grRules));
0130 genesToAdd.genes = regexp(genesToAdd.genes,' |)|(|and|or','split');
0131 genesToAdd.genes = genesToAdd.genes(~cellfun(@isempty,genesToAdd.genes));
0132 genesToAdd.genes = setdiff(unique(genesToAdd.genes),model.genes);
0133 if isfield(model,'geneComps')
0134 genesToAdd.geneComps(1:numel(genesToAdd.genes)) = repmat(11,numel(genesToAdd.genes),1);
0135 end
0136 if ~isempty(genesToAdd.genes)
0137 fprintf('\nNew genes added to the model:\n')
0138 fprintf([strjoin(genesToAdd.genes,'\n') '\n'])
0139 newModel=addGenesRaven(model,genesToAdd);
0140 else
0141 newModel=model;
0142 end
0143 else
0144 newModel=model;
0145 end
0146
0147
0148 if isempty(rxnsToAdd)
0149 return;
0150 end
0151
0152 if eqnType==2 || (eqnType==1 && allowNewMets==true)
0153 compartment=char(compartment);
0154 if ~ismember(compartment,model.comps)
0155 EM='compartment must match one of the compartments in model.comps';
0156 dispEM(EM);
0157 end
0158 end
0159
0160 if ~isfield(rxnsToAdd,'rxns')
0161 EM='rxns is a required field in rxnsToAdd';
0162 dispEM(EM);
0163 else
0164 rxnsToAdd.rxns=convertCharArray(rxnsToAdd.rxns);
0165
0166 rxnsToAdd.rxns=rxnsToAdd.rxns(:);
0167 end
0168 if ~(isfield(rxnsToAdd,'equations') || (isfield(rxnsToAdd,'mets') && isfield(rxnsToAdd,'stoichCoeffs')))
0169 EM='Either "equations" or "mets"+"stoichCoeffs" are required fields in rxnsToAdd';
0170 dispEM(EM);
0171 end
0172
0173 if any(ismember(rxnsToAdd.rxns,model.rxns))
0174 EM='One or more reaction id was already present in the model. Use changeRxns or remove the existing reactions before adding new ones';
0175 dispEM(EM);
0176 end
0177
0178
0179 if isfield(rxnsToAdd,'equations')
0180 rxnsToAdd.equations=convertCharArray(rxnsToAdd.equations);
0181
0182
0183 else
0184
0185
0186 if iscellstr(rxnsToAdd.mets)
0187 rxnsToAdd.mets={rxnsToAdd.mets};
0188 end
0189 if isnumeric(rxnsToAdd.stoichCoeffs)
0190 rxnsToAdd.stoichCoeffs = {rxnsToAdd.stoichCoeffs};
0191 end
0192
0193
0194 if ~iscell(rxnsToAdd.mets) || ~iscell(rxnsToAdd.stoichCoeffs)
0195 EM='rxnsToAdd.mets & rxnsToAdd.stoichCoeffs must be cell arrays';
0196 dispEM(EM);
0197 elseif length(rxnsToAdd.stoichCoeffs) ~= length(rxnsToAdd.mets)
0198 EM = 'rxnsToAdd.stoichCoeffs must have the same number of elements as rxnsToAdd.mets';
0199 dispEM(EM);
0200 end
0201
0202 if ~isfield(rxnsToAdd,'lb')
0203
0204 rxnsToAdd.lb=-inf(size(rxnsToAdd.mets));
0205 elseif ~isnumeric(rxnsToAdd.lb)
0206 EM = 'rxnsToAdd.lb must be a vector';
0207 dispEM(EM);
0208 elseif length(rxnsToAdd.lb) ~= length(rxnsToAdd.mets)
0209 EM = 'rxnsToAdd.lb must have the same number of elements as rxnsToAdd.mets';
0210 dispEM(EM);
0211 end
0212
0213 rxnsToAdd.equations = cell(size(rxnsToAdd.mets));
0214 for i = 1:length(rxnsToAdd.mets)
0215 mets = rxnsToAdd.mets{i};
0216 stoichCoeffs = rxnsToAdd.stoichCoeffs{i};
0217 if isfield(rxnsToAdd,'ub')
0218 isrev = rxnsToAdd.lb(i) < 0 && rxnsToAdd.ub(i) > 0;
0219 else
0220 isrev = rxnsToAdd.lb(i) < 0;
0221 end
0222 rxnsToAdd.equations{i} = buildEquation(mets,stoichCoeffs,isrev);
0223 end
0224 end
0225
0226 nRxns=numel(rxnsToAdd.rxns);
0227 nOldRxns=numel(model.rxns);
0228 filler=cell(nRxns,1);
0229 filler(:)={''};
0230 cellfiller=cellfun(@(x) cell(0,0),filler,'UniformOutput',false);
0231 largeFiller=cell(nOldRxns,1);
0232 largeFiller(:)={''};
0233 celllargefiller=cellfun(@(x) cell(0,0),largeFiller,'UniformOutput',false);
0234
0235
0236 if numel(rxnsToAdd.equations)~=nRxns
0237 EM='rxnsToAdd.equations must have the same number of elements as rxnsToAdd.rxns';
0238 dispEM(EM);
0239 end
0240
0241
0242
0243 [S, mets, badRxns, reversible]=constructS(rxnsToAdd.equations);
0244 EM='The following equations have one or more metabolites both as substrate and product. Only the net equations will be added:';
0245 dispEM(EM,false,rxnsToAdd.rxns(badRxns));
0246
0247 newModel.rev=[newModel.rev;reversible];
0248 newModel.rxns=[newModel.rxns;rxnsToAdd.rxns(:)];
0249
0250 if isfield(rxnsToAdd,'rxnNames')
0251 rxnsToAdd.rxnNames=convertCharArray(rxnsToAdd.rxnNames);
0252 if numel(rxnsToAdd.rxnNames)~=nRxns
0253 EM='rxnsToAdd.rxnNames must have the same number of elements as rxnsToAdd.rxns';
0254 dispEM(EM);
0255 end
0256
0257 if ~isfield(newModel,'rxnNames')
0258 newModel.rxnNames=largeFiller;
0259 end
0260 newModel.rxnNames=[newModel.rxnNames;rxnsToAdd.rxnNames(:)];
0261 else
0262
0263 if isfield(newModel,'rxnNames')
0264 newModel.rxnNames=[newModel.rxnNames;filler];
0265 end
0266 end
0267
0268 if isfield(newModel,'annotation') & isfield(newModel.annotation,'defaultLB')
0269 newLb=newModel.annotation.defaultLB;
0270 else
0271 newLb=-inf;
0272 end
0273
0274 if isfield(rxnsToAdd,'lb')
0275 if numel(rxnsToAdd.lb)~=nRxns
0276 EM='rxnsToAdd.lb must have the same number of elements as rxnsToAdd.rxns';
0277 dispEM(EM);
0278 end
0279
0280 if ~isfield(newModel,'lb')
0281 newModel.lb=zeros(nOldRxns,1);
0282 newModel.lb(newModel.rev~=0)=newLb;
0283 end
0284 newModel.lb=[newModel.lb;rxnsToAdd.lb(:)];
0285 else
0286
0287 if isfield(newModel,'lb')
0288 I=zeros(nRxns,1);
0289 I(reversible~=0)=newLb;
0290 newModel.lb=[newModel.lb;I];
0291 end
0292 end
0293
0294 if isfield(newModel,'annotation') & isfield(newModel.annotation,'defaultUB')
0295 newUb=newModel.annotation.defaultUB;
0296 else
0297 newUb=inf;
0298 end
0299
0300 if isfield(rxnsToAdd,'ub')
0301 if numel(rxnsToAdd.ub)~=nRxns
0302 EM='rxnsToAdd.ub must have the same number of elements as rxnsToAdd.rxns';
0303 dispEM(EM);
0304 end
0305
0306 if ~isfield(newModel,'ub')
0307 newModel.ub=repmat(newUb,nOldRxns,1);
0308 end
0309 newModel.ub=[newModel.ub;rxnsToAdd.ub(:)];
0310 else
0311
0312 if isfield(newModel,'ub')
0313 newModel.ub=[newModel.ub;repmat(newUb,nRxns,1)];
0314 end
0315 end
0316
0317 if isfield(rxnsToAdd,'c')
0318 if numel(rxnsToAdd.c)~=nRxns
0319 EM='rxnsToAdd.c must have the same number of elements as rxnsToAdd.rxns';
0320 dispEM(EM);
0321 end
0322
0323 if ~isfield(newModel,'c')
0324 newModel.c=zeros(nOldRxns,1);
0325 end
0326 newModel.c=[newModel.c;rxnsToAdd.c(:)];
0327 else
0328
0329 if isfield(newModel,'c')
0330 newModel.c=[newModel.c;zeros(nRxns,1)];
0331 end
0332 end
0333
0334 if isfield(rxnsToAdd,'eccodes')
0335 rxnsToAdd.eccodes=convertCharArray(rxnsToAdd.eccodes);
0336 if numel(rxnsToAdd.eccodes)~=nRxns
0337 EM='rxnsToAdd.eccodes must have the same number of elements as rxnsToAdd.rxns';
0338 dispEM(EM);
0339 end
0340
0341 if ~isfield(newModel,'eccodes')
0342 newModel.eccodes=largeFiller;
0343 end
0344 newModel.eccodes=[newModel.eccodes;rxnsToAdd.eccodes(:)];
0345 else
0346
0347 if isfield(newModel,'eccodes')
0348 newModel.eccodes=[newModel.eccodes;filler];
0349 end
0350 end
0351
0352 if isfield(rxnsToAdd,'subSystems')
0353 if numel(rxnsToAdd.subSystems)~=nRxns
0354 EM='rxnsToAdd.subSystems must have the same number of elements as rxnsToAdd.rxns';
0355 dispEM(EM);
0356 end
0357 for i=1:numel(rxnsToAdd.subSystems)
0358 if ischar(rxnsToAdd.subSystems{i})
0359 rxnsToAdd.subSystems{i}=rxnsToAdd.subSystems(i);
0360 end
0361 end
0362
0363 if ~isfield(newModel,'subSystems')
0364 newModel.subSystems=celllargefiller;
0365 end
0366 newModel.subSystems=[newModel.subSystems;rxnsToAdd.subSystems(:)];
0367 else
0368
0369 if isfield(newModel,'subSystems')
0370 newModel.subSystems=[newModel.subSystems;cellfiller];
0371 end
0372 end
0373 if isfield(rxnsToAdd,'rxnMiriams')
0374 if numel(rxnsToAdd.rxnMiriams)~=nRxns
0375 EM='rxnsToAdd.rxnMiriams must have the same number of elements as rxnsToAdd.rxns';
0376 dispEM(EM);
0377 end
0378
0379 if ~isfield(newModel,'rxnMiriams')
0380 newModel.rxnMiriams=cell(nOldRxns,1);
0381 end
0382 newModel.rxnMiriams=[newModel.rxnMiriams;rxnsToAdd.rxnMiriams(:)];
0383 else
0384
0385 if isfield(newModel,'rxnMiriams')
0386 newModel.rxnMiriams=[newModel.rxnMiriams;cell(nRxns,1)];
0387 end
0388 end
0389 if isfield(rxnsToAdd,'rxnComps')
0390 if numel(rxnsToAdd.rxnComps)~=nRxns
0391 EM='rxnsToAdd.rxnComps must have the same number of elements as rxnsToAdd.rxns';
0392 dispEM(EM);
0393 end
0394
0395 if ~isfield(newModel,'rxnComps')
0396 newModel.rxnComps=ones(nOldRxns,1);
0397 EM='Adding reactions with compartment information to a model without such information. All existing reactions will be assigned to the first compartment';
0398 dispEM(EM,false);
0399 end
0400 newModel.rxnComps=[newModel.rxnComps;rxnsToAdd.rxnComps(:)];
0401 else
0402
0403 if isfield(newModel,'rxnComps')
0404 newModel.rxnComps=[newModel.rxnComps;ones(nRxns,1)];
0405
0406
0407 end
0408 end
0409 if isfield(rxnsToAdd,'grRules')
0410 rxnsToAdd.grRules=convertCharArray(rxnsToAdd.grRules);
0411 if numel(rxnsToAdd.grRules)~=nRxns
0412 EM='rxnsToAdd.grRules must have the same number of elements as rxnsToAdd.rxns';
0413 dispEM(EM);
0414 end
0415
0416 if ~isfield(newModel,'grRules')
0417 newModel.grRules=largeFiller;
0418 end
0419 newModel.grRules=[newModel.grRules;rxnsToAdd.grRules(:)];
0420 else
0421
0422 if isfield(newModel,'grRules')
0423 newModel.grRules=[newModel.grRules;filler];
0424 end
0425 end
0426
0427 if isfield(rxnsToAdd,'rxnFrom')
0428 rxnsToAdd.rxnFrom=convertCharArray(rxnsToAdd.rxnFrom);
0429 if numel(rxnsToAdd.rxnFrom)~=nRxns
0430 EM='rxnsToAdd.rxnFrom must have the same number of elements as rxnsToAdd.rxns';
0431 dispEM(EM);
0432 end
0433
0434 if ~isfield(newModel,'rxnFrom')
0435 newModel.rxnFrom=largeFiller;
0436 end
0437 newModel.rxnFrom=[newModel.rxnFrom;rxnsToAdd.rxnFrom(:)];
0438 else
0439
0440 if isfield(newModel,'rxnFrom')
0441 newModel.rxnFrom=[newModel.rxnFrom;filler];
0442 end
0443 end
0444
0445 if isfield(rxnsToAdd,'rxnNotes')
0446 rxnsToAdd.rxnNotes=convertCharArray(rxnsToAdd.rxnNotes);
0447 if numel(rxnsToAdd.rxnNotes)~=nRxns
0448 EM='rxnsToAdd.rxnNotes must have the same number of elements as rxnsToAdd.rxns';
0449 dispEM(EM);
0450 end
0451
0452 if ~isfield(newModel,'rxnNotes')
0453 newModel.rxnNotes=largeFiller;
0454 end
0455 newModel.rxnNotes=[newModel.rxnNotes;rxnsToAdd.rxnNotes(:)];
0456 else
0457
0458 if isfield(newModel,'rxnNotes')
0459 newModel.rxnNotes=[newModel.rxnNotes;filler];
0460 end
0461 end
0462
0463 if isfield(rxnsToAdd,'rxnReferences')
0464 rxnsToAdd.rxnReferences=convertCharArray(rxnsToAdd.rxnReferences);
0465 if numel(rxnsToAdd.rxnReferences)~=nRxns
0466 EM='rxnsToAdd.rxnReferences must have the same number of elements as rxnsToAdd.rxns';
0467 dispEM(EM);
0468 end
0469
0470 if ~isfield(newModel,'rxnReferences')
0471 newModel.rxnReferences=largeFiller;
0472 end
0473 newModel.rxnReferences=[newModel.rxnReferences;rxnsToAdd.rxnReferences(:)];
0474 else
0475
0476 if isfield(newModel,'rxnReferences')
0477 newModel.rxnReferences=[newModel.rxnReferences;filler];
0478 end
0479 end
0480
0481 if isfield(rxnsToAdd,'pwys')
0482 rxnsToAdd.pwys=convertCharArray(rxnsToAdd.pwys);
0483 if numel(rxnsToAdd.pwys)~=nRxns
0484 EM='rxnsToAdd.pwys must have the same number of elements as rxnsToAdd.rxns';
0485 dispEM(EM);
0486 end
0487
0488 if ~isfield(newModel,'pwys')
0489 newModel.pwys=largeFiller;
0490 end
0491 newModel.pwys=[newModel.pwys;rxnsToAdd.pwys(:)];
0492 else
0493
0494 if isfield(newModel,'pwys')
0495 newModel.pwys=[newModel.pwys;filler];
0496 end
0497 end
0498
0499 if isfield(rxnsToAdd,'rxnConfidenceScores')
0500 if numel(rxnsToAdd.rxnConfidenceScores)~=nRxns
0501 EM='rxnsToAdd.rxnConfidenceScores must have the same number of elements as rxnsToAdd.rxns';
0502 dispEM(EM);
0503 end
0504
0505 if ~isfield(newModel,'rxnConfidenceScores')
0506 newModel.rxnConfidenceScores=NaN(nOldRxns,1);
0507 end
0508 newModel.rxnConfidenceScores=[newModel.rxnConfidenceScores;rxnsToAdd.rxnConfidenceScores(:)];
0509 else
0510
0511 if isfield(newModel,'rxnConfidenceScores')
0512 newModel.rxnConfidenceScores=[newModel.rxnConfidenceScores;NaN(nRxns,1)];
0513 end
0514 end
0515
0516 if isfield(rxnsToAdd,'rxnDeltaG')
0517 if numel(rxnsToAdd.rxnDeltaG)~=nRxns
0518 EM='rxnsToAdd.rxnDeltaG must have the same number of elements as rxnsToAdd.rxns';
0519 dispEM(EM);
0520 end
0521
0522 if ~isfield(newModel,'rxnDeltaG')
0523 newModel.rxnDeltaG=NaN(nOldRxns,1);
0524 end
0525 newModel.rxnDeltaG=[newModel.rxnDeltaG;rxnsToAdd.rxnDeltaG(:)];
0526 else
0527
0528 if isfield(newModel,'rxnDeltaG')
0529 newModel.rxnDeltaG=[newModel.rxnDeltaG;NaN(nRxns,1)];
0530 end
0531 end
0532
0533
0534
0535
0536 if eqnType==1
0537 [I, J]=ismember(mets,model.mets);
0538 if ~all(I)
0539 if allowNewMets==true || ischar(allowNewMets)
0540
0541 metsToAdd.mets=mets(~I);
0542 metsToAdd.metNames=metsToAdd.mets;
0543 metsToAdd.compartments=compartment;
0544 if ischar(allowNewMets)
0545 newModel=addMets(newModel,metsToAdd,true,allowNewMets);
0546 else
0547 newModel=addMets(newModel,metsToAdd,true);
0548 end
0549 else
0550 EM='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. Are you sure that eqnType=1?';
0551 dispEM(EM);
0552 end
0553 end
0554
0555 metIndexes=J;
0556 metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0557 end
0558
0559
0560 if eqnType==2 || eqnType==3
0561
0562 t2=strcat(model.metNames,'***',model.comps(model.metComps));
0563 end
0564
0565
0566 if eqnType==2
0567
0568
0569
0570 t1=strcat(mets,'***',compartment);
0571 [I, J]=ismember(t1,t2);
0572
0573 if ~all(I)
0574 if allowNewMets==true || ischar(allowNewMets)
0575
0576 metsToAdd.metNames=mets(~I);
0577 metsToAdd.compartments=compartment;
0578 if ischar(allowNewMets)
0579 newModel=addMets(newModel,metsToAdd,true,allowNewMets);
0580 else
0581 newModel=addMets(newModel,metsToAdd,true);
0582 end
0583 else
0584 EM='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';
0585 dispEM(EM);
0586 end
0587 end
0588
0589
0590 metIndexes=J;
0591 metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0592 end
0593
0594
0595 if eqnType==3
0596
0597 metNames=cell(numel(mets),1);
0598 compartments=metNames;
0599 for i=1:numel(mets)
0600 starts=max(strfind(mets{i},'['));
0601 ends=max(strfind(mets{i},']'));
0602
0603
0604 if isempty(starts) || isempty(ends) || ends<numel(mets{i})
0605 EM=['The metabolite ' mets{i} ' is not correctly formatted for eqnType=3'];
0606 dispEM(EM);
0607 end
0608
0609
0610 compartments{i}=mets{i}(starts+1:ends-1);
0611 I=ismember(compartments(i),newModel.comps);
0612 if ~I
0613 EM=['The metabolite ' mets{i} ' has a compartment that is not in model.comps'];
0614 dispEM(EM);
0615 end
0616 metNames{i}=mets{i}(1:starts-1);
0617 end
0618
0619
0620 t1=strcat(metNames,'***',compartments);
0621 [I, J]=ismember(t1,t2);
0622
0623 if ~all(I)
0624 if allowNewMets==true | ischar(allowNewMets)
0625
0626 metsToAdd.metNames=metNames(~I);
0627 metsToAdd.compartments=compartments(~I);
0628 if ischar(allowNewMets)
0629 newModel=addMets(newModel,metsToAdd,true,allowNewMets);
0630 else
0631 newModel=addMets(newModel,metsToAdd,true);
0632 end
0633 else
0634 EM='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';
0635 dispEM(EM);
0636 end
0637 end
0638
0639
0640 metIndexes=J;
0641 metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0642 end
0643
0644
0645 newModel.S=[newModel.S sparse(size(newModel.S,1),nRxns)];
0646
0647 for i=1:nRxns
0648 newModel.S(metIndexes,nOldRxns+i)=S(:,i);
0649
0650
0651 if isfield(newModel,'grRules')
0652 rule=newModel.grRules{nOldRxns+i};
0653 rule=strrep(rule,'(','');
0654 rule=strrep(rule,')','');
0655 rule=strrep(rule,' or ',' ');
0656 rule=strrep(rule,' and ',' ');
0657 genes=regexp(rule,' ','split');
0658 [I, J]=ismember(genes,newModel.genes);
0659 if ~all(I) && any(rule)
0660 EM=['Not all genes for reaction ' rxnsToAdd.rxns{i} ' were found in model.genes. If needed, add genes with addGenesRaven before calling this function, or set the ''allowNewGenes'' flag to true'];
0661 dispEM(EM);
0662 end
0663 end
0664 end
0665
0666
0667
0668 newRxnsModel.genes=newModel.genes;
0669 newRxnsModel.grRules=newModel.grRules(length(model.rxns)+1:end);
0670 newRxnsModel.rxns=newModel.rxns(length(model.rxns)+1:end);
0671
0672
0673 [grRules,rxnGeneMat] = standardizeGrRules(newRxnsModel,true);
0674 newModel.rxnGeneMat = [newModel.rxnGeneMat; rxnGeneMat];
0675 newModel.grRules = [newModel.grRules(1:nOldRxns); grRules];
0676 end