addReactionGEM

PURPOSE ^

addReactionGEM manually adds reactions to a specified model, may add one or more reactions at a time

SYNOPSIS ^

function [newmodel, HTABLE] = addReactionGEM(model,rxns,rxnNames,rxnFormulas,rev,lb,ub,nRxn,subSystems,grRules,rules,genes, HTABLE)

DESCRIPTION ^

addReactionGEM manually adds reactions to a specified model, may add one or more reactions at a time 

   [newmodel] = addReactionSmiley(model,rxns,rxnNames,rxnFormulas,rev,lb,ub,subSystems,grRules,rules,genes)

 - Manually add reactions to a specified model, can either add one or
   multiple reactions at a time
 - All syntax standards must comply with the specified model
 - For reaction formulas, use: '-->' for irreversible or '<==>' for
   reversible
 e.g. [modelLB_NH3] = addReactionSmiley(modelLB,'NH3r','NH3 protonization',cellstr('1 NH3[c] + 1 H[c] <==> 1 NH4[c]'),1,-1000,1000,'Others');
 Inputs
     model
     rxns
     rxnNames
     rxnFormulas     see note above
     rev             0 = irrev, 1 = rev
     lb
     ub
     subSystems      default = ''
     grRules         default = ''
     rules           default = ''
     genes           default = ''
 Output
     newmodel

 based on AddRxn: Aarash Bordbar 11/2/07
 IT 11-19-08

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [newmodel, HTABLE] = addReactionGEM(model,rxns,rxnNames,rxnFormulas,rev,lb,ub,nRxn,subSystems,grRules,rules,genes, HTABLE)
0002 %addReactionGEM manually adds reactions to a specified model, may add one or more reactions at a time
0003 %
0004 %   [newmodel] = addReactionSmiley(model,rxns,rxnNames,rxnFormulas,rev,lb,ub,subSystems,grRules,rules,genes)
0005 %
0006 % - Manually add reactions to a specified model, can either add one or
0007 %   multiple reactions at a time
0008 % - All syntax standards must comply with the specified model
0009 % - For reaction formulas, use: '-->' for irreversible or '<==>' for
0010 %   reversible
0011 % e.g. [modelLB_NH3] = addReactionSmiley(modelLB,'NH3r','NH3 protonization',cellstr('1 NH3[c] + 1 H[c] <==> 1 NH4[c]'),1,-1000,1000,'Others');
0012 % Inputs
0013 %     model
0014 %     rxns
0015 %     rxnNames
0016 %     rxnFormulas     see note above
0017 %     rev             0 = irrev, 1 = rev
0018 %     lb
0019 %     ub
0020 %     subSystems      default = ''
0021 %     grRules         default = ''
0022 %     rules           default = ''
0023 %     genes           default = ''
0024 % Output
0025 %     newmodel
0026 %
0027 % based on AddRxn: Aarash Bordbar 11/2/07
0028 % IT 11-19-08
0029 
0030 if ~exist('subSystems', 'var') || isempty(subSystems)
0031     for i = 1:length(rev)
0032         subSystems(i,1) = {''};
0033     end
0034 end
0035 
0036 if ~exist('grRules', 'var') || isempty(grRules)
0037     for i = 1:length(rev)
0038         grRules(i,1) = {''};
0039     end
0040 end
0041 
0042 if ~exist('rules', 'var')|| isempty(rules)
0043     for i = 1:length(rev)
0044         rules(i,1) = {''};
0045     end
0046 end
0047 
0048 if ~exist('genes', 'var')|| isempty(genes)
0049     genes = {''};
0050 end
0051 
0052 if ~exist('nRxn', 'var') || isempty(nRxn)
0053     nRxn = length(model.lb)+1;
0054 end
0055 
0056 nMet = length(model.mets)+1;
0057 
0058 if (isfield(model,'genes'))
0059     nGenes = length(model.genes)+1;
0060 end
0061 newmodel = model;
0062 
0063 putNecessary = false;
0064 if exist('HTABLE', 'var') || length(rxnNames) > 2
0065     useHashTable = true;
0066     if exist('HTABLE', 'var') % no need to initialize HTABLE
0067         putNecessary = true; % will need it for output as well.
0068     else %create new hashtable
0069         HTABLE = java.util.Hashtable; % initialize HTABLE
0070         for i = 1:length(newmodel.mets)
0071             HTABLE.put(newmodel.mets{i}, i);
0072         end
0073     end
0074 else
0075     useHashTable = false;
0076 end
0077 
0078 
0079 
0080 
0081 %h = waitbar(0, 'Adding Rxns ...');
0082 for i = 1:length(rev)
0083     %IO indicates whether it is part of the reactants or part of the
0084     %product (ie. IO = -1 is reactant, IO = 1 products)
0085     IO = -1;
0086     %newmodel.rxns{nRxn} = char(rxns(i));
0087     newmodel.lb(nRxn,1) = lb(i,1);
0088     newmodel.ub(nRxn,1) = ub(i,1);
0089     newmodel.rev(nRxn,1) = rev(i,1);
0090     newmodel.subSystems{nRxn,1} = subSystems(i,1);
0091     newmodel.grRules(nRxn,1) = grRules(i,1);
0092     newmodel.rules(nRxn,1) = rules(i,1);
0093     newmodel.c(nRxn,1) = 0;
0094     newmodel.rxnNames{nRxn,1} = char(rxnNames(i,1));
0095     
0096     
0097     %newmodel.S(:, i) = zeros(length(newmodel.mets), 1);
0098     
0099     %parses reaction formula into components of string
0100     [parsing{1,1},parsing{2,1}] = strtok(rxnFormulas{i});
0101     for j = 2:100
0102         [parsing{j,1},parsing{j+1,1}] = strtok(parsing{j,1});
0103         if isempty(parsing{j+1,1})==1
0104             break
0105         end
0106     end
0107 
0108     nComponents = length(parsing)-1;
0109     j = 1;
0110     while j <= nComponents
0111         if strcmp(parsing{j},'+') == 1
0112             j = j+1;
0113         elseif strcmp(parsing{j},' ') == 1
0114             j = j+1;
0115         
0116         %if its '-->' or '<==>' then switches IO to 1 which indicates the
0117         %next compounds are part of the product
0118         elseif strcmp(parsing{j},'<==>') == 1
0119             IO = 1;
0120             j = j+1;
0121         elseif strcmp(parsing{j},'-->') == 1
0122             IO = 1;
0123             j = j+1;
0124         elseif str2double(parsing{j}) > 0
0125 
0126             j = j+1;
0127         else
0128             %if met already exists in newmodel then sets metLoc to
0129             %corresponding met location
0130             
0131             if useHashTable
0132                 c = HTABLE.get(parsing{j});
0133             else
0134                 c = strmatch(parsing{j},newmodel.mets,'exact');
0135                 display('slow method')
0136             end
0137 %             c2 = strmatch(parsing{j},newmodel.mets,'exact');
0138 %
0139 %             %c = HTABLE.get(parsing{j});
0140 %             if any(c2 ~= c)
0141 %                 c
0142 %                 c2
0143 %                 newmodel.mets
0144 %                 pause;
0145 %             end
0146             
0147             %c
0148 %             c = zeros(1,0);
0149 %             ptemp = parsing{j};
0150 %             for k2 = 1:length(newmodel.mets)
0151 %                if strcmp(newmodel.mets{k2}, ptemp)
0152 %                    c = k2;
0153 %                    break;
0154 %                end
0155 %             end
0156             %c
0157             %pause;
0158             %c = strcmp(parsing{j},newmodel.mets);
0159             if ~isempty(c)
0160                 metLoc = c;
0161                 parsing(j);
0162                 
0163             %if met doesn't exist then metLoc is set at nMet (the end of
0164             %list)
0165             else
0166                 metLoc = nMet;
0167                 newmodel.mets{metLoc} = parsing{j};
0168                 newmodel.S(metLoc,:) = 0;
0169                 newmodel.b(metLoc) = 0;
0170                 if putNecessary
0171                     HTABLE.put(parsing{j}, metLoc);
0172                 end
0173                 nMet = nMet+1;
0174             end
0175             %finds reaction coefficient of met
0176             rxnCoeff = str2double(parsing{j-1});
0177             if isnan(rxnCoeff)
0178                 rxnCoeff=1;
0179                 rxnNames(i,1);
0180             end
0181             s = size(model.S);
0182             if(metLoc <= s(1))
0183                 origCoeff = newmodel.S(metLoc,nRxn);
0184                 if(origCoeff < 0)
0185                     newmodel.rxns(nRxn);
0186                 end
0187             else
0188                 origCoeff = 0;
0189             end
0190             newmodel.S(metLoc,nRxn) = rxnCoeff*IO + origCoeff;
0191             
0192             
0193             j=j+1;
0194         end
0195     end
0196 
0197     clear parsing
0198 
0199     
0200     [parsing{1,1},parsing{2,1}] = strtok(grRules{i});
0201     if ~isempty(parsing{2,1}) %length(parsing{2,1}) ~= 0
0202         for j = 2:100
0203             [parsing{j,1},parsing{j+1,1}] = strtok(parsing{j,1});
0204             if isempty(parsing{j+1,1})==1
0205                 break
0206             end
0207         end
0208     end
0209     nRxn = nRxn + 1;
0210     clear parsing
0211     %    waitbar(i/length(rev),h);
0212 end
0213 %close(h);
0214 
0215 for i = 1:length(genes)
0216     if ~isempty(genes{i}) %length(genes{i}) ~= 0
0217         newmodel.genes(nGenes,1) = genes(i,1);
0218     end
0219 end

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