addSpontaneousRxns

 addSpontaneousRxns
   Retrieve spontaneous reactions based on the pathway-spontaneousRxn
   associations curated by the MetaCyc database

   Input:
   rxnList              query list of reactions in cell array
   metList              query list of metabolites in cell array

   spontaneousRxnList   reterieved spontaneous reactions associated to
                        the queried MetaCyc reactions and metabolites

   Output:
   pathwayID            the cell array of relevant pathways

   Usage: spontaneousRxnList=addSpontaneousRxns(rxnList, metList)

0001 function [spontaneousRxnList, pathwayID]=addSpontaneousRxns(rxnList, metList)
0002 % addSpontaneousRxns
0003 %   Retrieve spontaneous reactions based on the pathway-spontaneousRxn
0004 %   associations curated by the MetaCyc database
0005 %
0006 %   Input:
0007 %   rxnList              query list of reactions in cell array
0008 %   metList              query list of metabolites in cell array
0009 %
0010 %   spontaneousRxnList   reterieved spontaneous reactions associated to
0011 %                        the queried MetaCyc reactions and metabolites
0012 %
0013 %   Output:
0014 %   pathwayID            the cell array of relevant pathways
0015 %
0016 %   Usage: spontaneousRxnList=addSpontaneousRxns(rxnList, metList)
0017 
0018 metList=convertCharArray(metList);
0019 rxnList=convertCharArray(rxnList);
0020 rxnList=unique(rxnList);
0021 
0022 % Create the matrix of MetaCyc pathways and spontaneous reactions
0023 load('metaCycRxns.mat');
0024 pathways={};
0025 for i=1:numel(metaCycRxns.rxns)
0026     if ~isempty(metaCycRxns.pwys{i})
0027         pathways=[pathways;transpose(strsplit(metaCycRxns.pwys{i},';'))];
0028     end
0029 end
0030 pathways=unique(pathways);
0031 
0032 % Genearte the matirx, row: spontaneousRxns, column: pathways
0033 sprxnPwyMat=zeros(numel(isSpontaneous),numel(pathways));
0034 for i=1:numel(isSpontaneous)
0035     [a, b]=ismember(isSpontaneous{i},metaCycRxns.rxns);
0036     if ~isempty(metaCycRxns.pwys{b})
0037         [crap, indexes]=ismember(strsplit(metaCycRxns.pwys{b},';'),pathways);
0038         sprxnPwyMat(i,indexes)=1;
0039     end
0040 end
0041 
0042 % Go through the rxnList to obtain relevant pathways
0043 pwys={};
0044 for i=1:numel(rxnList)
0045     [a, b]=ismember(rxnList{i},metaCycRxns.rxns);
0046     if a && ~isempty(metaCycRxns.pwys{b})
0047         pwys=[pwys;transpose(strsplit(metaCycRxns.pwys{b},';'))];
0048     end
0049 end
0050 pwys=unique(pwys);
0051 
0052 % Get spontaneous reactions associated with the pathways
0053 hits=[];
0054 for i=1:numel(pwys)
0055     [a, b]=ismember(pwys{i},pathways);
0056     if ~isempty(find(sprxnPwyMat(:,b)))
0057         hits=[hits;find(sprxnPwyMat(:,b))];
0058     end
0059 end
0060 spontaneousRxnList=isSpontaneous(unique(hits));
0061 
0062 % Check if the reactants/products are included in metList
0063 for i=1:numel(spontaneousRxnList)
0064     [a, b]=ismember(spontaneousRxnList{i},metaCycRxns.rxns);
0065     
0066     %Obtain the reactants and products
0067     reactants=all(ismember(metaCycRxns.mets(find(metaCycRxns.S(:,b)==-1)),metList));
0068     products=all(ismember(metaCycRxns.mets(find(metaCycRxns.S(:,b)==1)),metList));
0069     
0070     %Either reactants or products should be present for reversible reaction
0071     if metaCycRxns.rev(b)==1
0072         if ~(reactants || products)
0073             spontaneousRxnList{i}=[];
0074         end
0075         %Reactants should be present for irreversible reactions
0076     else
0077         if ~reactants
0078             spontaneousRxnList{i}=[];
0079         end
0080     end
0081 end
0082 spontaneousRxnList=spontaneousRxnList(~cellfun(@isempty, spontaneousRxnList));
0083 
0084 % Remove the spontaneousRxns that have already been included
0085 spontaneousRxnList=setdiff(spontaneousRxnList,rxnList);
0086 
0087 % Generate the cell array of relevant pathways
0088 pathwayID=cell(numel(spontaneousRxnList),1);
0089 for i=1:numel(spontaneousRxnList)
0090     index=find(strcmp(spontaneousRxnList{i},isSpontaneous));
0091     % locate the relevant pathways
0092     relevantpwys=intersect(pathways(find(sprxnPwyMat(index,:))),pwys);
0093     pathwayID{i}=strjoin(relevantpwys,';');
0094 end
0095 
0096 end