0001 function [model,rxnIDexists] = addReaction(model,rxnName,metaboliteList,stoichCoeffList,revFlag,lowerBound,upperBound,objCoeff,subSystem,grRule,geneNameList,systNameList,checkDuplicate)
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 parseFormulaFlag = false;
0054 rxnIDexists = [];
0055
0056 if iscell(rxnName)&&length(rxnName)>1
0057 rxnNameFull = rxnName{2};
0058 rxnName = rxnName{1};
0059 end
0060
0061
0062 nRxns = length(model.rxns);
0063 if (sum(strcmp(rxnName,model.rxns)) > 0)
0064 warning('Reaction with the same name already exists in the model');
0065 [tmp,rxnID] = ismember(rxnName,model.rxns);
0066 oldRxnFlag = true;
0067 else
0068 rxnID = nRxns+1;
0069 oldRxnFlag = false;
0070 end
0071
0072
0073 if (nargin < 4)
0074 if (~iscell(metaboliteList))
0075 parseFormulaFlag = true;
0076 else
0077 error('Missing stoichiometry information');
0078 end
0079 else
0080 if isempty(stoichCoeffList)
0081 parseFormulaFlag = true;
0082 else
0083 if (length(metaboliteList) ~= length(stoichCoeffList))
0084 error('Incorrect number of stoichiometric coefficients provided');
0085 end
0086 end
0087 end
0088
0089
0090 if (nargin < 5 | isempty(revFlag))
0091 if (oldRxnFlag)
0092 revFlag = model.rev(rxnID);
0093 else
0094 revFlag = true;
0095 end
0096 end
0097
0098
0099 if (parseFormulaFlag)
0100 rxnFormula = metaboliteList;
0101 [metaboliteList,stoichCoeffList,revFlag] = parseRxnFormula(rxnFormula);
0102 end
0103
0104
0105 if (nargin < 6 | isempty(lowerBound))
0106 if (oldRxnFlag)
0107 lowerBound = model.lb(rxnID);
0108 else
0109 if (revFlag)
0110 lowerBound = min(model.lb);
0111 if isempty(lowerBound)
0112 lowerBound=-1000;
0113 end
0114 else
0115 lowerBound = 0;
0116 end
0117 end
0118 end
0119 if (nargin < 7 | isempty(upperBound))
0120 if (oldRxnFlag)
0121 upperBound = model.ub(rxnID);
0122 else
0123 upperBound = max(model.ub);
0124 if isempty(upperBound)
0125 upperBound=1000;
0126 end
0127 end
0128 end
0129 if (nargin < 8 | isempty(objCoeff))
0130 if (oldRxnFlag)
0131 objCoeff = model.c(rxnID);
0132 else
0133 objCoeff = 0;
0134 end
0135 end
0136 if (nargin < 9)
0137 if (oldRxnFlag) && (isfield(model,'subSystems'))
0138 subSystem = model.subSystems{rxnID};
0139 else
0140 subSystem = '';
0141 end
0142 end
0143 if isempty(subSystem)
0144 if (oldRxnFlag) && (isfield(model,'subSystems'))
0145 subSystem = model.subSystems{rxnID};
0146 else
0147 subSystem = '';
0148 end
0149 end
0150 if (nargin < 10) && (isfield(model,'grRules'))
0151 if (oldRxnFlag)
0152 grRule = model.grRules{rxnID};
0153 else
0154 grRule = '';
0155 end
0156 end
0157
0158 if (~exist('checkDuplicate','var'))
0159 checkDuplicate=true;
0160 end
0161
0162 nMets = length(model.mets);
0163 Scolumn = sparse(nMets,1);
0164
0165 modelOrig = model;
0166
0167
0168 model.rxns{rxnID,1} = rxnName;
0169 if (revFlag)
0170 model.rev(rxnID,1) = 1;
0171 else
0172 model.rev(rxnID,1) = 0;
0173 end
0174 model.lb(rxnID,1) = lowerBound;
0175 model.ub(rxnID,1) = upperBound;
0176 model.c(rxnID,1) = objCoeff;
0177
0178 if (isfield(model,'rxnNames'))
0179 if exist('rxnNameFull','var')
0180 model.rxnNames{rxnID,1} = rxnNameFull;
0181 else
0182 model.rxnNames{rxnID,1} = model.rxns{rxnID};
0183 end
0184 end
0185 if (isfield(model,'subSystems'))
0186 model.subSystems{rxnID,1} = subSystem;
0187 end
0188 if isfield(model,'rxnNotes')
0189 model.rxnNotes{rxnID,1} = '';
0190 end
0191 if isfield(model,'confidenceScores')
0192 model.confidenceScores{rxnID,1} = '';
0193 end
0194 if isfield(model,'rxnReferences')
0195 model.rxnReferences{rxnID,1} = '';
0196 end
0197 if isfield(model,'rxnECNumbers')
0198 model.rxnECNumbers{rxnID,1} = '';
0199 end
0200
0201
0202
0203 [isInModel,metID] = ismember(metaboliteList,model.mets);
0204
0205 nNewMets = sum(~isInModel);
0206
0207
0208 newMetsCoefs=zeros(0);
0209 for i = 1:length(metaboliteList)
0210 if (isInModel(i))
0211 Scolumn(metID(i),1) = stoichCoeffList(i);
0212 else
0213 warning(['Metabolite ' metaboliteList{i} ' not in model - added to the model']);
0214 Scolumn(end+1,1) = stoichCoeffList(i);
0215 model.mets{end+1,1} = metaboliteList{i};
0216 newMetsCoefs(end+1) = stoichCoeffList(i);
0217 if (isfield(model,'metNames'))
0218 model.metNames{end+1,1} = regexprep(metaboliteList{i},'(\[.+\]) | (\(.+\))','') ;
0219 warning(['Metabolite name for ' metaboliteList{i} ' set to ' model.metNames{end}]);
0220
0221 end
0222 if (isfield(model,'metFormulas'))
0223 model.metFormulas{end+1,1} = '';
0224 warning(['Metabolite formula for ' metaboliteList{i} ' set to ''''']);
0225
0226 end
0227 if isfield(model,'metChEBIID')
0228 model.metChEBIID{end+1,1} = '';
0229 end
0230 if isfield(model,'metKEGGID')
0231 model.metKEGGID{end+1,1} = '';
0232 end
0233 if isfield(model,'metPubChemID')
0234 model.metPubChemID{end+1,1} = '';
0235 end
0236 if isfield(model,'metInChIString')
0237 model.metInChIString{end+1,1} = '';
0238 end
0239 if isfield(model,'metCharge')
0240 model.metCharge(end+1,1) = 0;
0241 end
0242 end
0243 end
0244
0245
0246
0247 if isfield(model,'b')
0248 model.b = [model.b;zeros(length(model.mets)-length(model.b),1)];
0249 end
0250
0251
0252
0253 if (isfield(model,'genes'))
0254 if (nargin < 11)
0255 model = changeGeneAssociation(model,rxnName,grRule);
0256 else
0257 model = changeGeneAssociation(model,rxnName,grRule,geneNameList,systNameList);
0258 end
0259 end
0260
0261
0262 rxnInModel=false;
0263 if (nNewMets > 0) && isempty(find(newMetsCoefs == 0, 1))
0264 Stmp = [model.S;sparse(nNewMets,nRxns)];
0265 else
0266 Stmp = model.S;
0267 if (checkDuplicate)
0268 if size(Stmp,2)<6000
0269 tmpSel = all(repmat((Scolumn),1,size(Stmp,2)) == (Stmp));
0270 rxnIDexists = full(find(tmpSel));
0271 if (~isempty(rxnIDexists))
0272 rxnIDexists=rxnIDexists(1);
0273 rxnInModel = true;
0274 end
0275 else
0276 for i=1:size(Stmp,2)
0277 if(Scolumn==Stmp(:,i))
0278 rxnInModel=true;
0279 rxnIDexists=i;
0280 break
0281 end
0282 end
0283 end
0284 end
0285 end
0286
0287 if (rxnInModel)
0288 warning(['Model already has the same reaction you tried to add: ' modelOrig.rxns{rxnIDexists}]);
0289 model = modelOrig;
0290 else
0291 if (oldRxnFlag)
0292 model.S = Stmp;
0293 model.S(:,rxnID) = Scolumn;
0294 else
0295 model.S = [Stmp Scolumn];
0296 end
0297 printRxnFormula(model,rxnName);
0298 end