function [Metkegg,Metbigg] = KEGG2BiGGMets(D,CbModel)
% KEGG2BiGGMets
% is a subfunction of MapAdjuster for displaying compound's
% data on created customized metabolic map by NetDraw. The function
% extracts all metabolite correspondences between KEGG and BiGG based on
% the current COBRA model and data in BiGG API generated by
% AllKEGG2BiGGmet. The structure is very similar to that of used in
% Bigg2Kegg with minor adjustments.
% 
% Inputs:
% D: XML data of current COBRA model from BiGG
% CbModel: current COBRA model structure
% 
% Outputs:
% Metkegg: All KEGG compounds' IDs
% Metbigg: All BiGG compounds' IDs

% O. Jamialahmadi
% TMU, Chem. Eng. Dept., Biotech. Group 
% July 2016
% -------------------------------------------------------------------------

Pth1 = which ('Bigg2Kegg.m');
tind = find(Pth1=='\',1,'last');
Pth = Pth1(1:tind-1);
% Reading xml
% Find listOfSpecies which contains metabolite annotations 
SpeciesLoci1 = strfind(D{1},'listOfSpecies');
SpeciesLoci = find(~cellfun('isempty', SpeciesLoci1));
% Discard other data in the model:
SpeciesData = D{1}(SpeciesLoci(1):SpeciesLoci(2));
clear D
% Find each metabolite-------------------------------------------------
SpeciesStart1 = strfind(SpeciesData,'<species');
SpeciesStart = find(~cellfun('isempty', SpeciesStart1));
SpeciesEnd1 = strfind(SpeciesData,'</species>');
SpeciesEnd = find(~cellfun('isempty', SpeciesEnd1));
if numel(SpeciesStart) ~= numel(SpeciesEnd)
    msgbox('SBML format error:Species have wrong format!','Error',...
        'error');
    return
end
% ---------------------------------------------------------------------
% Search for KEGG compounds for each BiGG metabolite. Note that each BiGG
% metabolite may have more than one equivalent KEGG compund, and although
% HMR and RECON2 (Thiele lab.) do not provide more than one KEGG compound
% for each metabolite, here, we followed the format of BiGG models (i.e.
% several equivalent KEGG compounds for a specific metabolite).
Metkegg = {''}; % KEGG compounds will fill this cell!
for count = 1:numel(SpeciesStart)
    SpeciesTempStr = SpeciesData(SpeciesStart(count):SpeciesEnd(count));
    SpeciesTempStr = [SpeciesTempStr{:}];
    KeggCpdLoci = regexp(SpeciesTempStr,...
        '(?<=http://identifiers.org/kegg.compound/)[CG]\d{5}','match');
    if isempty(KeggCpdLoci)
        Metkegg{count} = '';
        continue
    end
    for count1 = 1:numel(KeggCpdLoci)
        if strcmp(KeggCpdLoci{count1},'C00001') % H2O
            Metkegg{count}= {'C00001'};
            break
        elseif strcmp(KeggCpdLoci{count1},'C00027') % H2O2
            Metkegg{count}= {'C00027'};
            break
        else
            Metkegg{count}{count1}= KeggCpdLoci{count1};
        end
    end
end
% Model modification: some metabolites have wrong/absent KEGG annotations.
% Therefore, on the basis of modified KEGG compound IDs for RECON1.xml,
% these metabolites will be added to Metkegg separately.
D = CbModel;
ModelName = D.description;
tind1 = find(ModelName=='\',1,'last');
ModelName1 = ModelName(tind1+1:end-4);
if ~strcmp(ModelName1,'RECON1')
    MetAbr = CbModel.mets;
    Metkegg = ModelMoidfy (Metkegg,MetAbr);
end

load(fullfile(Pth,'AllKEGG2BiGGmet.mat'))
Metkegg1 = AllKEGG2BiGGmet.Metkegg; Metkegg1 = [Metkegg1{:}];
Metbigg1 = AllKEGG2BiGGmet.Metbigg; Metbigg1 = [Metbigg1{:}];
MetAbr = strrep(MetAbr,'[','_'); MetAbr = strrep(MetAbr,']','');
Metkegg2 = ({}); ct = 1; MetAbr1 = ({});
for i1 = 1:numel(Metkegg)
    if isempty(Metkegg{i1})
        continue
    end
    for i2 = 1:numel(Metkegg{i1})
        Metkegg2{ct} = Metkegg{i1}{i2};
        MetAbr1{ct} = MetAbr{i1};
        ct = ct + 1;
    end
end
clear Metkegg 
Metbigg = [Metbigg1,MetAbr1]; Metkegg = [Metkegg1,Metkegg2];
Metbigg = regexp(Metbigg,'\w*(?=_\w{1}\>)','match'); % Remove compartments
Metbigg = [Metbigg{:}];
Mixdat = strcat(Metbigg,Metkegg);
[~, N1] = unique(Mixdat);
Metbigg = Metbigg(N1);
Metkegg = Metkegg(N1);

% Consensus Metkegg: Metkegg1 from BiGG API and Metkegg from COBRA models
% are mixed to generate the most complete set of KEGG metabolites for a
% SBML model.
% for cp1 = 1:numel(Metkegg1)
%     if isempty(Metkegg1{cp1}) || (numel(Metkegg1{cp1}) ~= numel(Metkegg{cp1}))
%         for cp2 = 1:numel(Metkegg{cp1})
%             Metkegg1{cp1}{cp2} = Metkegg{cp1}{cp2};
%         end
%     end
% end


function Metkegg = ModelMoidfy (Metkegg,MetAbr)
% Read Metkegg data for RECON1.xml
Mdat = load('RECON1Metkegg.mat');
RECON1k = Mdat.RECON1meta.K;
RECON1m = Mdat.RECON1meta.B;
[Fx1,Fx2] = ismember(RECON1m,MetAbr);
Fx2(Fx2==0) = [];
% Metkegg(Fx2) = RECON1k(Fx1);
Fx3 = find(Fx1);
for i1 = 1:numel(Fx2)
    if numel(Metkegg{Fx2(i1)}) > numel(RECON1k{Fx3(i1)})
       Tempmet = Metkegg{Fx2(i1)};
       whrc = ismember(Tempmet,RECON1k{Fx3(i1)});
       if isempty(whrc)
        Metkegg{Fx2(i1)} = RECON1k(Fx3(i1));
        end
    else
        Metkegg{Fx2(i1)} = RECON1k{Fx3(i1)};
    end
end

