linkMetaCycKEGGRxns

 linkMetaCycKEGGRxns
   Link additional MetaCyc and KEGG reactions through metabolite mapping
   This function only need to run once when the MetaCyc database is updated

   NOTE: No arguments are required

   Usage: linkMetaCycKEGGRxns

0001 function linkMetaCycKEGGRxns
0002 % linkMetaCycKEGGRxns
0003 %   Link additional MetaCyc and KEGG reactions through metabolite mapping
0004 %   This function only need to run once when the MetaCyc database is updated
0005 %
0006 %   NOTE: No arguments are required
0007 %
0008 %   Usage: linkMetaCycKEGGRxns
0009 
0010 load('metaCycRxns.mat'); %load MetaCyc reactions
0011 fprintf('NOTE: Importing MetaCyc reactions...\n');
0012 metaCycModel = metaCycRxns;
0013 metaCycModel.rxnFrom=cell(numel(metaCycModel.rxns),1);
0014 metaCycModel.rxnFrom(:)={'MetaCyc'};
0015 metaCycModel.grRules={};
0016 metaCycModel.genes={};
0017 metaCycModel.rxnGeneMat=sparse(numel(metaCycModel.rxns),1);
0018 
0019 
0020 keggModel=getRxnsFromKEGG(); %load KEGG reactions
0021 
0022 %Remove KEGG reactions with MetaCyc links based on original rxnLinks
0023 rxnToRemove=intersect(rxnLinks.kegg,keggModel.rxns);
0024 rxnToRemove=unique(rxnToRemove);
0025 shrinkedKeggModel=removeReactions(keggModel,rxnToRemove,true,true);
0026 
0027 fprintf('Mapping MetaCyc and KEGG reactions...\n');
0028 %Resolve the shared but unmapped reactions (through mapping the involved
0029 %metabolites). Replace mets information in KEGG model with the
0030 %corresponding ones in MetaCyc.
0031 load('metaCycMets.mat');
0032 for i=1:numel(shrinkedKeggModel.mets)
0033     [a, b]=ismember(shrinkedKeggModel.mets{i},metaCycMets.keggid);
0034     if a
0035         %Replace the met ids
0036         shrinkedKeggModel.mets{i}=metaCycMets.mets{b};
0037         shrinkedKeggModel.metNames{i}=metaCycMets.mets{b};
0038     end
0039 end
0040 
0041 %Prepare for the merge of KEGG and MetaCyc super models
0042 
0043 %Adding fields (comps, compNames, metNames, metComps)
0044 metaCycModel.comps={'s'};
0045 metaCycModel.compNames={'System'};
0046 metaCycModel.metNames=metaCycModel.mets;
0047 if ~isfield(metaCycModel,'metComps')
0048     metaCycModel.metComps=ones(numel(metaCycModel.mets),1);
0049 end
0050 
0051 shrinkedKeggModel.comps={'s'};
0052 shrinkedKeggModel.compNames={'System'};
0053 shrinkedKeggModel.metNames=shrinkedKeggModel.mets;
0054 if ~isfield(shrinkedKeggModel,'metComps')
0055     shrinkedKeggModel.metComps=ones(numel(shrinkedKeggModel.mets),1);
0056 end
0057 
0058 %Merge models
0059 mappingModel=mergeModels({shrinkedKeggModel metaCycModel},'metNames');
0060 
0061 %Remove compounds proton and water because KEGG reactions often miss them
0062 mappingModel=removeMets(mappingModel,{'PROTON','WATER'});
0063 
0064 %Find same/similiar reactions by check identical columns in S matrix Go
0065 %through the pure KEGG model and collect rxns id that have identical
0066 %reactions to the ones in MetaCyc
0067 I=find(ismember(mappingModel.rxnFrom,'KEGG'));
0068 for i=1:numel(I)
0069     testRxn=mappingModel.S(:,I(i));
0070     match=find(~any(bsxfun(@minus, mappingModel.S,testRxn))); %Figure out the matched columns in S matrix
0071     % if ~isequal(match,I(i))
0072     count=numel(match);
0073     if count>1 %Only consider one-to-one match here
0074         for j=2:count
0075             %Only consider the match between KEGG and MetaCyc
0076             if isequal(mappingModel.rxnFrom{match(j)},'MetaCyc')
0077                 %disp([mappingModel.rxns(match(1))
0078                 %mappingModel.rxns(match(j))]); %keep this for later
0079                 %checking
0080                 rxnLinks.kegg=[rxnLinks.kegg;mappingModel.rxns(match(1))];
0081                 rxnLinks.metacyc=[rxnLinks.metacyc;mappingModel.rxns(match(j))];
0082             end
0083         end
0084     end
0085     % end
0086 end
0087 
0088 % Load orginal rxnLinks
0089 numLink=numel(rxnLinks.kegg);
0090 rxnLinks.check=cell(numLink,1);
0091 for i=1:numLink
0092     rxnLinks.check{i}=strcat(rxnLinks.kegg{i},rxnLinks.metacyc{i});
0093 end
0094 
0095 [~, index]=unique(rxnLinks.check);
0096 rxnLinks.kegg=rxnLinks.kegg(index);
0097 rxnLinks.metacyc=rxnLinks.metacyc(index);
0098 rxnLinks=rmfield(rxnLinks,'check');
0099 
0100 %Get the MetaCyc path and update the metaCycRxns.mat
0101 ravenPath=findRAVENroot();
0102 rxnsFile=fullfile(ravenPath,'external','metacyc','metaCycRxns.mat');
0103 save(rxnsFile,'metaCycRxns','rxnLinks','TRANSPORT','UNBALANCED','UNDETERMINED','isSpontaneous');
0104 fprintf(['Reaction associations between MetaCyc and KEGG have been successfully updated!\n\n']);
0105 
0106 end