Description of getRxnsFromKEGG

0001 function [model,isSpontaneous,isUndefinedStoich,isIncomplete,...
0002     isGeneral]=getRxnsFromKEGG(keggPath)
0003 % getRxnsFromKEGG
0004 %   Retrieves information on all reactions stored in KEGG database
0005 %
0006 %   Input:
0007 %   keggPath            if keggRxns.mat is not in the RAVEN\external\kegg
0008 %                       directory, this function will attempt to read data
0009 %                       from a local FTP dump of the KEGG database.
0010 %                       keggPath is the path to the root of this database
0011 %
0012 %   Output:
0013 %   model               a model structure generated from the database. The
0014 %                       following fields are filled
0015 %       id                  'KEGG'
0016 %       name         'Automatically generated from KEGG database'
0017 %       rxns                KEGG reaction ids
0018 %       rxnNames            Name for each reaction entry
0019 %       mets                KEGG compound ids. If the equations use
0020 %                           stoichiometry such as ID(n+1) then the whole
0021 %                           expression is saved as the id
0022 %       eccodes             Corresponding ec-number if available
0023 %       rxnMiriams          Contains reaction specific information such as
0024 %                           KO id and pathways that the reaction is
0025 %                           associated to
0026 %       S                   Stoichiometric matrix
0027 %       lb                  -1000 for all reactions
0028 %       ub                  1000 for all reactions
0029 %       rev                 1 for reversible and 0 for irreversible. For
0030 %                           reactions present in pathway maps the
0031 %                           reversibility is taken from there
0032 %       b                   0 for all metabolites
0033 %   isSpontaneous       a cell array with the reactions labelled as
0034 %                       "spontaneous"
0035 %   isUndefinedStoich    a cell array with the reactions labelled as with
0036 %                       undefined stoichiometry
0037 %   isIncomplete        a cell array with the reactions labelled as
0038 %                       "incomplete", "erroneous" or "unclear"
0039 %   isGeneral           a cell array with the reactions labelled as
0040 %                       "general reaction"
0041 %
0042 %   NOTE: Reactions on the form A <=> A + B will not be loaded. If the file
0043 %   keggRxns.mat is in the RAVEN/external/kegg directory it will be loaded
0044 %   instead of parsing of the KEGG files. If it does not exist it will be
0045 %   saved after parsing of the KEGG files. In general, you should remove
0046 %   the keggRxns.mat file along with other KEGG mat files if you want to
0047 %   rebuild the model structure from a newer version of KEGG.
0048 %
0049 % Usage: [model,isSpontaneous,isUndefinedStoich,isIncomplete,...
0050 %    isGeneral]=getRxnsFromKEGG(keggPath)
0051 %
0052 % NOTE: This is how one entry looks in the file
0053 %
0054 % ENTRY       R00010                      Reaction
0055 % NAME        alpha,alpha-trehalose glucohydrolase
0056 % DEFINITION  alpha,alpha-Trehalose + H2O <=> 2 D-Glucose
0057 % EQUATION    C01083 + C00001 <=> 2 C00031
0058 % REMARK      Same as: R06103
0059 % RCLASS      RC00049  C00031_C01083
0060 % ENZYME      3.2.1.28
0061 % PATHWAY     rn00500  Starch and sucrose metabolism
0062 %             rn01100  Metabolic pathways
0063 % ORTHOLOGY   K01194  alpha,alpha-trehalase [EC:3.2.1.28]
0064 % DBLINKS     RHEA: 32678
0065 % ///
0066 %
0067 % The file is not tab-delimited. Instead each label is 12 characters
0068 % (except for '///')
0069 %
0070 
0071 ravenPath=findRAVENroot();
0072 
0073 if nargin<1
0074     keggPath=fullfile(ravenPath,'external','kegg');
0075 else
0076     keggPath=char(keggPath);
0077 end
0078 
0079 %Check if the reactions have been parsed before and saved. If so, load the
0080 %model
0081 rxnsFile=fullfile(ravenPath,'external','kegg','keggRxns.mat');
0082 if exist(rxnsFile, 'file')
0083     fprintf(['Importing KEGG reactions from ' strrep(rxnsFile,'\','/') '... ']);
0084     load(rxnsFile);
0085 else
0086     fprintf(['NOTE: Cannot locate ' strrep(rxnsFile,'\','/') ', it will therefore be generated from the local KEGG database\n']);
0087     if ~isfile(fullfile(keggPath,'reaction')) || ~isfile(fullfile(keggPath,'reaction.lst')) || ~isfile(fullfile(keggPath,'reaction_mapformula.lst'))
0088         EM=fprintf(['The files ''reaction'', ''reaction.lst'' and ''reaction_mapformula.lst'' cannot be located at ' strrep(keggPath,'\','/') '/ and should be downloaded from the KEGG FTP\n']);
0089         dispEM(EM);
0090     else
0091         fprintf('Generating keggRxns.mat file... ');
0092         %Add new functionality in the order specified in models
0093         model.id='KEGG';
0094         model.name='Automatically generated from KEGG database';
0095         
0096         %Preallocate memory for 15000 reactions
0097         model.rxns=cell(15000,1);
0098         model.rxnNames=cell(15000,1);
0099         model.eccodes=cell(15000,1);
0100         model.subSystems=cell(15000,1);
0101         model.rxnMiriams=cell(15000,1);
0102         model.rxnNotes=cell(15000,1);
0103         equations=cell(15000,1);
0104         %Temporarily store the equations
0105         
0106         isSpontaneous=false(15000,1);
0107         isIncomplete=false(15000,1);
0108         isGeneral=false(15000,1);
0109 
0110         %First load information on reaction ID, reaction name, KO, pathway,
0111         %and ec-number
0112         fid = fopen(fullfile(keggPath,'reaction'), 'r');
0113         
0114         %Keep track of how many reactions have been added
0115         rxnCounter=0;
0116         
0117         %Loop through the file
0118         orthology=false;
0119         pathway=false;
0120         module=false;
0121         while 1
0122             %Get the next line
0123             tline = fgetl(fid);
0124             
0125             %Abort at end of file
0126             if ~ischar(tline)
0127                 break;
0128             end
0129             
0130             %Skip '///'
0131             if numel(tline)<12
0132                 continue;
0133             end
0134             % Skip other lines with unused information
0135             if strcmp(tline(1:5),'BRITE')
0136                 pathway = false;
0137                 continue;
0138             end
0139 
0140             %Check if it's a new reaction
0141             if strcmp(tline(1:12),'ENTRY       ')
0142                 rxnCounter=rxnCounter+1;
0143                 
0144                 %Add empty strings where there should be such
0145                 model.rxnNames{rxnCounter}='';
0146                 model.eccodes{rxnCounter}='';
0147                 %model.subSystems{rxnCounter}={''}; %remain empty cell
0148                 model.rxnNotes{rxnCounter}='';
0149                 equations{rxnCounter}='';
0150                 
0151                 %Add reaction ID (always 6 characters)
0152                 model.rxns{rxnCounter}=tline(13:18);
0153                 orthology=false;
0154                 pathway=false;
0155                 module=false;
0156                 
0157                 %Add KEGG reaction ID miriam
0158                 tempStruct=model.rxnMiriams{rxnCounter};
0159                 tempStruct.name{1,1}='kegg.reaction';
0160                 tempStruct.value{1,1}=tline(13:18);
0161                 model.rxnMiriams{rxnCounter}=tempStruct;
0162             end
0163             
0164             %Add name
0165             if strcmp(tline(1:12),'NAME        ')
0166                 model.rxnNames{rxnCounter}=tline(13:end);
0167             end
0168             
0169             %Add whether the comment includes "incomplete", "erroneous" or
0170             %"unclear"
0171             if strcmp(tline(1:12),'COMMENT     ')
0172                 %Read all text until '///', 'RPAIR', 'ENZYME', 'PATHWAY' or 'RCLASS'
0173                 commentText=tline(13:end);
0174                 while 1
0175                     tline = fgetl(fid);
0176                     if ~strcmp(tline(1:3),'///') && ~strcmp(tline(1:3),'RPA') && ~strcmp(tline(1:3),'ENZ') && ~strcmp(tline(1:3),'PAT') && ~strcmp(tline(1:3),'RCL')
0177                         commentText=[commentText ' ' strtrim(tline)];
0178                     else
0179                         break;
0180                     end
0181                 end
0182                 if any(regexpi(commentText,'SPONTANEOUS'))
0183                     %It should start this way
0184                     isSpontaneous(rxnCounter)=true;
0185                 end
0186                 if any(regexpi(commentText,'INCOMPLETE')) || any(regexpi(commentText,'ERRONEOUS')) || any(regexpi(commentText,'UNCLEAR'))
0187                     isIncomplete(rxnCounter)=true;
0188                 end
0189                 if any(regexpi(commentText,'GENERAL REACTION'))
0190                     %It should start this way
0191                     isGeneral(rxnCounter)=true;
0192                 end
0193                 
0194                 %Go to next iteration if it is '///'
0195                 if numel(tline)<12
0196                     continue;
0197                 end
0198             end
0199             
0200             %Add ec-number
0201             if strcmp(tline(1:12),'ENZYME      ')
0202                 model.eccodes{rxnCounter}=tline(13:end);
0203                 model.eccodes{rxnCounter}=deblank(model.eccodes{rxnCounter});
0204                 model.eccodes{rxnCounter}=regexprep(model.eccodes{rxnCounter},'\s+',';');
0205             end
0206             if numel(tline)>8
0207                 if strcmp(tline(1:9),'REFERENCE')
0208                     pathway=false;
0209                     orthology=false;
0210                 end
0211             end
0212             
0213             %Add module ids
0214             if numel(tline)>18
0215                 if strcmp(tline(1:12),'MODULE      ') || module==true
0216                     pathway=false;
0217                     orthology=false;
0218                     if isstruct(model.rxnMiriams{rxnCounter})
0219                         addToIndex=numel(model.rxnMiriams{rxnCounter}.name)+1;
0220                     else
0221                         addToIndex=1;
0222                     end
0223                     tempStruct=model.rxnMiriams{rxnCounter};
0224                     tempStruct.name{addToIndex,1}='kegg.module';
0225                     tempStruct.value{addToIndex,1}=tline(13:18);
0226                     model.rxnMiriams{rxnCounter}=tempStruct;
0227                 end
0228             end
0229             
0230             %Add RHEA id
0231             if numel(tline)>18
0232                 if strcmp(tline(1:18),'DBLINKS     RHEA: ')
0233                     pathway=false;
0234                     orthology=false;
0235                     module=false;
0236                     if isstruct(model.rxnMiriams{rxnCounter})
0237                         addToIndex=numel(model.rxnMiriams{rxnCounter}.name)+1;
0238                     else
0239                         addToIndex=1;
0240                     end
0241                     tempStruct=model.rxnMiriams{rxnCounter};
0242                     tempStruct.name{addToIndex,1}='rhea';
0243                     tempStruct.value{addToIndex,1}=tline(19:end);
0244                     model.rxnMiriams{rxnCounter}=tempStruct;
0245                 end
0246             end
0247             
0248             %Add KO-ids
0249             if numel(tline)>16
0250                 if strcmp(tline(1:16),'ORTHOLOGY   KO: ') || strcmp(tline(1:16),'            KO: ') || strcmp(tline(1:12),'ORTHOLOGY   ') || orthology==true
0251                     pathway=false;
0252                     module=false;
0253                     %Check if KO has been added already (each reaction may
0254                     %belong to several)
0255                     if isstruct(model.rxnMiriams{rxnCounter})
0256                         addToIndex=numel(model.rxnMiriams{rxnCounter}.name)+1;
0257                     else
0258                         addToIndex=1;
0259                     end
0260                     
0261                     tempStruct=model.rxnMiriams{rxnCounter};
0262                     tempStruct.name{addToIndex,1}='kegg.orthology';
0263                     if strcmp(tline(13:16),'KO:')
0264                         %This is in the old version
0265                         tempStruct.value{addToIndex,1}=tline(17:22);
0266                     else
0267                         %This means that it found one KO in the new format
0268                         %and that subsequent lines might be other KOs
0269                         orthology=true;
0270                         tempStruct.value{addToIndex,1}=tline(13:18);
0271                     end
0272                     model.rxnMiriams{rxnCounter}=tempStruct;
0273                 end
0274             end
0275             
0276             %Add pathways
0277             if numel(tline)>18
0278                 if strcmp(tline(1:18),'PATHWAY     PATH: ') || strcmp(tline(1:18),'            PATH: ') || strcmp(tline(1:12),'PATHWAY     ') || pathway==true
0279                     orthology=false;
0280                     module=false;
0281                     %Check if annotation has been added already
0282                     if isstruct(model.rxnMiriams{rxnCounter})
0283                         addToIndex=numel(model.rxnMiriams{rxnCounter}.name)+1;
0284                     else
0285                         addToIndex=1;
0286                     end
0287                     
0288                     tempStruct=model.rxnMiriams{rxnCounter};
0289                     tempStruct.name{addToIndex,1}='kegg.pathway';
0290                     %If it is the old version
0291                     if strcmp(tline(14:17),'PATH:')
0292                         tempStruct.value{addToIndex,1}=tline(19:25);
0293                     else
0294                         %If it is the new version
0295                         tempStruct.value{addToIndex,1}=tline(13:19);
0296                         pathway=true;
0297                     end
0298                     
0299                     %Do not save global or overview pathways. The ids for
0300                     %such pathways begin with rn011 or rn012
0301                     if ~strcmp('rn011',tempStruct.value{addToIndex,1}(1:5)) && ~strcmp('rn012',tempStruct.value{addToIndex,1}(1:5))
0302                         model.rxnMiriams{rxnCounter}=tempStruct;
0303                         
0304                         %Also save the subSystems names. For the old KEGG
0305                         %format, only the first mentioned subsystem is
0306                         %picked. Use the newer KEGG format to fetch all the
0307                         %subsystems
0308                         if strcmp(tline(14:17),'PATH:')
0309                             %The old format
0310                             model.subSystems{rxnCounter}=tline(28:end);
0311                         else
0312                             %The new format
0313                             model.subSystems{rxnCounter,1}{numel(model.subSystems{rxnCounter,1})+1,1}=tline(22:end);
0314                         end
0315                     end
0316                 end
0317             end
0318         end
0319         
0320         %Close the file
0321         fclose(fid);
0322         
0323         %This is done here since the the indexes won't match since some
0324         %reactions are removed along the way
0325         isIncomplete=model.rxns(isIncomplete);
0326         isGeneral=model.rxns(isGeneral);
0327         isSpontaneous=model.rxns(isSpontaneous);
0328 
0329         %If too much space was allocated, shrink the model
0330         model.rxns=model.rxns(1:rxnCounter);
0331         model.rxnNames=model.rxnNames(1:rxnCounter);
0332         model.eccodes=model.eccodes(1:rxnCounter);
0333         equations=equations(1:rxnCounter);
0334         model.rxnMiriams=model.rxnMiriams(1:rxnCounter);
0335         model.rxnNotes=model.rxnNotes(1:rxnCounter);
0336         model.subSystems=model.subSystems(1:rxnCounter);
0337         
0338         emptySubSys = cellfun(@isempty,model.subSystems);
0339         model.subSystems(emptySubSys) = {{''}};
0340 
0341         %Then load the equations from another file. This is because the
0342         %equations are easier to retrieve from there
0343         
0344         %The format is rxnID: equation The reactions should have been
0345         %loaded in the exact same order
0346         fid = fopen(fullfile(keggPath,'reaction.lst'), 'r');
0347         
0348         %Loop through the file
0349         for i=1:rxnCounter
0350             %Get the next line
0351             tline = fgetl(fid);
0352             
0353             equations{i}=tline(9:end);
0354         end
0355         
0356         %Close the file
0357         fclose(fid);
0358         
0359         %Several equations may have two whitespaces between the last
0360         %reactant and the reversible arrow sign. The number of whitespaces
0361         %is thus reduced to one
0362         equations = regexprep(equations,'  <=>', ' <=>');
0363         
0364         %Construct the S matrix and list of metabolites
0365         [S, mets, badRxns]=constructS(equations);
0366         model.S=S;
0367         model.mets=mets;
0368         
0369         %There is some limited evidence for directionality in
0370         %reaction_mapformula.lst. The information there concerns how the
0371         %reactions are drawn in the KEGG maps. If a reaction is
0372         %irreversible in the same direction for all maps, then I consider
0373         %is irreversible, otherwise reversible. Also, not all reactions are
0374         %present in the maps, so not all will have directionality. They
0375         %will be considered to be reversible
0376         
0377         %The format is R00005: 00330: C01010 => C00011 Generate a
0378         %reversibility structure with the fields: *rxns: reaction ids
0379         %*product: one met id that is a product. This is because the
0380         %*reactions might be written in another direction compared to in
0381         % the reactions.lst file
0382         %*rev: 1 if reversible, otherwise 0
0383         reversibility.rxns={};
0384         reversibility.product={};
0385         reversibility.rev=[];
0386         
0387         fid = fopen(fullfile(keggPath,'reaction_mapformula.lst'), 'r');
0388         while 1
0389             %Get the next line
0390             tline = fgetl(fid);
0391             
0392             %Abort at end of file
0393             if ~ischar(tline)
0394                 break;
0395             end
0396             
0397             rxn=tline(1:6);
0398             prod=tline(end-5:end);
0399             rev=any(strfind(tline,'<=>'));
0400             if isempty(reversibility.rxns)
0401                 reversibility.rxns{1}=rxn;
0402                 reversibility.product{1}=prod;
0403                 reversibility.rev(1)=rev;
0404             elseif strcmp(reversibility.rxns(end),rxn)
0405                 %Check if the reaction was added before. It's an ordered
0406                 %list, so only check the last element If it's reversible in
0407                 %the new reaction or reversible in the old reaction then
0408                 %set (keep) to be reversible
0409                 if rev==1 || reversibility.rev(end)==1
0410                     reversibility.rev(end)=1;
0411                 else
0412                     %This means that the reaction was already loaded, that
0413                     %it was irreversible before and irreversible in the new
0414                     %reaction. However, it could be that they are written
0415                     %in diferent directions. If the product differ, then
0416                     %set to be reversible. This assumes that the reactions
0417                     %are written with the same metabolite as the last one
0418                     %if they are in the same direction
0419                     if ~strcmp(prod,reversibility.product(end))
0420                         reversibility.rev(end)=1;
0421                     end
0422                 end
0423             else
0424                 reversibility.rxns=[reversibility.rxns;rxn];
0425                 reversibility.product=[reversibility.product;prod];
0426                 reversibility.rev=[reversibility.rev;rev];
0427             end
0428         end
0429         fclose(fid);
0430         
0431         %Update the reversibility
0432         model.rev=ones(rxnCounter,1);
0433         %Match the reaction ids
0434         irrevIDs=find(reversibility.rev==0);
0435         [~, I]=ismember(reversibility.rxns(irrevIDs),model.rxns);
0436         [~, prodMetIDs]=ismember(reversibility.product(irrevIDs),model.mets);
0437         model.rev(I)=0;
0438         
0439         %See if the reactions are written in the same order in model.S
0440         linearInd=sub2ind(size(model.S), prodMetIDs, I);
0441         changeOrder=I(model.S(linearInd)<0);
0442         %Change the order of these reactions
0443         model.S(:,changeOrder)=model.S(:,changeOrder).*-1;
0444         
0445         %Add some stuff to get a correct model structure
0446         model.ub=ones(rxnCounter,1)*1000;
0447         model.lb=model.rev*-1000;
0448         model.c=zeros(rxnCounter,1);
0449         model.b=zeros(numel(model.mets),1);
0450         model=removeReactions(model,badRxns,true,true);
0451         
0452         %Identify reactions with undefined stoichiometry. Such
0453         %reactions involve metabolites with an ID containing the letter "n"
0454         %or "m"
0455         I=cellfun(@any,strfind(model.mets,'n')) | cellfun(@any,strfind(model.mets,'m'));
0456         [~, J]=find(model.S(I,:));
0457         isUndefinedStoich=model.rxns(unique(J));   
0458         %Sort model that metabolites with undefined stoichiometry would
0459         %appear in the end of metabolites list
0460         metList=[model.mets(~I);model.mets(I)];
0461         [~,metIndexes]=ismember(metList,model.mets);
0462         model=permuteModel(model,metIndexes,'mets');
0463         
0464         %Sort model that i) spontaneous, ii) with undefined
0465         %stoichiometry, iii) incomplete and iv) general reactions would bve
0466         %ranked in the end of the model
0467         endRxnList=unique([model.rxns(ismember(model.rxns,isSpontaneous));model.rxns(ismember(model.rxns,isUndefinedStoich));model.rxns(ismember(model.rxns,isIncomplete));model.rxns(ismember(model.rxns,isGeneral))],'stable');
0468         rxnList=[model.rxns(~ismember(model.rxns,endRxnList));endRxnList];
0469         [~,rxnIndexes]=ismember(rxnList,model.rxns);
0470         model=permuteModel(model,rxnIndexes,'rxns');
0471         
0472         %Add information in rxnNotes, whether reaction belongs to any of
0473         %type i-iv mentioned a few lines above
0474         for i=(numel(rxnList)-numel(endRxnList)+1):numel(model.rxns)
0475             if ismember(model.rxns(i),isSpontaneous)
0476                 model.rxnNotes(i)=strcat(model.rxnNotes(i),'Spontaneous');
0477             end
0478             if ismember(model.rxns(i),isUndefinedStoich)
0479                 if isempty(model.rxnNotes{i})
0480                     model.rxnNotes(i)=strcat(model.rxnNotes(i),'With undefined stoichiometry');
0481                 else
0482                     model.rxnNotes(i)=strcat(model.rxnNotes(i),', with undefined stoichiometry');
0483                 end
0484             end
0485             if ismember(model.rxns(i),isIncomplete)
0486                 if isempty(model.rxnNotes{i})
0487                     model.rxnNotes(i)=strcat(model.rxnNotes(i),'Incomplete');
0488                 else
0489                     model.rxnNotes(i)=strcat(model.rxnNotes(i),', incomplete');
0490                 end
0491             end
0492             if ismember(model.rxns(i),isGeneral)
0493                 if isempty(model.rxnNotes{i})
0494                     model.rxnNotes(i)=strcat(model.rxnNotes(i),'General');
0495                 else
0496                     model.rxnNotes(i)=strcat(model.rxnNotes(i),', general');
0497                 end
0498             end
0499             model.rxnNotes(i)=strcat(model.rxnNotes(i),' reaction');
0500         end
0501         
0502         %Save the model structure
0503         save(rxnsFile,'model','isGeneral','isIncomplete','isUndefinedStoich','isSpontaneous');
0504     end
0505 end
0506 fprintf('COMPLETE\n');
0507 
0508 end
getRxnsFromKEGG

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
