Home > RAVEN > addRxns.m

addRxns

PURPOSE ^

addRxns

SYNOPSIS ^

function newModel=addRxns(model,rxnsToAdd,eqnType,compartment,allowNewMets)

DESCRIPTION ^

 addRxns
   Adds reactions to a model

   model            a model structure
   rxnsToAdd        the reaction structure can have the following fields:
                    rxns       cell array with unique strings that 
                               identifies each reaction
                    equations  cell array with equation strings. Decimal 
                               coefficients are expressed as “1.2”. 
                               Reversibility is indicated by “<=>” or “=>”
                    rxnNames   cell array with the names of each reaction
                               (opt, default '')
                    lb         vector with the lower bounds (opt, default
                               -inf for reversible reactions and 0 for
                               irreversible)
                    ub         vector with the upper bounds (opt, default
                               inf)
                    c          vector with the objective function
                               coefficients (opt, default 0)
                    eccodes    cell array with the EC-numbers for each
                               reactions. Delimit several EC-numbers with
                               ";" (opt, default '')
                    subSystems cell array with the subsystems for each
                               reaction (opt, default '')
                    grRules    cell array with the gene-reaction
                               relationship for each reaction. For example
                               "(A and B) or (C)" means that the reaction 
                               could be catalyzed by a complex between 
                               A & B or by C on its own. All the genes 
                               have to be present in model.genes. Add 
                               genes with addGenes before calling this
                               function if needed (opt, default '')
   eqnType          double describing how the equation string should be
                    interpreted
                    1 - The metabolites are matched to model.mets. New
                        metabolites (if allowed) are added to
                        "compartment"
                    2 - The metabolites are matched to model.metNames and
                        all metabolites are assigned to "compartment". Any
                        new metabolites that are added will be assigned
                        IDs "m1", "m2"... If IDs on the same form are 
                        already used in the model then the numbering will
                        start from the highest used integer+1
                    3 - The metabolites are written as 
                        "metNames[comps]". Only compartments in
                        model.comps are allowed. Any
                        new metabolites that are added will be assigned
                        IDs "m1", "m2"... If IDs on the same form are 
                        already used in the model then the numbering will
                        start from the highest used integer+1  
   compartment      a string with the compartment the metabolites should
                    be placed in when using eqnType=2. Must match 
                    model.comps (opt when eqnType=1 or eqnType=3) 
   allowNewMets     true if the function is allowed to add new
                    metabolites. It is highly recommended to first add
                    any new metabolites with addMets rather than
                    automatically through this function. addMets supports
                    more annotation of metabolites, allows for the use of
                    exchange metabolites, and using it reduces the risk
                    of parsing errors (opt, default false)
                     
   newModel         an updated model structure

   NOTE: This function does not make extensive checks about formatting of
   gene-reaction rules. 

   NOTE: When adding metabolites to a compartment where it previously
   doesn't exist, the function will copy any available information from
   the metabolite in another compartment.

   Usage: newModel=addRxns(model,rxnsToAdd,eqnType,compartment,allowNewMets)

   Rasmus Agren, 2013-02-06

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function newModel=addRxns(model,rxnsToAdd,eqnType,compartment,allowNewMets)
0002 % addRxns
0003 %   Adds reactions to a model
0004 %
0005 %   model            a model structure
0006 %   rxnsToAdd        the reaction structure can have the following fields:
0007 %                    rxns       cell array with unique strings that
0008 %                               identifies each reaction
0009 %                    equations  cell array with equation strings. Decimal
0010 %                               coefficients are expressed as “1.2”.
0011 %                               Reversibility is indicated by “<=>” or “=>”
0012 %                    rxnNames   cell array with the names of each reaction
0013 %                               (opt, default '')
0014 %                    lb         vector with the lower bounds (opt, default
0015 %                               -inf for reversible reactions and 0 for
0016 %                               irreversible)
0017 %                    ub         vector with the upper bounds (opt, default
0018 %                               inf)
0019 %                    c          vector with the objective function
0020 %                               coefficients (opt, default 0)
0021 %                    eccodes    cell array with the EC-numbers for each
0022 %                               reactions. Delimit several EC-numbers with
0023 %                               ";" (opt, default '')
0024 %                    subSystems cell array with the subsystems for each
0025 %                               reaction (opt, default '')
0026 %                    grRules    cell array with the gene-reaction
0027 %                               relationship for each reaction. For example
0028 %                               "(A and B) or (C)" means that the reaction
0029 %                               could be catalyzed by a complex between
0030 %                               A & B or by C on its own. All the genes
0031 %                               have to be present in model.genes. Add
0032 %                               genes with addGenes before calling this
0033 %                               function if needed (opt, default '')
0034 %   eqnType          double describing how the equation string should be
0035 %                    interpreted
0036 %                    1 - The metabolites are matched to model.mets. New
0037 %                        metabolites (if allowed) are added to
0038 %                        "compartment"
0039 %                    2 - The metabolites are matched to model.metNames and
0040 %                        all metabolites are assigned to "compartment". Any
0041 %                        new metabolites that are added will be assigned
0042 %                        IDs "m1", "m2"... If IDs on the same form are
0043 %                        already used in the model then the numbering will
0044 %                        start from the highest used integer+1
0045 %                    3 - The metabolites are written as
0046 %                        "metNames[comps]". Only compartments in
0047 %                        model.comps are allowed. Any
0048 %                        new metabolites that are added will be assigned
0049 %                        IDs "m1", "m2"... If IDs on the same form are
0050 %                        already used in the model then the numbering will
0051 %                        start from the highest used integer+1
0052 %   compartment      a string with the compartment the metabolites should
0053 %                    be placed in when using eqnType=2. Must match
0054 %                    model.comps (opt when eqnType=1 or eqnType=3)
0055 %   allowNewMets     true if the function is allowed to add new
0056 %                    metabolites. It is highly recommended to first add
0057 %                    any new metabolites with addMets rather than
0058 %                    automatically through this function. addMets supports
0059 %                    more annotation of metabolites, allows for the use of
0060 %                    exchange metabolites, and using it reduces the risk
0061 %                    of parsing errors (opt, default false)
0062 %
0063 %   newModel         an updated model structure
0064 %
0065 %   NOTE: This function does not make extensive checks about formatting of
0066 %   gene-reaction rules.
0067 %
0068 %   NOTE: When adding metabolites to a compartment where it previously
0069 %   doesn't exist, the function will copy any available information from
0070 %   the metabolite in another compartment.
0071 %
0072 %   Usage: newModel=addRxns(model,rxnsToAdd,eqnType,compartment,allowNewMets)
0073 %
0074 %   Rasmus Agren, 2013-02-06
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 %If no reactions should be added
0087 if isempty(rxnsToAdd)
0088     return;
0089 end
0090 
0091 %Check the input
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     %To fit with some later printing
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     %It could also be a string, but it's not encouraged
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     %It could also be a string, but it's not encouraged
0131     throw(MException('','rxnsToAdd.equations must be a cell array of strings'));
0132 else
0133     rxnsToAdd.equations=cellstr(rxnsToAdd.equations);
0134 end
0135 
0136 %Check some formatting
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     %If the mets are metNames
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 %***Add everything to the model except for the equations.
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 %Parse the equations. This is done at this early stage since I need the
0174 %reversibility info
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    %Fill with standard if it doesn't exist
0191    if ~isfield(newModel,'rxnNames')
0192        newModel.rxnNames=largeFiller;
0193    end
0194    newModel.rxnNames=[newModel.rxnNames;rxnsToAdd.rxnNames(:)];
0195 else
0196     %Fill with standard if it doesn't exist
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    %Fill with standard if it doesn't exist
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     %Fill with standard if it doesn't exist
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    %Fill with standard if it doesn't exist
0226    if ~isfield(newModel,'ub')
0227        newModel.ub=inf(nOldRxns,1);
0228    end
0229    newModel.ub=[newModel.ub;rxnsToAdd.ub(:)];
0230 else
0231     %Fill with standard if it doesn't exist
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    %Fill with standard if it doesn't exist
0242    if ~isfield(newModel,'c')
0243        newModel.c=zeros(nOldRxns,1);
0244    end
0245    newModel.c=[newModel.c;rxnsToAdd.c(:)];
0246 else
0247     %Fill with standard if it doesn't exist
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    %Fill with standard if it doesn't exist
0258    if ~isfield(newModel,'eccodes')
0259        newModel.eccodes=largeFiller;
0260    end
0261    newModel.eccodes=[newModel.eccodes;rxnsToAdd.eccodes(:)];
0262 else
0263     %Fill with standard if it doesn't exist
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    %Fill with standard if it doesn't exist
0274    if ~isfield(newModel,'subSystems')
0275        newModel.subSystems=largeFiller;
0276    end
0277    newModel.subSystems=[newModel.subSystems;rxnsToAdd.subSystems(:)];
0278 else
0279     %Fill with standard if it doesn't exist
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    %Fill with standard if it doesn't exist
0289    if ~isfield(newModel,'rxnMiriams')
0290        newModel.rxnMiriams=cell(nOldRxns,1);
0291    end
0292    newModel.rxnMiriams=[newModel.rxnMiriams;rxnsToAdd.rxnMiriams(:)];
0293 else
0294     %Fill with standard if it doesn't exist
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    %Fill with standard if it doesn't exist
0304    if ~isfield(newModel,'grRules')
0305        newModel.grRules=largeFiller;
0306    end
0307    newModel.grRules=[newModel.grRules;rxnsToAdd.grRules(:)];
0308 else
0309     %Fill with standard if it doesn't exist
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    %Fill with standard if it doesn't exist
0320    if ~isfield(newModel,'rxnFrom')
0321        newModel.rxnFrom=largeFiller;
0322    end
0323    newModel.rxnFrom=[newModel.rxnFrom;rxnsToAdd.rxnFrom(:)];
0324 else
0325     %Fill with standard if it doesn't exist
0326    if isfield(newModel,'rxnFrom')
0327        newModel.rxnFrom=[newModel.rxnFrom;filler];
0328    end
0329 end
0330 
0331 %Check that ids contain no weird characters. This is only done if the
0332 %equations are with metabolite ids
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     %If the mets are metNames
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 %***Start parsing the equations and adding the info to the S matrix
0359 %The mets are matched to model.mets
0360 if eqnType==1
0361     [I J]=ismember(mets,model.mets);
0362     if ~all(I)
0363         if allowNewMets==true
0364             %Add the new mets
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     %Calculate the indexes of the metabolites and add the info
0374     metIndexes=J;
0375     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0376 end
0377 
0378 %Do some stuff that is the same for eqnType=2 and eqnType=3
0379 if eqnType==2 || eqnType==3
0380     %For later..
0381     t2=strcat(model.metNames,'¤¤¤',model.comps(model.metComps));
0382 end
0383 
0384 %The mets are matched to model.metNames and assigned to "compartment"
0385 if eqnType==2
0386     %%Check that the metabolite names aren't present in the same compartment.
0387     %Not the neatest way maybe..
0388     t1=strcat(mets,'¤¤¤',compartment);
0389     [I J]=ismember(t1,t2);
0390  
0391     if ~all(I)
0392         if allowNewMets==true
0393             %Add the new mets
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     %Calculate the indexes of the metabolites
0403     metIndexes=J;
0404     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0405 end
0406 
0407 %The equations are on the form metNames[compName]
0408 if eqnType==3
0409     %Parse the metabolite names
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         %Check that the formatting is correct
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         %Check that the compartment is correct
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     %Check if the metabolite exists already
0431     t1=strcat(metNames,'¤¤¤',compartments);
0432     [I J]=ismember(t1,t2);
0433  
0434     if ~all(I)
0435         if allowNewMets==true
0436             %Add the new mets
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     %Calculate the indexes of the metabolites
0446     metIndexes=J;
0447     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0448 end
0449 
0450 %Add the info to the stoichiometric matrix and to the rxnGeneMat. I do this
0451 %in a loop, but maybe not necessary
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     %Parse the grRules and add to rxnGeneMat
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             %If there are no genes for the reaction, the rxnGeneMat should
0472             %still be resized
0473             newModel.rxnGeneMat=[newModel.rxnGeneMat;sparse(1,numel(newModel.genes))];
0474        end
0475     end
0476 end
0477 end

Generated on Tue 16-Jul-2013 21:50:02 by m2html © 2005