MultiRxns

 MultiRxns
 gets a set of reactions and based on data provided in KEGG
 database for each reaction, searches for multi-step reactions and saves
 each multi-step reaction along with all of its single-step reactions
 (Child reactions).
 Inputs:
 K = A cell array of KEGG reactions identifiers.

0001 function [Po,Hout1] = MultiRxns(K)
0002 % MultiRxns
0003 % gets a set of reactions and based on data provided in KEGG
0004 % database for each reaction, searches for multi-step reactions and saves
0005 % each multi-step reaction along with all of its single-step reactions
0006 % (Child reactions).
0007 % Inputs:
0008 % K = A cell array of KEGG reactions identifiers.
0009 
0010 % Output:
0011 % Po = A cell array of multi-step reactions found in K
0012 % Hout1 = A cell array of decomposed single-step reactions corresponding to
0013 %         each reaction in Po
0014 %
0015 % O. Jamialahmadi
0016 % TMU, Chem. Eng. Dept., Biotech. Group
0017 % Nov. 2015
0018 % EDITED:
0019 % Jan. 2016 : Offline mode was included for reducing the cpu time
0020 
0021 NetStat = 0; % Offline mode
0022 Pth1 = which ('Bigg2Kegg.m');
0023 tind = find(Pth1=='\',1,'last');
0024 Pth = Pth1(1:tind-1);
0025 Pth3 = fullfile(Pth,'BiGG2KEGG\Multirxns.mat');
0026 Pth2 = fullfile(Pth,'OfflineRxnData');
0027 if exist(Pth3,'file')
0028     load (Pth3)
0029     if exist('Multirxns','var')
0030         Crxns = Multirxns.C;
0031         Prxns = Multirxns.P;
0032         PoT = Prxns(ismember(Prxns,K));
0033         Hout1T = Crxns(ismember(Prxns,K));
0034         K = K(~ismember(K,Prxns));
0035     end
0036 end
0037 if isempty(K)
0038     Po = PoT;
0039     Hout1 = Hout1T;
0040     return
0041 end
0042 % Check internet connection
0043 if NetStat
0044     TestLink='http://rest.kegg.jp';
0045     [~,Stat1]=urlread(TestLink);
0046     if ~Stat1
0047         msgbox({'Function: MultiRxns needs a stable internet connection';...
0048             'Seemingly you are offline!'});
0049         return
0050     end
0051 end
0052 
0053 ct=1;
0054 for count = 1:numel(K)
0055     cts=1; setappdata(0,'cts',cts)
0056     if exist('Po','var') && any(strcmp(K{count},Po))
0057         continue % No need to check again
0058     end
0059     if NetStat
0060         A = urlread(['http://rest.kegg.jp/get/',K{count}]);
0061     else
0062         load([Pth2,'\',K{count},'.mat']);
0063         A = RnDat;
0064         clear RnDat
0065     end
0066     [~,c2]=regexp(A,'COMMENT\s');
0067     if ~isempty(c2)
0068         [c3,~]=regexp(A,'RPAIR');
0069         if isempty(c3)
0070             [c3,~]=regexp(A,'ENZYME');
0071         end
0072         Temp = A(c2:c3-1);
0073         Tip1 = regexp(Temp,'\w*-step\>');
0074         if ~isempty(Tip1)
0075             Temp2 = Temp(Tip1:end);
0076             c4 = regexp(Temp2,'R\d{4}');
0077             if ~isempty(c4)
0078                 c5 = regexp(Temp2,'+');
0079                 c54 = sort([c4,c5]);
0080                 c5Loci = find(ismember(c54,c5));
0081                 c5LociP = c5Loci+1; c5LociM = c5Loci-1;
0082                 c5Loci1 = unique([c5LociP,c5LociM]);
0083                 c6 = regexp(A,'part of', 'once');
0084                 c7 = regexp(A,'second step of', 'once');
0085                 c8 = regexp(A,'first step of', 'once');
0086                 if ~isempty(c5) && isempty(c6) && isempty(c7) && isempty(c8)
0087                     Po{ct}=K{count};
0088                     for lop = 1:numel(c5Loci1)
0089                         Pol{ct}{lop} = Temp2(c54(c5Loci1(lop)):c54(c5Loci1(lop))+5);
0090                     end
0091                     Hout = SecLump(Pol{ct});
0092                     if ~isempty(Hout)
0093                         Hout1{ct} = SecCheck(Hout);
0094                         if ~isempty(Hout1{ct})
0095                             fprintf('Multi-step reaction was found:%d of %d:%s\n',...
0096                                 count,numel(K),K{count});
0097                         end
0098                         ct=ct+1;
0099                     end
0100                 end
0101             end
0102         end
0103         clear c2 c3 c4 c5 c6 c7 c8
0104     end
0105 end
0106 if exist('Po','var')
0107     Po = [Po,PoT];
0108 else
0109     Po = PoT;
0110 end
0111 if exist('Hout1','var')
0112     Hout1 = [Hout1,Hout1T];
0113 else
0114     Hout1 = Hout1T;
0115 end
0116 
0117 function Hout = SecLump(Hin)
0118 Pth1 = which ('Bigg2Kegg.m');
0119 tind = find(Pth1=='\',1,'last');
0120 Pth = Pth1(1:tind-1);
0121 Pth2 = fullfile(Pth,'OfflineRxnData');
0122 NetStat = 0;
0123 ct=1;
0124 for count = 1:length(Hin)
0125     if NetStat
0126         A = urlread(['http://rest.kegg.jp/get/',Hin{count}]);
0127     else
0128         load([Pth2,'\',Hin{count},'.mat']);
0129         A = RnDat;
0130         clear RnDat
0131     end
0132     [~,c2]=regexp(A,'COMMENT\s');
0133     if ~isempty(c2)
0134         [c3,~]=regexp(A,'RPAIR');
0135         if isempty(c3)
0136             [c3,~]=regexp(A,'ENZYME');
0137         end
0138         Temp = A(c2:c3-1);
0139         Tip1 = regexp(Temp,'\w*-dependent enzyme\>');
0140         if ~isempty(Tip1)
0141             Temp2 = Temp(Tip1:end);
0142             c4 = regexp(Temp2,'R\d{4}');
0143             for lop = 1:numel(c4)
0144                 Pol{ct}{lop} = Temp2(c4(lop):c4(lop)+5);
0145             end
0146             Hout = SecLump(Pol{ct});
0147             ct = ct+1;
0148         else
0149             
0150             cts=getappdata(0,'cts');
0151             Hout{cts} = Hin{count};
0152             cts=cts+1;
0153             setappdata(0,'cts',cts) 
0154         end
0155     else
0156         cts=getappdata(0,'cts');
0157         Hout{cts} = Hin{count};
0158         cts=cts+1;
0159         setappdata(0,'cts',cts)
0160     end
0161 end
0162 
0163 function Hout1 = SecCheck (Hin1)
0164 Fileid1 = fopen('rxn2map.txt','r');
0165 rxn2map = textscan(Fileid1,'%s %s');
0166 RawRxns = rxn2map{1};
0167 RawMaps = rxn2map{2};
0168 [Loci1,~] = regexp(RawMaps, 'path:rn');
0169 Loci2 = ~cellfun('isempty', Loci1);
0170 RawMaps(Loci2) = [];
0171 RawRxns(Loci2) = [];
0172 Maps = strrep(RawMaps,'path:map','');
0173 Rxns = strrep(RawRxns,'rn:','');
0174 MapsTemp={0};
0175 for count = 1:length(Hin1)
0176     MapsTemp{count} = Maps(ismember(Rxns,Hin1{count}))';
0177     if isempty(MapsTemp{count})
0178         Hout1 = [];
0179         break
0180     end
0181 end
0182 if ~exist('Hout1','var')
0183     MapsT = [MapsTemp{:}];
0184     MapsTU = unique(MapsT);
0185     if numel(MapsT)>numel(MapsTU)
0186         Hout1 = Hin1;
0187     end
0188 end