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 '')
                    rxnMiriams cell array with Miriam structures (opt,
                               default [])
                    rxnComps   cell array with compartments (as in
                               model.comps) (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, 2014-01-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 %                    rxnMiriams cell array with Miriam structures (opt,
0035 %                               default [])
0036 %                    rxnComps   cell array with compartments (as in
0037 %                               model.comps) (opt, default {})
0038 %   eqnType          double describing how the equation string should be
0039 %                    interpreted
0040 %                    1 - The metabolites are matched to model.mets. New
0041 %                        metabolites (if allowed) are added to
0042 %                        "compartment"
0043 %                    2 - The metabolites are matched to model.metNames and
0044 %                        all metabolites are assigned to "compartment". Any
0045 %                        new metabolites that are added will be assigned
0046 %                        IDs "m1", "m2"... If IDs on the same form are
0047 %                        already used in the model then the numbering will
0048 %                        start from the highest used integer+1
0049 %                    3 - The metabolites are written as
0050 %                        "metNames[comps]". Only compartments in
0051 %                        model.comps are allowed. Any
0052 %                        new metabolites that are added will be assigned
0053 %                        IDs "m1", "m2"... If IDs on the same form are
0054 %                        already used in the model then the numbering will
0055 %                        start from the highest used integer+1
0056 %   compartment      a string with the compartment the metabolites should
0057 %                    be placed in when using eqnType=2. Must match
0058 %                    model.comps (opt when eqnType=1 or eqnType=3)
0059 %   allowNewMets     true if the function is allowed to add new
0060 %                    metabolites. It is highly recommended to first add
0061 %                    any new metabolites with addMets rather than
0062 %                    automatically through this function. addMets supports
0063 %                    more annotation of metabolites, allows for the use of
0064 %                    exchange metabolites, and using it reduces the risk
0065 %                    of parsing errors (opt, default false)
0066 %
0067 %   newModel         an updated model structure
0068 %
0069 %   NOTE: This function does not make extensive checks about formatting of
0070 %   gene-reaction rules.
0071 %
0072 %   NOTE: When adding metabolites to a compartment where it previously
0073 %   doesn't exist, the function will copy any available information from
0074 %   the metabolite in another compartment.
0075 %
0076 %   Usage: newModel=addRxns(model,rxnsToAdd,eqnType,compartment,allowNewMets)
0077 %
0078 %   Rasmus Agren, 2014-01-06
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 %If no reactions should be added
0091 if isempty(rxnsToAdd)
0092     return;
0093 end
0094 
0095 %Check the input
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     %To fit with some later printing
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     %It could also be a string, but it's not encouraged
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     %It could also be a string, but it's not encouraged
0135     dispEM('rxnsToAdd.equations must be a cell array of strings');
0136 else
0137     rxnsToAdd.equations=cellstr(rxnsToAdd.equations);
0138 end
0139 
0140 %Check some formatting
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 %***Add everything to the model except for the equations.
0157 if numel(rxnsToAdd.equations)~=nRxns
0158    dispEM('rxnsToAdd.equations must have the same number of elements as rxnsToAdd.rxns');
0159 end
0160 
0161 %Parse the equations. This is done at this early stage since I need the
0162 %reversibility info
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    %Fill with standard if it doesn't exist
0174    if ~isfield(newModel,'rxnNames')
0175        newModel.rxnNames=largeFiller;
0176    end
0177    newModel.rxnNames=[newModel.rxnNames;rxnsToAdd.rxnNames(:)];
0178 else
0179     %Fill with standard if it doesn't exist
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    %Fill with standard if it doesn't exist
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     %Fill with standard if it doesn't exist
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    %Fill with standard if it doesn't exist
0209    if ~isfield(newModel,'ub')
0210        newModel.ub=inf(nOldRxns,1);
0211    end
0212    newModel.ub=[newModel.ub;rxnsToAdd.ub(:)];
0213 else
0214     %Fill with standard if it doesn't exist
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    %Fill with standard if it doesn't exist
0225    if ~isfield(newModel,'c')
0226        newModel.c=zeros(nOldRxns,1);
0227    end
0228    newModel.c=[newModel.c;rxnsToAdd.c(:)];
0229 else
0230     %Fill with standard if it doesn't exist
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    %Fill with standard if it doesn't exist
0241    if ~isfield(newModel,'eccodes')
0242        newModel.eccodes=largeFiller;
0243    end
0244    newModel.eccodes=[newModel.eccodes;rxnsToAdd.eccodes(:)];
0245 else
0246     %Fill with standard if it doesn't exist
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    %Fill with standard if it doesn't exist
0257    if ~isfield(newModel,'subSystems')
0258        newModel.subSystems=largeFiller;
0259    end
0260    newModel.subSystems=[newModel.subSystems;rxnsToAdd.subSystems(:)];
0261 else
0262     %Fill with standard if it doesn't exist
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    %Fill with standard if it doesn't exist
0272    if ~isfield(newModel,'rxnMiriams')
0273        newModel.rxnMiriams=cell(nOldRxns,1);
0274    end
0275    newModel.rxnMiriams=[newModel.rxnMiriams;rxnsToAdd.rxnMiriams(:)];
0276 else
0277     %Fill with standard if it doesn't exist
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    %Fill with standard if it doesn't exist
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     %Fill with standard if it doesn't exist
0294    if isfield(newModel,'rxnComps')
0295        newModel.rxnComps=[newModel.rxnComps;ones(nRxns,1)];
0296        %fprintf('NOTE: The added reactions will be assigned to the first compartment\n');
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    %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        dispEM('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     dispEM('Illegal character(s) in metabolite IDs:',true,mets(~cellfun(@isempty,illegalCells)));
0336 else
0337     %If the mets are metNames
0338     illegalCells=regexp(mets,'["%<>\\]', 'once');
0339     dispEM('Illegal character(s) in metabolite names:',true,mets(~cellfun(@isempty,illegalCells)));
0340 end
0341 
0342 %***Start parsing the equations and adding the info to the S matrix
0343 %The mets are matched to model.mets
0344 if eqnType==1
0345     [I J]=ismember(mets,model.mets);
0346     if ~all(I)
0347         if allowNewMets==true
0348             %Add the new mets
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     %Calculate the indexes of the metabolites and add the info
0358     metIndexes=J;
0359     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0360 end
0361 
0362 %Do some stuff that is the same for eqnType=2 and eqnType=3
0363 if eqnType==2 || eqnType==3
0364     %For later..
0365     t2=strcat(model.metNames,'¤¤¤',model.comps(model.metComps));
0366 end
0367 
0368 %The mets are matched to model.metNames and assigned to "compartment"
0369 if eqnType==2
0370     %%Check that the metabolite names aren't present in the same compartment.
0371     %Not the neatest way maybe..
0372     t1=strcat(mets,'¤¤¤',compartment);
0373     [I J]=ismember(t1,t2);
0374  
0375     if ~all(I)
0376         if allowNewMets==true
0377             %Add the new mets
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     %Calculate the indexes of the metabolites
0387     metIndexes=J;
0388     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0389 end
0390 
0391 %The equations are on the form metNames[compName]
0392 if eqnType==3
0393     %Parse the metabolite names
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         %Check that the formatting is correct
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         %Check that the compartment is correct
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     %Check if the metabolite exists already
0415     t1=strcat(metNames,'¤¤¤',compartments);
0416     [I J]=ismember(t1,t2);
0417  
0418     if ~all(I)
0419         if allowNewMets==true
0420             %Add the new mets
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     %Calculate the indexes of the metabolites
0430     metIndexes=J;
0431     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0432 end
0433 
0434 %Add the info to the stoichiometric matrix and to the rxnGeneMat. I do this
0435 %in a loop, but maybe not necessary
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     %Parse the grRules and add to rxnGeneMat
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             %If there are no genes for the reaction, the rxnGeneMat should
0456             %still be resized
0457             newModel.rxnGeneMat=[newModel.rxnGeneMat;sparse(1,numel(newModel.genes))];
0458        end
0459     end
0460 end
0461 end

Generated on Mon 06-Jan-2014 14:58:12 by m2html © 2005