addReaction

PURPOSE ^

addReaction Add a reaction to the model or modify an existing reaction

SYNOPSIS ^

function [model,rxnIDexists] = addReaction(model,rxnName,metaboliteList,stoichCoeffList,revFlag,lowerBound,upperBound,objCoeff,subSystem,grRule,geneNameList,systNameList,checkDuplicate)

DESCRIPTION ^

addReaction Add a reaction to the model or modify an existing reaction

 model = addReaction(model,rxnName,metaboliteList,stoichCoeffList,revFlag,lowerBound,upperBound,objCoeff,subSystem,grRule,geneNameList,systNameList,checkDuplicate)
 model = addReaction(model,rxnName,rxnFormula)

INPUTS
 model             COBRA model structure
 rxnName           Reaction name abbreviation (i.e. 'ACALD')
                   (Note: can also be a cell array {'abbr','name'}
 metaboliteList    Cell array of metabolite names or alternatively the
                   reaction formula for the reaction
 stoichCoeffList   List of stoichiometric coefficients (reactants -ve,
                   products +ve), empty if reaction formula is provided

OPTIONAL INPUTS
 revFlag           Reversibility flag (Default = true)
 lowerBound        Lower bound (Default = 0 or -vMax)
 upperBound        Upper bound (Default = vMax)
 objCoeff          Objective coefficient (Default = 0)
 subSystem         Subsystem (Default = '')
 grRule            Gene-reaction rule in boolean format (and/or allowed)
                   (Default = '');
 geneNameList      List of gene names (used only for translation from 
                   common gene names to systematic gene names)
 systNameList      List of systematic names
 checkDuplicate    Check S matrix too see if a duplicate reaction is
                   already in the model (Deafult true)

OUTPUTS
 model             COBRA model structure with new reaction
 rxnIDexists       Empty if the reaction did not exist previously, or if
                   checkDuplicate is false. Otherwise it contains the ID
                   of an identical reaction already present in the model.

 Examples:

 1) Add a new irreversible reaction using the formula approach

    model = addReaction(model,'newRxn1','A -> B + 2 C')

 2) Add a the same reaction using the list approach

    model = addReaction(model,'newRxn1',{'A','B','C'},[-1 1 2],false);

 Markus Herrgard 1/12/07

 Modified the check to see if duplicate reaction already is in model by
 using S matrix coefficients to be able to handle larger matricies
 Richard Que 11/13/2008

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [model,rxnIDexists] = addReaction(model,rxnName,metaboliteList,stoichCoeffList,revFlag,lowerBound,upperBound,objCoeff,subSystem,grRule,geneNameList,systNameList,checkDuplicate)
0002 %addReaction Add a reaction to the model or modify an existing reaction
0003 %
0004 % model = addReaction(model,rxnName,metaboliteList,stoichCoeffList,revFlag,lowerBound,upperBound,objCoeff,subSystem,grRule,geneNameList,systNameList,checkDuplicate)
0005 % model = addReaction(model,rxnName,rxnFormula)
0006 %
0007 %INPUTS
0008 % model             COBRA model structure
0009 % rxnName           Reaction name abbreviation (i.e. 'ACALD')
0010 %                   (Note: can also be a cell array {'abbr','name'}
0011 % metaboliteList    Cell array of metabolite names or alternatively the
0012 %                   reaction formula for the reaction
0013 % stoichCoeffList   List of stoichiometric coefficients (reactants -ve,
0014 %                   products +ve), empty if reaction formula is provided
0015 %
0016 %OPTIONAL INPUTS
0017 % revFlag           Reversibility flag (Default = true)
0018 % lowerBound        Lower bound (Default = 0 or -vMax)
0019 % upperBound        Upper bound (Default = vMax)
0020 % objCoeff          Objective coefficient (Default = 0)
0021 % subSystem         Subsystem (Default = '')
0022 % grRule            Gene-reaction rule in boolean format (and/or allowed)
0023 %                   (Default = '');
0024 % geneNameList      List of gene names (used only for translation from
0025 %                   common gene names to systematic gene names)
0026 % systNameList      List of systematic names
0027 % checkDuplicate    Check S matrix too see if a duplicate reaction is
0028 %                   already in the model (Deafult true)
0029 %
0030 %OUTPUTS
0031 % model             COBRA model structure with new reaction
0032 % rxnIDexists       Empty if the reaction did not exist previously, or if
0033 %                   checkDuplicate is false. Otherwise it contains the ID
0034 %                   of an identical reaction already present in the model.
0035 %
0036 % Examples:
0037 %
0038 % 1) Add a new irreversible reaction using the formula approach
0039 %
0040 %    model = addReaction(model,'newRxn1','A -> B + 2 C')
0041 %
0042 % 2) Add a the same reaction using the list approach
0043 %
0044 %    model = addReaction(model,'newRxn1',{'A','B','C'},[-1 1 2],false);
0045 %
0046 % Markus Herrgard 1/12/07
0047 %
0048 % Modified the check to see if duplicate reaction already is in model by
0049 % using S matrix coefficients to be able to handle larger matricies
0050 % Richard Que 11/13/2008
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 % Figure out if reaction already exists
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 % Figure out what input format is used
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 % Reversibility
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 % Parse formula
0099 if (parseFormulaFlag)
0100     rxnFormula = metaboliteList;
0101     [metaboliteList,stoichCoeffList,revFlag] = parseRxnFormula(rxnFormula);
0102 end
0103 
0104 % Missing arguments
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 % Update model fields
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 % Figure out which metabolites are already in the model
0203 [isInModel,metID] = ismember(metaboliteList,model.mets);
0204 
0205 nNewMets = sum(~isInModel);
0206 
0207 % Construct S-matrix column
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'))      %Prompts to add missing info if desired
0218             model.metNames{end+1,1} = regexprep(metaboliteList{i},'(\[.+\]) | (\(.+\))','') ;
0219             warning(['Metabolite name for ' metaboliteList{i} ' set to ' model.metNames{end}]);
0220 %             model.metNames(end) = cellstr(input('Enter complete metabolite name, if available:', 's'));
0221         end
0222         if (isfield(model,'metFormulas'))
0223             model.metFormulas{end+1,1} = '';
0224             warning(['Metabolite formula for ' metaboliteList{i} ' set to ''''']);
0225 %             model.metFormulas(end) = cellstr(input('Enter metabolite chemical formula, if available:', 's'));
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 %printLabeledData(model.mets,Scolumn,1);
0246 
0247 if isfield(model,'b')
0248     model.b = [model.b;zeros(length(model.mets)-length(model.b),1)];
0249 end
0250 
0251 % if ~oldRxnFlag, model.rxnGeneMat(rxnID,:)=0; end
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 % Figure out if the new reaction already exists
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

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003