Description of Bigg2Kegg

Home > BiKEGG > Bigg2Kegg.m

# Bigg2Kegg

## PURPOSE

BiKEGG

## SYNOPSIS

function Bigg2Kegg(Idfier)

## DESCRIPTION

     BiKEGG
     generates corresponding reactions in KEGG and BiGG/HMR
     databases and saves the output as a *.mat file in BiGG2KEGG folder.
     Input:
         - Idfier : a string of either 'bigg' or 'hmr'
     Full documentation can be found in the user manual

     O. Jamialahmadi
     TMU, Chem. Eng. Dept., Biotech. Group
     Oct. 2015

## CROSS-REFERENCE INFORMATION

This function calls:

- Bigg2KeggRestricted Bigg2KeggRestricted
- Bigg2KeggTable BIGG2KEGGTABLE MATLAB code for Bigg2KeggTable.fig
- MultiRxns MultiRxns

This function is called by:

## SUBFUNCTIONS

- function Cbs = MetCombs(Metkegg,Prd)
- function [Bg,Ar,Bg2,Ar2] = RxnPrune(BiGGID, AllRxn, RxnNames)
- function Keggrm = KeggApiMatch(TempMet)
- function Metkegg = ModelMoidfy (Metkegg,MetAbr)
- function [BiGGID,AllRxn] = AddLumpRxns(Steprxns,Insteprxns,BiGGID1,AllRxn1)
- function [TransKeggT,TransBiggT] = CheckTransRxns(TransKegg,TransBiGG)
- function[boutref,koutref,bout,kout] = biggref(bin,kin)

## SOURCE CODE

    0001 function Bigg2Kegg(Idfier)
    0002 % BiKEGG
    0003 % generates corresponding reactions in KEGG and BiGG/HMR
    0004 % databases and saves the output as a *.mat file in BiGG2KEGG folder.
    0005 % Input:
    0006 %     - Idfier : a string of either 'bigg' or 'hmr'
    0007 % Full documentation can be found in the user manual
    0008 %
    0009 % O. Jamialahmadi
    0010 % TMU, Chem. Eng. Dept., Biotech. Group
    0011 % Oct. 2015
    0012
    0013 Pth1 = which ('Bigg2Kegg.m');
    0014 tind = find(Pth1=='\',1,'last');
    0015 Pth = Pth1(1:tind-1);
    0016 if strcmp(Idfier,'hmr')
    0017     [pname,fname]=uigetfile('*.xlsx','Select the HMRdatabase xlsx file');
    0018     [~,dat2]=xlsread([fname,pname],'RXNS');
    0019     ind=dat2(:,1);
    0020     ind1=ismember(ind,'#');
    0021     ind2=find(ind1,3,'first');
    0022     hdrs=dat2(1,:);
    0023     % Reading the BIGG IDs:
    0024     BIGG=dat2(2:ind2(3)-1,ismember(hdrs,'BIGG DATABASE ID'));
    0025     % Reading the KEGG rxn IDs:
    0026     KEGG=dat2(2:ind2(3)-1,ismember(hdrs,'KEGG ID'));
    0027     % Reading the rxn equations:
    0028     Eqn=dat2(2:ind2(3)-1,ismember(hdrs,'EQUATION'));
    0029     % Locating empty cells
    0030     r1= ~ismember(KEGG,'');
    0031     r2= ~ismember(BIGG,'');
    0032     r3 = find(r1 & r2);
    0033     KEGG=KEGG(r3);BIGG=BIGG(r3);Eqn=Eqn(r3);
    0034     B2Kegg.Kegg=KEGG;B2Kegg.Bigg=BIGG;B2Kegg.Eqn=Eqn;
    0035     fname = fullfile(Pth,'BiGG2KEGG','BiGG2KEGG_HMRbased-RECON1.mat');
    0036     save (fname,'B2Kegg')
    0037     return
    0038 end
    0039
    0040 if ~strcmp(Idfier,'bigg')
    0041     msgbox('Unknown identifier: bigg/hmr','Error','error');
    0042     return
    0043 end
    0044 % Read text file extracted from http://rest.kegg.jp/link/rn/cpd, which
    0045 % contains corresponding rxns for each metabolite.
    0046 Fileid = fopen(fullfile(Pth,'cpd2rxn.txt'),'r');
    0047 cpd2rxn = textscan(Fileid,'%s %[^\n]');
    0048 rxnData = cell2mat(cpd2rxn{2}); % KEGG rxn IDs
    0049 cpds = strrep(cpd2rxn{1},'cpd:',''); % KEGG Cpds
    0050 fclose(Fileid);
    0051 % Read the model file======================================================
    0052 [pname,fname1]=uigetfile({'*.*'},'Select the input model file (mat|xml)');
    0053 Testname=pname(strfind(pname,'.')+1:end);
    0054 if strcmp(Testname,'mat')
    0055     D = importdata([fname1,pname]);
    0056     keggloc1 = fieldnames(D);
    0057     keggid = strcmpi(keggloc1,'metkeggid');
    0058     if ~any(keggid)
    0059         msgbox('There is no metKeggID field in the model!','Error','error');
    0060         return
    0061     end
    0062     Metkegg = D.(keggloc1{keggid});
    0063 %     NOTE: following commented lines can be activated to retrieve the
    0064 %     Metkegg online form BiGG database
    0065 %     Mets = D.mets;
    0066 %     for counter = 1:length(Mets)
    0067 %         BiggMet = urlread(['http://bigg.ucsd.edu/api/v2/models/',...
    0068 %             D.description,'/metabolites/',Mets{counter}]);
    0069 %         KeggLoci = strfind(BiggMet,'KEGG');
    0070 %         if ~isempty(KeggLoci)
    0071 %             Metkegg{counter} = BiggMet(KeggLoci+62:KeggLoci+67);
    0072 %         else
    0073 %             Metkegg{counter} = '';
    0074 %         end
    0075 %     end
    0076 elseif strcmp(Testname,'xml')
    0077     Idf = fopen([fname1,pname]);
    0078     D = textscan(Idf,'%s');
    0079     CbModel = readCbModel([fname1,pname]);
    0080 else
    0081     msgbox('Unknown filetype!','Error','error');
    0082     return
    0083 end
    0084 % =========================================================================
    0085 if strcmp(Testname,'xml')
    0086     % Reading xml
    0087     % Find listOfSpecies which contains metabolite annotations
    0088     SpeciesLoci1 = strfind(D{1},'listOfSpecies');
    0089     SpeciesLoci = find(~cellfun('isempty', SpeciesLoci1));
    0090     % Discard other data in the model:
    0091     SpeciesData = D{1}(SpeciesLoci(1):SpeciesLoci(2));
    0092     clear D
    0093     % Find each metabolite-------------------------------------------------
    0094     SpeciesStart1 = strfind(SpeciesData,'<species');
    0095     SpeciesStart = find(~cellfun('isempty', SpeciesStart1));
    0096     SpeciesEnd1 = strfind(SpeciesData,'</species>');
    0097     SpeciesEnd = find(~cellfun('isempty', SpeciesEnd1));
    0098     if numel(SpeciesStart) ~= numel(SpeciesEnd)
    0099         msgbox('SBML format error:Species have wrong format!','Error',...
    0100             'error');
    0101         return
    0102     end
    0103     % ---------------------------------------------------------------------
    0104     % Search for KEGG compounds for each BiGG metabolite. Note that each BiGG
    0105     % metabolite may have more than one equivalent KEGG compund, and although
    0106     % HMR and RECON2 (Thiele lab.) do not provide more than one KEGG compound
    0107     % for each metabolite, here, we followed the format of BiGG models (i.e.
    0108     % several equivalent KEGG compounds for a specific metabolite).
    0109     Metkegg = {''}; % KEGG compounds will fill this cell!
    0110     for count = 1:numel(SpeciesStart)
    0111         SpeciesTempStr = SpeciesData(SpeciesStart(count):SpeciesEnd(count));
    0112         SpeciesTempStr = [SpeciesTempStr{:}];
    0113         KeggCpdLoci = regexp(SpeciesTempStr,...
    0114             '(?<=http://identifiers.org/kegg.compound/)[CG]\d{5}','match');
    0115         if isempty(KeggCpdLoci)
    0116             Metkegg{count} = '';
    0117             continue
    0118         end
    0119         for count1 = 1:numel(KeggCpdLoci)
    0120             if strcmp(KeggCpdLoci{count1},'C00001') % H2O
    0121                 Metkegg{count}= {'C00001'};
    0122                 break
    0123             elseif strcmp(KeggCpdLoci{count1},'C00027') % H2O2
    0124                 Metkegg{count}= {'C00027'};
    0125                 break
    0126             else
    0127                 Metkegg{count}{count1}= KeggCpdLoci{count1};
    0128             end
    0129         end
    0130     end
    0131 end
    0132 % Model modification: some metabolites have wrong/absent KEGG annotations.
    0133 % Therefore, on the basis of modified KEGG compound IDs for RECON1.xml,
    0134 % these metabolites will be added to Metkegg separately.
    0135 D = CbModel;
    0136 ModelName = D.description;
    0137 tind1 = find(ModelName=='\',1,'last');
    0138 ModelName1 = ModelName(tind1+1:end-4);
    0139 if ~strcmp(ModelName1,'RECON1')
    0140     MetAbr = CbModel.mets;
    0141     Metkegg = ModelMoidfy (Metkegg,MetAbr);
    0142 end
    0143
    0144 % Search for Metkegg components with numel > 1, compare their names to
    0145 % BiGG metabolites, and remove other components with different names if a
    0146 % unique result is obtained. KEGGcpd.txt is used which contains all KEGG
    0147 % compounds.
    0148 Fileid1 = fopen(fullfile(Pth,'KEGGcpd.txt'),'r');
    0149 Keggcpd = textscan(Fileid1,'%s %[^\n]');
    0150 CpdData = strrep(Keggcpd{1},'cpd:','');
    0151 CpdNames = regexp(Keggcpd{2},';','split');
    0152 BiGGmets = D.metNames;
    0153 for count = 1:numel(Metkegg)
    0154     if all(ismember({'C01326','C00177'},Metkegg{count})) % For cyan
    0155         continue
    0156     end
    0157 % Second condition: L-carnitine definition discrepancy between KEGG and BiGG
    0158     if numel(Metkegg{count})>1 && ~all(ismember({'C00318','C00487'},...
    0159             Metkegg{count}))
    0160         MatchLoci = find(ismember(CpdData,Metkegg{count}));
    0161         for incount = 1:numel(MatchLoci)
    0162             CpdNames{MatchLoci(incount)}=strtrim(CpdNames{MatchLoci(incount)});
    0163             FoundL = find(ismember(BiGGmets,CpdNames{MatchLoci(incount)}));
    0164             if FoundL
    0165                 % Check if found cpd participates in any rxn
    0166                 if any(strcmp(CpdData(MatchLoci(incount)),cpds))
    0167                     Metkegg{count} = CpdData(MatchLoci(incount));
    0168                     break
    0169                 end
    0170             end
    0171         end
    0172     end
    0173 end
    0174
    0175 % Extract KEGG metabolites from BIGG API ==================================
    0176 Fileid3 = fopen(fullfile(Pth,'bigg_models_metabolites.txt'),'r');
    0177 biggmets = textscan(Fileid3,'%s %s %[^\n]');
    0178 fclose(Fileid3);
    0179 biggmetid = biggmets{1}; biggmetdes = biggmets{3}; clear biggmets
    0180 Metkegg1 = {''};
    0181 modelm = D.mets;
    0182 modelm=strrep(modelm,'[','_');modelm=strrep(modelm,']','');
    0183 for ct1 = 1:numel(modelm)
    0184     tempdes = biggmetdes(strcmp(biggmetid,modelm{ct1}));
    0185     temploci = regexp(tempdes,...
    0186         'http://identifiers.org/kegg.compound/C\S{5}','match');
    0187     temploci = temploci{1};
    0188     if isempty(temploci)
    0189         Metkegg1{ct1} = '';
    0190     else
    0191         for ct2 = 1:numel(temploci)
    0192             regtemp = regexp(temploci{ct2},'C\S{5}','match');
    0193             Metkegg1{ct1}{ct2} = regtemp{1};
    0194         end
    0195     end
    0196 end
    0197 % Consensus Metkegg: Metkegg1 from BiGG API and Metkegg from COBRA models
    0198 % are mixed to generate the most complete set of KEGG metabolites for a
    0199 % SBML model.
    0200 for cp1 = 1:numel(Metkegg1)
    0201     if isempty(Metkegg1{cp1}) || (numel(Metkegg1{cp1}) ~= numel(Metkegg{cp1}))
    0202         for cp2 = 1:numel(Metkegg{cp1})
    0203             Metkegg1{cp1}{cp2} = Metkegg{cp1}{cp2};
    0204         end
    0205     end
    0206 end
    0207 clear Metkegg
    0208 Metkegg = Metkegg1;
    0209 % ==========================================================================
    0210 % Check the existence of metabolites participating in Bigg rxns and map the
    0211 % Bigg rxns to KEGG rxns if possible. To do so, metabolites in Bigg rxns are
    0212 % compared to KEGG metabolites, and KEGG metabolite IDs are located in KEGG
    0213 % databse to identify the KEGG rxn IDs they involve in KEGG databse. If all
    0214 % KEGG metabolites share a unique KEGG rxn ID, then that rxn ID corresponds
    0215 % to the Bigg rxn from which the metabolites were extracted (restricted
    0216 % condition).
    0217 % To reduce the CPU time, rxns are first compared to rxns in PrunedRxns.txt
    0218 % and those in UniModelKEGG to exclude them from search process.
    0219 load('UniModelKEGG.mat')
    0220 UniB = B2Kegg.B; UniK = B2Kegg.K;
    0221 Umulti = find(ismember(UniB,'MULTIR'));
    0222 UniB = UniB(1:Umulti-1); UniK = UniK(1:Umulti-1);
    0223 clear B2Kegg
    0224 Fileid12 = fopen('PrunedRxns.txt','r');
    0225 Prunedata12 = textscan(Fileid12,'%s %s %d %[^\n]');
    0226 Testpruned = Prunedata12{1};
    0227 fclose(Fileid12);
    0228 RxnNames = D.rxns; % BiGG rxn names
    0229 Rxns = D.S; % It is assumed that this field is universal.
    0230 Keggrm = (0);
    0231 BiGGID = {0};
    0232 AllRxn = {0};
    0233 Uctr = 1; ChkTrans = 1; TransKegg = ({}); TransBiGG = ({});
    0234 for rn = 1:size(Rxns,2) % for all rxns
    0235     fprintf('Rxn %d of %d\n',rn,size(Rxns,2))
    0236     % ---------------------------------------------------------------------
    0237     %if current reaction exists in pruned reactions, or UniModelKEGG there
    0238     %is no need to further search for potential equivalent reactions.
    0239     if any(ismember(Testpruned,RxnNames{rn}))
    0240         continue
    0241     end
    0242     if any(ismember(UniB,RxnNames{rn}))
    0243         UniBloci = find(ismember(UniB,RxnNames{rn}));
    0244         for u1 = 1:numel(UniBloci)
    0245             BiGGID{Uctr} = UniB{UniBloci(u1)};
    0246             AllRxn{Uctr} = UniK{UniBloci(u1)};
    0247             Uctr = Uctr + 1;
    0248         end
    0249         continue
    0250     end
    0251     %----------------------------------------------------------------------
    0252     Prd = find(Rxns(:,rn)); % Find substrates and products for each rxn
    0253     CheckTransport = [Metkegg{Prd}]; % For transport rxns
    0254     if (numel(CheckTransport)/numel(unique(CheckTransport)) == 2) && ...
    0255             numel(Prd) > 2
    0256         continue
    0257     end
    0258     % 1st condition: All reactants and products have metKeggIDs.
    0259     % 2nd condition: Exclude all one element rxns
    0260     if (~sum(strcmp (Metkegg(Prd),'')) && (length(Prd) > 1))
    0261         Cbs = MetCombs(Metkegg,Prd); % All combinations (permutations) of
    0262         % KEGG compounds, regarding the fact that there is more than one
    0263         % KEGG cpd for some metabolites in the model.
    0264         for Outctr = 1:size(Cbs,1)
    0265             % Contains all KEGG cpds in a certain rxn
    0266             TempMet = unique(Cbs(Outctr,:));
    0267             ChkTempMet = Cbs(Outctr,:);
    0268             % Transport rxns with 2 cpds
    0269             % 1st cond: finding a transport rxn; 3 is as a result of Cbs
    0270             % structure.
    0271             % 2nd cond: Transport reactions with 2 metabolites; the reason
    0272             % for 1 is similar to the 1st cond.
    0273             % 3rd cond: No redundancy
    0274             if (numel(ChkTempMet)/numel(TempMet) == 3) ...
    0275                     && (numel(TempMet) == 1) ...
    0276                    && ~any(ismember(TransBiGG,RxnNames{rn}))
    0277                 fprintf('Transport rxn: %s\n',RxnNames{rn})
    0278                 TransKegg{ChkTrans}{1} = TempMet{1};
    0279                 TransKegg{ChkTrans}{2} = TempMet{1};
    0280                 TransBiGG{ChkTrans} = RxnNames{rn};
    0281                 ChkTrans = ChkTrans + 1;
    0282                 continue
    0283             end
    0284             if (numel(TempMet) == 2) && any(strcmp(TempMet,'C00080'))
    0285                 continue
    0286             end
    0287
    0288             if length(TempMet) > 1 % Exclude all exchange rxns
    0289                 KeggrmMethod = 'off'; %Offline
    0290                 switch KeggrmMethod
    0291                     case 'off'
    0292                         Id1 = ismember(cpds,TempMet);
    0293                         Keggrm = str2num(rxnData(Id1,5:end));
    0294                     case 'on'
    0295                         Keggrm = KeggApiMatch (TempMet);
    0296                 end
    0297                 % Identify unique KEGG rxn ID
    0298                 Uq = unique(Keggrm);
    0299                 N = histc(Keggrm, Uq);
    0300                 if isempty(max(N))
    0301                     continue
    0302                 end
    0303                 if (max(N) == numel(TempMet)) || ...
    0304                         (max(N) == numel(TempMet)-1) % Found repeated rxns
    0305                     UqTemp = Uq(N==max(N));
    0306                     % Convert rxn codes to string
    0307                     for count = 1:numel(UqTemp)
    0308                         UqTemp1{count} = strcat('rn:R',repmat('0',...
    0309                             [1,4- floor(log10(UqTemp(count)+eps))])...
    0310                             ,num2str(UqTemp(count)));
    0311                     end
    0312                     if numel(UqTemp) > 1
    0313                         for count = 1:length(UqTemp1)
    0314                             UqTempSum(count) = sum(ismember(cpd2rxn{2},...
    0315                                 UqTemp1{count}));
    0316                         end
    0317                         % Restrict condition: If only the number of metabolites
    0318                         % in min(UqTempSum)reaction are identical to those
    0319                         % of TempMet
    0320                         if min(UqTempSum)==numel(TempMet) || ...
    0321                                min(UqTempSum)+1 == numel(TempMet) || ...
    0322                                min(UqTempSum)-1 == numel(TempMet)
    0323                             In1 = find(UqTempSum == min(UqTempSum(:)));
    0324                             for Inctr = 1:numel(In1)
    0325                                 AllRxn{Uctr} = UqTemp1{In1(Inctr)}(4:end);
    0326                                 BiGGID {Uctr} = D.rxns{rn};
    0327                                 Uctr = Uctr + 1;
    0328                             end
    0329                         end
    0330                     else
    0331                         UqTempSum = sum(ismember(cpd2rxn{2},UqTemp1{1}));
    0332                         if min(UqTempSum)==numel(TempMet) || ...
    0333                                 min(UqTempSum)+1 == numel(TempMet) || ...
    0334                                 min(UqTempSum)-1 == numel(TempMet)
    0335                             AllRxn{Uctr} = UqTemp1{1}(4:end);
    0336                             BiGGID {Uctr} = D.rxns{rn};
    0337                             Uctr = Uctr + 1;
    0338                         end
    0339
    0340                     end
    0341                     UqTempSum = 0;
    0342                     UqTemp1 = {0};
    0343
    0344                 end
    0345             end
    0346             Keggrm = 0;
    0347         end
    0348     end
    0349 end
    0350 % Identify transport reactions (diffusion)
    0351 [TransKeggT,TransBiggT] = CheckTransRxns(TransKegg,TransBiGG);
    0352
    0353 if ~isempty(TransBiggT{1})
    0354     BiGGID = [BiGGID,TransBiggT];
    0355     AllRxn = [AllRxn,TransKeggT];
    0356 end
    0357
    0358 % Identify BiGG rxns for which there is more than one KEGG rxns, and prune
    0359 % BiGGIDs and AllRxns
    0360 [BiGGID,AllRxn,MultiBID,MultiRxn] = RxnPrune(BiGGID, AllRxn, RxnNames);
    0361 AllRxnTemp = {[]}; BiGGIDTemp = {[]}; Tctr = 1;
    0362 if any ([MultiBID{:}])
    0363     for ct = 1:numel(MultiBID)
    0364         for ci = 1:numel(MultiRxn{ct})
    0365             BiGGIDTemp{Tctr} = MultiBID{ct};
    0366             AllRxnTemp{Tctr} = MultiRxn{ct}{ci};
    0367             Tctr = Tctr + 1;
    0368         end
    0369     end
    0370     BiGGID = [BiGGID,BiGGIDTemp];
    0371     AllRxn = [AllRxn,AllRxnTemp];
    0372 end
    0373 % Check with restricted condition
    0374 [B1,K1] = Bigg2KeggRestricted(D,Metkegg,rxnData,cpds,cpd2rxn);
    0375 Bdiff = BiGGID(~ismember(BiGGID,B1));
    0376 Kdiff = AllRxn(~ismember(BiGGID,B1));
    0377 if any ([MultiBID{:}])
    0378     Fileid = fopen(fullfile(Pth,'PrunedRxns.txt'),'r');
    0379     Prunedata = textscan(Fileid,'%s %s %d %[^\n]');
    0380     Prunedrxns = Prunedata{1};
    0381     fclose(Fileid);
    0382     MultiBIDtemp = MultiBID(~ismember(MultiBID,Prunedrxns));
    0383     MultiRxntemp = MultiRxn(~ismember(MultiBID,Prunedrxns));
    0384     Bmultidiff = MultiBIDtemp(~ismember(MultiBIDtemp,Bdiff));
    0385     Kmultidiff = MultiRxntemp(~ismember(MultiBIDtemp,Bdiff));
    0386     Tctr = 1; BmultidiffT = {[]}; KmultidiffT = {[]};
    0387     for ct = 1:numel(Bmultidiff)
    0388         for ci = 1:numel(Kmultidiff{ct})
    0389             BmultidiffT{Tctr} = Bmultidiff{ct};
    0390             KmultidiffT{Tctr} = Kmultidiff{ct}{ci};
    0391             Tctr = Tctr + 1;
    0392         end
    0393     end
    0394     if ~isempty(BmultidiffT{1})
    0395         Bdiff = [Bdiff,BmultidiffT];
    0396         Kdiff = [Kdiff,KmultidiffT];
    0397     end
    0398 end
    0399 % Check with bigg rxn ref.
    0400 [bref,kref,Bdiff,Kdiff] = biggref(Bdiff,Kdiff);
    0401 % -----------------------------------
    0402 if ~isempty(Bdiff{1})
    0403     Diffdata = cell(numel(Bdiff),3);
    0404     Diffdata(:,1) = Bdiff; Diffdata(:,2) = Kdiff;
    0405     for ct = 1:numel(Bdiff)
    0406         Diffdata{ct,3} = 0;
    0407     end
    0408     setappdata(0,'Diffdata',Diffdata);
    0409     waitfor(Bigg2KeggTable)
    0410     Finalv=getappdata(0,'Finalvalidity');
    0411     Fileid = fopen(fullfile(Pth,'PrunedRxns.txt'),'a');
    0412     checkb = Finalv(:,1);
    0413     checkv = Finalv(:,3);
    0414     checkv = [checkv{:}]; checkb1 = ({});
    0415     for ct = 1:size(Finalv,1)
    0416         if any(strcmp(Finalv{ct,1},checkb1))
    0417             continue
    0418         end
    0419         checkb1{ct} = Finalv{ct,1};
    0420         checkloci = find(ismember(checkb,Finalv{ct,1}));
    0421         checkfinal = checkv(checkloci);
    0422         if ~any(checkfinal)
    0423             for ct1 = 1:numel(checkfinal)
    0424                 fprintf(Fileid,'\n%s\t%s\t%d',Finalv{checkloci(ct1),1},...
    0425                     Finalv{checkloci(ct1),2},...
    0426                     0);
    0427             end
    0428         else
    0429             for ct1 = 1:numel(checkfinal)
    0430                 if checkv(checkloci(ct1))
    0431                     fprintf(Fileid,'\n%s\t%s\t%d',Finalv{checkloci(ct1),1},...
    0432                         Finalv{checkloci(ct1),2},1);
    0433                 end
    0434             end
    0435         end
    0436     end
    0437     fclose(Fileid);
    0438     AT=getappdata(0);
    0439     AT1 = fieldnames(AT);
    0440     for co = 1:length(AT1)
    0441         if strcmp(AT1{co},'Diffdata') || strcmp(AT1{co},'Finalvalidity')
    0442             rmappdata(0,AT1{co});
    0443         end
    0444     end
    0445 end
    0446 if ~isempty(bref{1})
    0447     rmvtemp = ismember(BiGGID,bref);
    0448     BiGGID(rmvtemp) = [];
    0449     AllRxn(rmvtemp) = [];
    0450     BiGGID = [BiGGID,bref];
    0451     AllRxn = [AllRxn,kref];
    0452 end
    0453 % Identify BiGG rxns for which there is more than one KEGG rxns, and prune
    0454 % BiGGIDs and AllRxns
    0455 [BiGGID,AllRxn,MultiBID,MultiRxn] = RxnPrune(BiGGID, AllRxn, RxnNames);
    0456 AllRxnTemp = {[]}; BiGGIDTemp = {[]}; Tctr = 1;
    0457 if any ([MultiBID{:}])
    0458     for ct = 1:numel(MultiBID)
    0459         for ci = 1:numel(MultiRxn{ct})
    0460             BiGGIDTemp{Tctr} = MultiBID{ct};
    0461             AllRxnTemp{Tctr} = MultiRxn{ct}{ci};
    0462             Tctr = Tctr + 1;
    0463         end
    0464     end
    0465     BiGGID = [BiGGID,BiGGIDTemp];
    0466     AllRxn = [AllRxn,AllRxnTemp];
    0467 end
    0468
    0469 % Extract multi-step reactions from KEGG. This step may need internet
    0470 % connection
    0471 [Steprxns,Insteprxns] = MultiRxns(AllRxn);
    0472 [BiGGID,AllRxn] = AddLumpRxns(Steprxns,Insteprxns,BiGGID,AllRxn);
    0473
    0474 % Save output data
    0475 B2Kegg.B=BiGGID; B2Kegg.K=AllRxn;
    0476 FileName = [ModelName1,'KEGG.mat'];
    0477 fname = fullfile(Pth,'BiGG2KEGG',FileName);
    0478 save (fname,'B2Kegg')
    0479
    0480 if any ([MultiBID{:}])
    0481     fprintf('There exists %.0f BiGG rxns with multiple KEGG IDs:\n',...
    0482         numel(MultiBID));
    0483     for ct = 1:numel(MultiBID)
    0484         fprintf('%s:\t',MultiBID{ct});
    0485         fprintf('%s  ',MultiRxn{ct}{:});
    0486         fprintf('\n')
    0487     end
    0488 else
    0489     disp('There exists no BiGG rxns with multiple KEGG IDs for this model')
    0490 end
    0491
    0492 %==========================================================================
    0493 %%===========================Subfunctions==================================
    0494 %==========================================================================
    0495 function Cbs = MetCombs(Metkegg,Prd)
    0496 % OutMet = A 2D cell with rows and columns corresponding to combinations
    0497 % and compounds, respectively.
    0498 % OutMet generates all combinations among Metkegg components. This is done
    0499 % to search through KEGG REACTIONs for all possibilities.
    0500 SizeCell = cell(1,numel(Prd));
    0501 for count = 1:numel(Prd)
    0502     TempCounter = numel(Metkegg{Prd(count)});
    0503     SizeCell{count} = 1:TempCounter;
    0504 end
    0505 MetkeggT = Metkegg(Prd);
    0506 MetkeggTemp = [MetkeggT{:}];
    0507 OutMet = combvec(SizeCell{:}).';
    0508 for count = 1:size(OutMet,2)-1 % Corresponds to linearized Metkegg(Prd)
    0509     OutMet(:,count+1)=OutMet(:,count+1)+max(OutMet(:,count));
    0510 end
    0511 CbsT = MetkeggTemp(OutMet);
    0512 % Check if H+ is peresent in Cbs. As KEGG reactions are not charge
    0513 % balanced, several BiGG reactions will not be found. To resolve this
    0514 % issue, new rows withouth H+ will be appended to Cbs.
    0515 if any(strcmp(CbsT(:),'C00080'))
    0516     Cbs1 = CbsT;
    0517     [N1,N2] = find(strcmp(Cbs1,'C00080'));
    0518     LenCbs = 1:numel(Cbs1);
    0519     for count = 1:numel(N1)
    0520         NonH2Cpd = find(~ismember(LenCbs,N2(count)));
    0521         Cbs1{N1(count),N2(count)} = CbsT{N1(count),NonH2Cpd(1)};
    0522     end
    0523     Cbs = [CbsT;Cbs1];
    0524 else
    0525     % On the contrary to above-mentioned case, the opposite case is also
    0526     % observed. In such a case, H+ will be added as an extra metabolite to
    0527     % each row in CbsT, and final Cbs will contain previous set of
    0528     % metabolites plus H+.
    0529     Cbs1 = CbsT;
    0530     Cbs2 = Cbs1;
    0531     for count = 1:size(CbsT,1)
    0532         Cbs1{count,size(CbsT,2)+1} = 'C00080'; % Add H+
    0533         % Keep dimensional consistency:
    0534         Cbs2{count,size(CbsT,2)+1} = Cbs2{count,size(CbsT,2)};
    0535     end
    0536     Cbs = [Cbs1;Cbs2];
    0537 end
    0538 %--------------------------------------------------------------------------
    0539 function [Bg,Ar,Bg2,Ar2] = RxnPrune(BiGGID, AllRxn, RxnNames)
    0540 % This subfunction uses PrunedRxns.txt, a set of manually found equivalence
    0541 % KEGG rxn IDs, to prune current rxn mapping. This set of reactions are
    0542 % generated based on univeral reactions in BiGG databse, and can
    0543 % be utilized for every desired model.
    0544 Fileid = fopen('PrunedRxns.txt','r');
    0545 PrAll = textscan(Fileid,'%s %s %d %[^\n]');
    0546 PrB = PrAll{3};
    0547 PrBY = find(PrB); % To be added
    0548 PrBN = ~PrB; % To be removed
    0549 TBrxns = PrAll{1}; TKrxns = PrAll{2};
    0550 BrxnY = TBrxns(PrBY); KrxnY = TKrxns(PrBY);
    0551 BrxnN = TBrxns(PrBN);
    0552 % Which BiGGID elements are in BrxnY?
    0553 % [Fx1,Fx2] = ismember(BrxnY,BiGGID);
    0554 % Fx2(Fx2==0) = [];
    0555 % AllRxn(Fx2) = KrxnY(Fx1);
    0556 BiGGInBrxnY = BiGGID(ismember(BiGGID,BrxnY));
    0557 KeggOfBrxnY = KrxnY(ismember(BrxnY,BiGGID));
    0558 BrxnYTemp = BrxnY(ismember(BrxnY,BiGGID));
    0559 SharedRxns = cell(1,numel(BiGGInBrxnY));
    0560 Tcount = 1; TempSharedRxns = ({}); TempBiGGID = ({});
    0561 for count = 1:numel(BrxnYTemp)
    0562     T= find(ismember(BiGGInBrxnY,BrxnYTemp{count}));
    0563     T1 = find(ismember(BrxnYTemp,BiGGInBrxnY(T(1))));
    0564
    0565     for count2 = 1:numel(T)
    0566         % Extract only the first element:: T1(1); others are stored in
    0567         % TempSharedRxns and will be added to AllRxn
    0568         SharedRxns(T(count2)) = KeggOfBrxnY(T1(1));
    0569         if numel(T1)>1 && ~any(ismember(TempBiGGID,BiGGInBrxnY(T(1))))
    0570             for count3 = 2:numel(T1)
    0571                 TempSharedRxns{Tcount} = KeggOfBrxnY{T1(count3)};
    0572                 TempBiGGID{Tcount} = BrxnYTemp{T1(count3)};
    0573                 Tcount = Tcount + 1;
    0574             end
    0575         end
    0576     end
    0577
    0578 end
    0579 AllRxn(ismember(BiGGID,BrxnY)) = SharedRxns;
    0580 if ~isempty(TempSharedRxns)
    0581     Addlen = numel(AllRxn)+1:numel(AllRxn)+numel(TempSharedRxns);
    0582     AllRxn(Addlen) = TempSharedRxns;
    0583     BiGGID(Addlen) = TempBiGGID;
    0584 end
    0585 % Which BrxnY elements are not in BiGGID and should be added separately?
    0586 NotInBiGG = BrxnY(~ismember(BrxnY,BiGGID));
    0587 RxnNot = KrxnY(~ismember(BrxnY,BiGGID));
    0588 % Are NotInBiGG elements exist in the current model (Note: BrxnY contains
    0589 % universal BiGG reactions) ?
    0590 YaLoci = find(ismember(NotInBiGG,RxnNames));
    0591 Yf1 = NotInBiGG(YaLoci);
    0592 AddedRxns = RxnNot(YaLoci);
    0593 ModLength = (1:numel(Yf1))+numel(BiGGID);
    0594 BiGGID(ModLength) = Yf1;
    0595 AllRxn(ModLength) = AddedRxns;
    0596 % Remove other found rxns in BiGGID which their absence is manually
    0597 % confirmed.
    0598 TobeRemovedRxns = ismember(BiGGID,BrxnN);
    0599 BiGGID(TobeRemovedRxns)=[]; AllRxn(TobeRemovedRxns)=[]; % Remove other rxns
    0600 % Find BiGG IDs with more than one rxn id
    0601 ct=1;cr=1;
    0602 Bg2={0};Ar2={0};Bg={0};Ar={0};
    0603 BiGGT = BiGGID; RxnT = AllRxn;
    0604 while numel(BiGGT)
    0605     T = find(ismember(BiGGT,BiGGT(1))); R = RxnT(T);
    0606     if numel(unique(R))~=1 % Reactions with more than one KEGG ID
    0607         Bg2{ct} = BiGGT{1};
    0608         Ar2{ct} = unique(R);
    0609         ct = ct+1;
    0610     else
    0611         Bg{cr}=BiGGT{1};
    0612         Ar{cr}=cell2mat(unique(R));
    0613         cr=cr+1;
    0614     end
    0615     BiGGT(T)=[]; %Remove
    0616     RxnT(T)=[];
    0617 end
    0618 %--------------------------------------------------------------------------
    0619 function Keggrm = KeggApiMatch(TempMet)
    0620 ctr=1;
    0621 Keggrm=(0);
    0622 for mt = 1:length(TempMet)
    0623     Metstr = urlread(['http://rest.kegg.jp/link/rn/',TempMet{mt}]);
    0624     Id1 =  regexp(Metstr,'rn:R'); % All rxns containing a cpd
    0625     if Id1 % KEGG rxn IDs are provided for this cpd
    0626         for di = 1:length(Id1) % Extracts rxns
    0627             Keggrm(ctr) = str2double(Metstr(Id1(di)+4:Id1(di)+9));
    0628             ctr = ctr+1;
    0629         end
    0630     else
    0631         break % Don't generate Metstr for other compounds
    0632     end
    0633 end
    0634 % -------------------------------------------------------------------------
    0635 function Metkegg = ModelMoidfy (Metkegg,MetAbr)
    0636 % Read Metkegg data for RECON1.xml
    0637 Mdat = load('RECON1Metkegg.mat');
    0638 RECON1k = Mdat.RECON1meta.K;
    0639 RECON1m = Mdat.RECON1meta.B;
    0640 [Fx1,Fx2] = ismember(RECON1m,MetAbr);
    0641 Fx2(Fx2==0) = [];
    0642 % Metkegg(Fx2) = RECON1k(Fx1);
    0643 Fx3 = find(Fx1);
    0644 for i1 = 1:numel(Fx2)
    0645     if numel(Metkegg{Fx2(i1)}) > numel(RECON1k{Fx3(i1)})
    0646        Tempmet = Metkegg{Fx2(i1)};
    0647        whrc = ismember(Tempmet,RECON1k{Fx3(i1)});
    0648        if isempty(whrc)
    0649         Metkegg{Fx2(i1)} = RECON1k(Fx3(i1));
    0650         end
    0651     else
    0652         Metkegg{Fx2(i1)} = RECON1k{Fx3(i1)};
    0653     end
    0654 end
    0655 % -------------------------------------------------------------------------
    0656 function [BiGGID,AllRxn] = AddLumpRxns(Steprxns,Insteprxns,BiGGID1,AllRxn1)
    0657 % ExcRxns contains rxns which need to be removed manually from Steprxns
    0658 ExcRxns = {'R03082','R01651','R01210','R08767','R01651','R03171'};
    0659 FindLocs = find(ismember(Steprxns,ExcRxns));
    0660 Steprxns(ismember(Steprxns,ExcRxns)) = [];
    0661 Insteprxns(FindLocs)=[];
    0662 PrevLen = numel(BiGGID1)+1;
    0663 BiGGID = BiGGID1; AllRxn = AllRxn1;
    0664 BiGGID{PrevLen} = 'MULTIR'; AllRxn{PrevLen} = 'MULTIR';
    0665 CurrLen = numel(BiGGID1)+2;
    0666 for count = 1:numel(Steprxns)
    0667     if ~isempty(Insteprxns{count})
    0668         if ~any(ismember(Insteprxns{count},AllRxn1)) %No common rxns
    0669             BiggTemp = BiGGID1(ismember(AllRxn1,Steprxns{count}));
    0670             for count1 = 1:numel(BiggTemp)
    0671                 for count2 = 1:numel(Insteprxns{count})
    0672                     BiGGID{CurrLen} = BiggTemp{count1};
    0673                     AllRxn{CurrLen} = Insteprxns{count}{count2};
    0674                     CurrLen = CurrLen+1;
    0675                 end
    0676             end
    0677         else
    0678             TempLoc = find(~ismember(Insteprxns{count},AllRxn1));
    0679             BiggTemp = BiGGID1(ismember(AllRxn1,Steprxns{count}));
    0680             for count1 = 1:numel(BiggTemp)
    0681                 for count2 = 1:numel(TempLoc)
    0682                     BiGGID{CurrLen} = BiggTemp{count1};
    0683                     AllRxn{CurrLen} = Insteprxns{count}{TempLoc(count2)};
    0684                     CurrLen = CurrLen+1;
    0685                 end
    0686             end
    0687         end
    0688     end
    0689 end
    0690 %--------------------------------------------------------------------------
    0691 function [TransKeggT,TransBiggT] = CheckTransRxns(TransKegg,TransBiGG)
    0692 % Read cpd2weight ===================
    0693 Fileid = fopen('cpd2weight.txt','r');
    0694 cpd2weight = textscan(Fileid,'%s %f');
    0695 cpd4weight = strrep(cpd2weight{1},'cpd:','');
    0696 Weight = cpd2weight{2};
    0697 fclose(Fileid);
    0698 % Read rxns with 2, 4 and 6 components ========
    0699 Fileid = fopen('W4TranRxns2.txt','r');
    0700 rxn2t = textscan(Fileid,'%s %f %f');
    0701 rxn2 = rxn2t{1};
    0702 W21T = rxn2t{2};
    0703 W22T = rxn2t{3};
    0704 W2 = zeros(numel(W22T),2);
    0705 for ct = 1:numel(W21T)
    0706     W2(ct,1) = W21T(ct);
    0707     W2(ct,2) = W22T(ct);
    0708 end
    0709 fclose(Fileid);
    0710 Fileid = fopen('W4TranRxns4.txt','r');
    0711 rxn4t = textscan(Fileid,'%s %f %f %f %f');
    0712 rxn4 = rxn4t{1};
    0713 W41T = rxn4t{2};
    0714 W42T = rxn4t{3};
    0715 W43T = rxn4t{4};
    0716 W44T = rxn4t{5};
    0717 W4 = zeros(numel(W44T),2);
    0718 for ct = 1:numel(W44T)
    0719     W4(ct,1) = W41T(ct);
    0720     W4(ct,2) = W42T(ct);
    0721     W4(ct,3) = W43T(ct);
    0722     W4(ct,4) = W44T(ct);
    0723 end
    0724 fclose(Fileid);
    0725 CurLen = 1;
    0726 TransKeggT = {[]}; TransBiggT = {[]};
    0727 for ct = 1:numel(TransKegg)
    0728     Text1 = ['Checking equivalent transport reaction for ',TransBiGG{ct},'\n'];
    0729     fprintf(Text1)
    0730     RxnComNo = numel(TransKegg{ct});
    0731     RxnWeight = Weight(ismember(cpd4weight,TransKegg{ct}));
    0732     if RxnComNo == 2
    0733         rxn24 = rxn2;
    0734         AllW = W2;
    0735     elseif RxnComNo == 4
    0736         rxn24 = rxn4;
    0737         AllW = W4;
    0738     end
    0739     if isempty(RxnWeight)
    0740         continue
    0741     end
    0742     Reslt = bsxfun(@eq,AllW,RxnWeight);
    0743     [mID,~]=find(Reslt);
    0744     if ~isempty(mID) && (RxnComNo.*numel(unique(mID)) == numel(mID))
    0745         CandRxns1 = rxn24(unique(mID));
    0746         for ct1 = 1:numel(CandRxns1)
    0747             TransBiggT{CurLen} = TransBiGG{ct};
    0748             TransKeggT{CurLen} = CandRxns1{ct1};
    0749             CurLen = CurLen + 1;
    0750         end
    0751     end
    0752 end
    0753 %--------------------------------------------------------------------------
    0754 function[boutref,koutref,bout,kout] = biggref(bin,kin)
    0755 % Read BiGG data
    0756 bout = {[]}; kout = {[]};  boutref = {[]}; koutref = {[]};
    0757 Pth1 = which ('Bigg2Kegg.m');
    0758 tind = find(Pth1=='\',1,'last');
    0759 Pth = Pth1(1:tind-1);
    0760 Fileid = fopen(fullfile(Pth,'bigg_models_reactions.txt'),'r');
    0761 biggrxns = textscan(Fileid,'%s %[^\n]');
    0762 fclose(Fileid);
    0763 biggid = biggrxns{1}; biggdes = biggrxns{2};
    0764 clear biggrxns
    0765 ct1 = 1; ct11 = 1;
    0766 while numel(bin)
    0767     bintemp = bin{1};
    0768     kintemp = kin(ismember(bin,bin{1}));
    0769     bin1loc = find(ismember(bin,bin{1}));
    0770     destemp = biggdes(ismember(biggid,bintemp));
    0771     temploci = regexp(destemp,...
    0772         'http://identifiers.org/kegg.reaction/R\S{5}','match');
    0773     if isempty(temploci) || isempty(temploci{1})
    0774         for ct3 = 1:numel(bin1loc)
    0775             bout{ct1} = bintemp;
    0776             kout{ct1} = kintemp{ct3};
    0777             ct1 = ct1 + 1;
    0778         end
    0779     else
    0780         temploci = temploci{1};
    0781         for ct2 = 1:numel(temploci)
    0782             regtemp = regexp(temploci{ct2},'R\S{5}','match');
    0783             boutref{ct11} = bintemp;
    0784             koutref{ct11} = regtemp{1};
    0785             ct11 = ct11 + 1;
    0786         end
    0787     end
    0788     bin(bin1loc) = [];
    0789     kin(bin1loc) = [];
    0790 end