getRxnsFromMetaCyc

 getRxnsFromMetaCyc
   Retrieves reactions information from MetaCyc database

   Input:
   metacycPath         if metaCycRxns.mat is not in the RAVEN\external\metacyc
                       directory, this function will attempt to build it by
                       reading info from a local dump of MetaCyc database,
                       metacycPath is the path to the MetaCyc data files
   keepTransportRxns   include transportation reactions, which often have identical
                       reactants and products that turn to be all-zero columns in
                       the S matrix (opt, default false)
   keepUnbalanced      include reactions cannot be balanced, usually
                       because they are polymeric reactions or with
                       specific difficulty in balancing class structures
                       (opt, default false)
   keepUndetermined    include reactions that have substrates lack chemical
                       structures or with non-numerical coefficients (e.g. n+1)
                       (opt, default false)

   Output:
   model     a model structure generated from the database. The following
             fields are filled
             id:             'MetaCyc'
             name:    'Automatically generated from MetaCyc database'
             rxns:           MetaCyc reaction ids
             rxnNames:       Name for each reaction entry
             mets:           MetaCyc compound ids. If the equations use
                             stoichiometry such as ID(n+1) then the whole
                             expression is saved as the id
             pwys:           MetaCyc pathway id
             eccodes:        Corresponding ec-number if available
             rxnMiriams:     Contains reaction specific information such as
                             associated RHEA and KEGG reaction ids
             S:              Stoichiometric matrix
             lb:             -1000 for all reactions
             ub:             1000 for all reactions
             rev:            1 for reversible and 0 for irreversible. For
                             reactions present in pathway maps the reversibility
                             is taken from there
             b:              0 for all metabolites
             version:        MetaCyc database version

   If the file metaCycRxns.mat is in the RAVEN/external/metacyc directory, it
   will be directly loaded instead of parsing the MetaCyc data files and
   pre-prepared lists of MetaCyc transport and undetermined reactions.

   Usage: model=getRxnsFromMetaCyc(metacycPath,keepTransportRxns,keepUnbalanced,keepUndetermined)

0001 function model=getRxnsFromMetaCyc(metacycPath,keepTransportRxns,keepUnbalanced,keepUndetermined)
0002 % getRxnsFromMetaCyc
0003 %   Retrieves reactions information from MetaCyc database
0004 %
0005 %   Input:
0006 %   metacycPath         if metaCycRxns.mat is not in the RAVEN\external\metacyc
0007 %                       directory, this function will attempt to build it by
0008 %                       reading info from a local dump of MetaCyc database,
0009 %                       metacycPath is the path to the MetaCyc data files
0010 %   keepTransportRxns   include transportation reactions, which often have identical
0011 %                       reactants and products that turn to be all-zero columns in
0012 %                       the S matrix (opt, default false)
0013 %   keepUnbalanced      include reactions cannot be balanced, usually
0014 %                       because they are polymeric reactions or with
0015 %                       specific difficulty in balancing class structures
0016 %                       (opt, default false)
0017 %   keepUndetermined    include reactions that have substrates lack chemical
0018 %                       structures or with non-numerical coefficients (e.g. n+1)
0019 %                       (opt, default false)
0020 %
0021 %   Output:
0022 %   model     a model structure generated from the database. The following
0023 %             fields are filled
0024 %             id:             'MetaCyc'
0025 %             name:    'Automatically generated from MetaCyc database'
0026 %             rxns:           MetaCyc reaction ids
0027 %             rxnNames:       Name for each reaction entry
0028 %             mets:           MetaCyc compound ids. If the equations use
0029 %                             stoichiometry such as ID(n+1) then the whole
0030 %                             expression is saved as the id
0031 %             pwys:           MetaCyc pathway id
0032 %             eccodes:        Corresponding ec-number if available
0033 %             rxnMiriams:     Contains reaction specific information such as
0034 %                             associated RHEA and KEGG reaction ids
0035 %             S:              Stoichiometric matrix
0036 %             lb:             -1000 for all reactions
0037 %             ub:             1000 for all reactions
0038 %             rev:            1 for reversible and 0 for irreversible. For
0039 %                             reactions present in pathway maps the reversibility
0040 %                             is taken from there
0041 %             b:              0 for all metabolites
0042 %             version:        MetaCyc database version
0043 %
0044 %   If the file metaCycRxns.mat is in the RAVEN/external/metacyc directory, it
0045 %   will be directly loaded instead of parsing the MetaCyc data files and
0046 %   pre-prepared lists of MetaCyc transport and undetermined reactions.
0047 %
0048 %   Usage: model=getRxnsFromMetaCyc(metacycPath,keepTransportRxns,keepUnbalanced,keepUndetermined)
0049 
0050 %NOTE: This is how one entry looks in the file
0051 
0052 % UNIQUE-ID - 4-HYDROXY-2-KETOPIMELATE-LYSIS-RXN
0053 % TYPES - Small-Molecule-Reactions
0054 % TYPES - Chemical-Reactions
0055 % ATOM-MAPPINGS - (:NO-HYDROGEN-ENCODING (6 7 0 8 1 9 2 10 3 12 11 5 4) (((PYRUVATE 0 5) (SUCC-S-ALD 6 12)) ((CPD-804 0 12))))
0056 % DBLINKS - (LIGAND-RXN "R01647" NIL |kr| 3650824949 NIL NIL)
0057 % DBLINKS - (LIGAND-RXN "R01645" NIL |taltman| 3459474590 NIL NIL)
0058 % EC-NUMBER - EC-4.1.2.52
0059 % ^OFFICIAL? - T
0060 % ENZYMATIC-REACTION - ENZRXN-3009
0061 % IN-PATHWAY - 3-HYDROXYPHENYLACETATE-DEGRADATION-PWY
0062 % LEFT - CPD-804
0063 % ORPHAN? - :NO
0064 % PHYSIOLOGICALLY-RELEVANT? - T
0065 % REACTION-DIRECTION - LEFT-TO-RIGHT
0066 % RIGHT - SUCC-S-ALD
0067 % RIGHT - PYRUVATE
0068 % //
0069 
0070 % A line that contains only '//' separates each object.
0071 if nargin<1
0072     ravenPath=findRAVENroot();
0073     metacycPath=fullfile(ravenPath,'external','metacyc');
0074 else
0075     metacycPath=char(metacycPath);
0076 end
0077 if nargin<2
0078     keepTransportRxns=false;
0079 end
0080 if nargin<3
0081     keepUnbalanced=false;
0082 end
0083 if nargin<4
0084     keepUndetermined=false;
0085 end
0086 
0087 %Check if the reactions have been parsed before and saved. Directly load
0088 %the model if so.
0089 rxnsFile=fullfile(metacycPath,'metaCycRxns.mat');
0090 metaCycRxnFile='reactions.dat';
0091 metaCycPwyFile='pathway-links.dat';
0092 
0093 if exist(rxnsFile, 'file')
0094     fprintf(['Importing MetaCyc reactions from ' strrep(rxnsFile,'\','/') '... ']);
0095     load(rxnsFile);
0096     fprintf('done\n');
0097 else
0098     fprintf(['Cannot locate ' strrep(rxnsFile,'\','/') '\nNow try to generate it from local MetaCyc data files...\n']);
0099     if ~isfile(fullfile(metacycPath,metaCycRxnFile)) || ~isfile(fullfile(metacycPath,metaCycPwyFile))
0100         EM=fprintf(['The files of reactions or pathways cannot be located, and should be downloaded from MetaCyc.\n']);
0101         dispEM(EM);
0102     else
0103         metaCycRxns.id='MetaCyc';
0104         metaCycRxns.name='Automatically generated from MetaCyc database';
0105         
0106         %Get pathway names and add them to the field of subSystems
0107         fid = fopen(fullfile(metacycPath,metaCycPwyFile), 'r');
0108         
0109         %Keeps track of how many pathways that have been added
0110         pwyCounter=0;
0111         pwys=cell(10000,1);
0112         pwyNames=cell(10000,1);
0113         
0114         %Loop through the file
0115         while 1
0116             %Get the next line
0117             tline = fgetl(fid);
0118             
0119             %Abort at end of file
0120             if ~ischar(tline)
0121                 break;
0122             end
0123             
0124             %Read in pathway ids and names
0125             if ~strcmp(tline(1),'#')
0126                 pwyCounter=pwyCounter+1;
0127                 vars = regexp(tline, '\t', 'split');
0128                 pwys{pwyCounter}=vars{1};
0129                 pwyNames{pwyCounter}=vars{2};
0130                 
0131                 %Romve HTML symbols
0132                 pwyNames{pwyCounter}=removeHTMLcodes(pwyNames{pwyCounter});
0133                 
0134             end
0135         end
0136         fclose(fid);
0137         pwys=pwys(1:pwyCounter);
0138         pwyNames=pwyNames(1:pwyCounter);
0139         
0140         %Preallocate memory for 50000 reactions
0141         metaCycRxns.rxns=cell(50000,1);
0142         metaCycRxns.rxnNames=cell(50000,1);
0143         metaCycRxns.eccodes=cell(50000,1);
0144         metaCycRxns.subSystems=cell(50000,1);
0145         metaCycRxns.pwys=cell(50000,1);
0146         metaCycRxns.rxnMiriams=cell(50000,1);
0147         metaCycRxns.rxnReferences=cell(50000,1);
0148         metaCycRxns.rev=ones(50000,1); %reversibility;
0149         
0150         rxnLinks.metacyc=cell(10000,1);
0151         rxnLinks.kegg=cell(10000,1);
0152         rxnLinks.check=cell(10000,1);
0153         isSpontaneous=false(10000,1); %spontaneous;
0154         UNBALANCED=false(10000,1);
0155         UNDETERMINED=false(10000,1);
0156         TRANSPORT={};                 %transport reactions;
0157         
0158         metaCycRxns.equations=cell(50000,1); %reaction equations
0159         left=cell(50000,1); %Temporarily stores the equations
0160         right=cell(50000,1); %Temporarily stores the equations
0161         
0162         %First load information on reaction ID, reaction name,, pathway,
0163         %and ec-number
0164         fid = fopen(fullfile(metacycPath,metaCycRxnFile), 'r');
0165         
0166         %Keeps track of how many reactions that have been added
0167         rxnCounter=0;
0168         addSpont=false;
0169         dbLinkCounter=0;
0170         
0171         %Loop through the file
0172         while 1
0173             %Get the next line
0174             tline = fgetl(fid);
0175             
0176             %Abort at end of file
0177             if ~ischar(tline)
0178                 break;
0179             end
0180 
0181             % Get the version of MetaCyc database
0182             if numel(tline)>11 && strcmp(tline(1:11),'# Version: ')
0183                 version=tline(12:end);
0184             end
0185 
0186             %Check if it is a new reaction
0187             if numel(tline)>12 && strcmp(tline(1:12),'UNIQUE-ID - ')
0188                 rxnCounter=rxnCounter+1;
0189                 nPwys=0;
0190                 
0191                 %Add empty strings where there should be such
0192                 metaCycRxns.rxnNames{rxnCounter}='';
0193                 metaCycRxns.eccodes{rxnCounter}='';
0194                 metaCycRxns.subSystems{rxnCounter}='';
0195                 metaCycRxns.pwys{rxnCounter}='';
0196                 metaCycRxns.equations{rxnCounter}='';
0197                 metaCycRxns.rxnReferences{rxnCounter}='';
0198                 reverse=0;   %Value for reversing the equation when necessary
0199                 
0200                 % A complex evaluation system for generating the equations
0201                 left{rxnCounter}='';
0202                 right{rxnCounter}='';
0203                 coefficient='';
0204                 templeft='';
0205                 tempright='';
0206                 % might be simplified with a better algorithem
0207                 
0208                 %Add reaction ID
0209                 metaCycRxns.rxns{rxnCounter}=tline(13:end);
0210                 
0211             end
0212             
0213             %Add name
0214             if numel(tline)>14 && strcmp(tline(1:14),'COMMON-NAME - ')
0215                 metaCycRxns.rxnNames{rxnCounter}=tline(15:end);
0216                 
0217                 %Romve HTML symbols
0218                 metaCycRxns.rxnNames{rxnCounter}=removeHTMLcodes(metaCycRxns.rxnNames{rxnCounter});
0219             end
0220             
0221             %Add eccodes
0222             if numel(tline)>15 && strcmp(tline(1:15),'EC-NUMBER - EC-')
0223                 if isempty(metaCycRxns.eccodes{rxnCounter})
0224                     metaCycRxns.eccodes{rxnCounter}=strcat('ec-code/',tline(16:end));
0225                 else
0226                     metaCycRxns.eccodes{rxnCounter}=strcat(metaCycRxns.eccodes{rxnCounter},';ec-code/',tline(16:end));
0227                 end
0228             end
0229             if numel(tline)>16 && strcmp(tline(1:16),'EC-NUMBER - |EC-')
0230                 if isempty(metaCycRxns.eccodes{rxnCounter})
0231                     metaCycRxns.eccodes{rxnCounter}=strcat('ec-code/',tline(17:end-1));
0232                 else
0233                     metaCycRxns.eccodes{rxnCounter}=strcat(metaCycRxns.eccodes{rxnCounter},';ec-code/',tline(17:end-1));
0234                 end
0235             end
0236             
0237             %Add pathway id and assign pathway name to subSystem field
0238             if numel(tline)>13 && strcmp(tline(1:13),'IN-PATHWAY - ')
0239                 if isempty(metaCycRxns.pwys{rxnCounter})
0240                     metaCycRxns.pwys{rxnCounter}=tline(14:end);
0241                 else
0242                     metaCycRxns.pwys{rxnCounter}=strcat(metaCycRxns.pwys{rxnCounter},';',tline(14:end));
0243                 end
0244                 
0245                 [x, y]=ismember(tline(14:end),pwys);
0246                 if x
0247                     metaCycRxns.subSystems{rxnCounter,1}{1,numel(metaCycRxns.subSystems{rxnCounter,1})+1}=pwyNames{y};
0248                 end
0249             end
0250             
0251             %Add references (pubmed ids)
0252             if numel(tline)>12 && strcmp(tline(1:12),'CITATIONS - ')
0253                 if isempty(metaCycRxns.rxnReferences{rxnCounter})
0254                     metaCycRxns.rxnReferences{rxnCounter}=strcat('pubmed/',tline(13:end));
0255                 else
0256                     metaCycRxns.rxnReferences{rxnCounter}=strcat(metaCycRxns.rxnReferences{rxnCounter},';pubmed/',tline(13:end));
0257                 end
0258             end
0259             
0260             %Add Miriam info (cross-links other database)
0261             if numel(tline)>11 && strcmp(tline(1:11),'DBLINKS - (')
0262                 dblink=tline(12:end);
0263                 %KEGG reaction id
0264                 if strcmp(dblink(1:12),'LIGAND-RXN "')
0265                     dblink=dblink(13:end);
0266                     s=strfind(dblink,'"');
0267                     if any(s)
0268                         dblink=dblink(1:s-1);
0269                     end
0270                     
0271                     if isstruct(metaCycRxns.rxnMiriams{rxnCounter})
0272                         addToIndex=numel(metaCycRxns.rxnMiriams{rxnCounter}.name)+1;
0273                     else
0274                         addToIndex=1;
0275                     end
0276                     tempStruct=metaCycRxns.rxnMiriams{rxnCounter};
0277                     tempStruct.name{addToIndex,1}='kegg.reaction';
0278                     tempStruct.value{addToIndex,1}=dblink;
0279                     metaCycRxns.rxnMiriams{rxnCounter}=tempStruct;
0280                     
0281                     %For generating the rxnLinks structure
0282                     dbLinkCounter=dbLinkCounter+1;
0283                     rxnLinks.metacyc{dbLinkCounter}=metaCycRxns.rxns{rxnCounter};
0284                     rxnLinks.kegg{dbLinkCounter}=dblink;
0285                     rxnLinks.check{dbLinkCounter}=strcat(metaCycRxns.rxns{rxnCounter},dblink);
0286                 end
0287                 
0288                 %RHEA reaction id
0289                 if strcmp(dblink(1:6),'RHEA "')
0290                     dblink=dblink(7:end);
0291                     s=strfind(dblink,'"');
0292                     if any(s)
0293                         dblink=dblink(1:s-1);
0294                     end
0295                     
0296                     if isstruct(metaCycRxns.rxnMiriams{rxnCounter})
0297                         addToIndex=numel(metaCycRxns.rxnMiriams{rxnCounter}.name)+1;
0298                     else
0299                         addToIndex=1;
0300                     end
0301                     tempStruct=metaCycRxns.rxnMiriams{rxnCounter};
0302                     tempStruct.name{addToIndex,1}='rhea';
0303                     tempStruct.value{addToIndex,1}=dblink;
0304                     metaCycRxns.rxnMiriams{rxnCounter}=tempStruct;
0305                 end
0306             end
0307             
0308             if numel(tline)>21 && strcmp(tline(1:21),'REACTION-DIRECTION - ')
0309                 rxnDirection=tline(22:end);
0310                 switch(rxnDirection)
0311                     case 'IRREVERSIBLE-LEFT-TO-RIGHT'
0312                         metaCycRxns.rev(rxnCounter,1)=0;
0313                     case 'LEFT-TO-RIGHT'
0314                         metaCycRxns.rev(rxnCounter,1)=0;
0315                     case 'PHYSIOL-LEFT-TO-RIGHT'
0316                         metaCycRxns.rev(rxnCounter,1)=0;
0317                     case 'IRREVERSIBLE-RIGHT-TO-LEFT'
0318                         metaCycRxns.rev(rxnCounter,1)=0;
0319                         reverse=1;
0320                     case 'RIGHT-TO-LEFT'
0321                         metaCycRxns.rev(rxnCounter,1)=0;
0322                         reverse=1;
0323                     case 'PHYSIOL-RIGHT-TO-LEFT'
0324                         metaCycRxns.rev(rxnCounter,1)=0;
0325                         reverse=1;
0326                 end
0327             end
0328             
0329             %Tag transport reactions
0330             if strcmp(tline,'TYPES - Transport-Reactions')
0331                 TRANSPORT=[TRANSPORT;metaCycRxns.rxns{rxnCounter}];
0332             end
0333             
0334             %Add spontaneous
0335             if strcmp(tline,'SPONTANEOUS? - T')
0336                 %metaCycRxns.spontaneous(rxnCounter,1)=1;
0337                 isSpontaneous(rxnCounter)=true;
0338             end
0339             
0340             %Extract mass-balance status
0341             if numel(tline)>27 && strcmp(tline(1:27),'REACTION-BALANCE-STATUS - :')
0342                 if isequal(tline(28:35), 'UNBALANC')
0343                     UNBALANCED(rxnCounter)=true;
0344                 elseif isequal(tline(28:35), 'UNDETERM')
0345                     UNDETERMINED(rxnCounter)=true;
0346                 end
0347             end
0348             
0349             %Add left side equation
0350             if numel(tline)>7 && strcmp(tline(1:7),'LEFT - ')
0351                 if strcmp(left{rxnCounter},'')
0352                     if strcmp(templeft,'')
0353                         templeft=tline(8:end); % this is the real first, save to a temp variable
0354                     else % count in the coefficient here
0355                         if strcmp(coefficient,'')
0356                             left{rxnCounter}=templeft;
0357                         else
0358                             left{rxnCounter}=strcat(coefficient,32,templeft);
0359                             coefficient='';
0360                         end
0361                         templeft=tline(8:end);
0362                     end
0363                 else
0364                     if strcmp(coefficient,'')
0365                         left{rxnCounter}=strcat(left{rxnCounter},' +',32,templeft);
0366                     else
0367                         left{rxnCounter}=strcat(left{rxnCounter},' +',32,coefficient,32,templeft);
0368                         coefficient='';
0369                     end
0370                     templeft=tline(8:end);
0371                 end
0372             end
0373             
0374             %Add right side equation
0375             if numel(tline)>8 && strcmp(tline(1:8),'RIGHT - ')
0376                 if strcmp(right{rxnCounter},'')
0377                     if strcmp(tempright,'')
0378                         
0379                         % a complicated process for the last left
0380                         % metabolite
0381                         if strcmp(coefficient,'')
0382                             if strcmp(left{rxnCounter},'')
0383                                 left{rxnCounter}=templeft;
0384                             else
0385                                 left{rxnCounter}=strcat(left{rxnCounter},' +',32,templeft);
0386                             end
0387                             
0388                         else
0389                             if strcmp(left{rxnCounter},'')
0390                                 left{rxnCounter}=strcat(coefficient,32,templeft);
0391                             else
0392                                 left{rxnCounter}=strcat(left{rxnCounter},' +',32,coefficient,32,templeft);
0393                             end
0394                             coefficient='';
0395                         end
0396                         % process end
0397                         
0398                         tempright=tline(9:end); % this is the real first, save to a temp variable
0399                     else
0400                         if strcmp(coefficient,'')
0401                             right{rxnCounter}=tempright;
0402                         else
0403                             right{rxnCounter}=strcat(coefficient,32,tempright);
0404                             coefficient='';
0405                         end
0406                         tempright=tline(9:end);
0407                     end
0408                 else
0409                     if strcmp(coefficient,'')
0410                         right{rxnCounter}=strcat(right{rxnCounter},' +',32,tempright);
0411                     else
0412                         right{rxnCounter}=strcat(right{rxnCounter},' +',32,coefficient,32,tempright);
0413                         coefficient='';
0414                     end
0415                     tempright=tline(9:end);
0416                 end
0417                 
0418             end
0419             
0420             if numel(tline)>15 && strcmp(tline(1:15),'^COEFFICIENT - ')
0421                 coefficient=tline(16:end);
0422             end
0423             
0424             %Generate equation at the end of each object section
0425             if strcmp(tline,'//')
0426                 
0427                 % a complicated process for the last right metabolite,
0428                 % sub-function is needed here
0429                 if strcmp(coefficient,'')
0430                     if strcmp(right{rxnCounter},'')
0431                         right{rxnCounter}=tempright;
0432                     else
0433                         right{rxnCounter}=strcat(right{rxnCounter},' +',32,tempright);
0434                     end
0435                     
0436                 else
0437                     if strcmp(right{rxnCounter},'')
0438                         right{rxnCounter}=strcat(coefficient,32,tempright);
0439                     else
0440                         right{rxnCounter}=strcat(right{rxnCounter},' +',32,coefficient,32,tempright);
0441                     end
0442                     coefficient='';
0443                 end
0444                 % process end
0445                 
0446                 %get the right direction symbol
0447                 if metaCycRxns.rev(rxnCounter,1)
0448                     symbol = ' <=>';
0449                 else
0450                     symbol = ' =>';
0451                 end
0452                 
0453                 if reverse
0454                     metaCycRxns.equations{rxnCounter}=strcat(right{rxnCounter},symbol,32,left{rxnCounter});
0455                 else
0456                     metaCycRxns.equations{rxnCounter}=strcat(left{rxnCounter},symbol,32,right{rxnCounter});
0457                 end
0458                 
0459                 %Final equation check
0460                 if strcmp(left{rxnCounter},'') || strcmp(right{rxnCounter},'')
0461                     rxnCounter=rxnCounter-1;
0462                     
0463                 end
0464                 
0465             end
0466             
0467         end
0468         %Close the file
0469         fclose(fid);
0470         
0471         %===
0472         UNBALANCED=metaCycRxns.rxns(UNBALANCED);
0473         UNDETERMINED=metaCycRxns.rxns(UNDETERMINED);
0474         isSpontaneous=metaCycRxns.rxns(isSpontaneous);
0475         
0476         %If too much space was allocated, shrink the model
0477         metaCycRxns.rxns=metaCycRxns.rxns(1:rxnCounter);
0478         metaCycRxns.rxnNames=metaCycRxns.rxnNames(1:rxnCounter);
0479         metaCycRxns.eccodes=metaCycRxns.eccodes(1:rxnCounter);
0480         metaCycRxns.equations=metaCycRxns.equations(1:rxnCounter);
0481         metaCycRxns.rxnMiriams=metaCycRxns.rxnMiriams(1:rxnCounter);
0482         metaCycRxns.rxnReferences=metaCycRxns.rxnReferences(1:rxnCounter);
0483         metaCycRxns.subSystems=metaCycRxns.subSystems(1:rxnCounter);
0484         metaCycRxns.pwys=metaCycRxns.pwys(1:rxnCounter);
0485         metaCycRxns.rev=metaCycRxns.rev(1:rxnCounter,:);
0486         
0487         rxnLinks.kegg=rxnLinks.kegg(1:dbLinkCounter);
0488         rxnLinks.metacyc=rxnLinks.metacyc(1:dbLinkCounter);
0489         rxnLinks.check=rxnLinks.check(1:dbLinkCounter);
0490         [~,index]=unique(rxnLinks.check);
0491         rxnLinks.kegg=rxnLinks.kegg(index);
0492         rxnLinks.metacyc=rxnLinks.metacyc(index);
0493         rxnLinks=rmfield(rxnLinks,'check');
0494         
0495         %Construct the S matrix and list of metabolites
0496         [S, mets, badRxns]=constructS(metaCycRxns.equations);
0497         metaCycRxns.S=S;
0498         metaCycRxns.mets=mets;
0499         
0500         %Add some stuff to get a correct model structure
0501         metaCycRxns.ub=ones(rxnCounter,1)*1000;
0502         metaCycRxns.lb=metaCycRxns.rev*-1000;
0503         metaCycRxns.c=zeros(rxnCounter,1);
0504         metaCycRxns.b=zeros(numel(metaCycRxns.mets),1);
0505         metaCycRxns.version=version;
0506         
0507         %Save the model structure
0508         save(rxnsFile,'metaCycRxns','rxnLinks','TRANSPORT','UNBALANCED','UNDETERMINED','isSpontaneous');
0509         fprintf(['New metaCycRxns.mat has been successfully updated!\n\n']);
0510     end
0511 end
0512 
0513 %Deal with reactions that are labeled as "TRANSPORT", "UNBALANCED", or
0514 %"UNDETERMINED" (depending on settings).
0515 model=metaCycRxns;
0516 if keepTransportRxns==false
0517     model=removeReactions(model,intersect(TRANSPORT,model.rxns),true,true);
0518 end
0519 if keepUnbalanced==false
0520     model=removeReactions(model,intersect(UNBALANCED,model.rxns),true,true);
0521 end
0522 if keepUndetermined==false
0523     model=removeReactions(model,intersect(UNDETERMINED,model.rxns),true,true);
0524 end
0525 end
0526 
0527 %Sub function for romving HTML symbols from the names of reaction and
0528 %pathway
0529 function newString=removeHTMLcodes(string)
0530 string=regexprep(string,'<(\w+)>','');
0531 string=regexprep(string,'</(\w+)>','');
0532 string=regexprep(string,'[&;]','');
0533 newString=string;
0534 end