function MetabolicEssentiality(PotenGenes,algorithmName,modelType)
% Identifying essential metabolic genes in generated models and compare
% them with results of Luo et al, (2007) and COLT-Cancer database

% Oveis Jamialahmadi
% Jan 2018

if modelType == 1
    load Recon1_genes
elseif modelType == 2
    load Recon2_genes
else
    error('Wrong modelType!')
end

%//////////////////////////////////////////////////////////////////////////
% Read raw data:
% 1- COLT-Cancer database -------------------------------------------------
% load all data from COLT-Cancer including p-value
load COLTCancerData
%--------------------------------------------------------------------------
% 2- Luo et al. results
load LuoData
%--------------------------------------------------------------------------
% 3- DisGeNet database
load DisGeNet_curated
clear diseaseName
load DisGeNet_diseaseMapping
neoplasia_idx1 = find(ismember(lower(SemanticType),'neoplastic process'));
neoplasia_idx2 = find(ismember(lower(DiseaseClass),'neoplasms'));
neoplasia_idx = union(neoplasia_idx1,neoplasia_idx2);
neoplasia_DisGeNet_dieaseID = diseaseid(neoplasia_idx);
neoplasia_sharedID = intersect(neoplasia_DisGeNet_dieaseID,diseaseId);

neoplasia_idx = [];
for i = 1:numel(neoplasia_sharedID)
    neoplasia_idx = [neoplasia_idx;find(ismember(diseaseId,neoplasia_sharedID{i}))];
end
geneId = geneId(neoplasia_idx);
geneId = unique(geneId);
%//////////////////////////////////////////////////////////////////////////

% Check with COLT-Cancer database
% fprintf('Checking with COLT-Cancer database...\n')
Allfound = ({});
for ct = 1:size(PotenGenes,1)
%     fprintf('Checking with COLT-Cancer: model:%d of %d\n',ct,size(PotenGenes,1))
    cntr = 1;  cntr1 = 1; ctrNci = 0;ctrNci1 = 0;
    found_cells = ({}); found_cancer = ({}); found_cells1 = ({}); found_cancer1 = ({});
    genelist = PotenGenes{ct,1}; % for each model
    genelist = unique(genelist);
    genelist1 = Recon_genes; % All unique genes in Recon1
    find_ncbi_id = ncbi_id(ismember(ncbi_id,genelist));
    find_allpval = coltpval(ismember(ncbi_id,genelist),:);
    find_ncbi_id1 = ncbi_id(ismember(ncbi_id,genelist1));
    find_allpval1 = coltpval(ismember(ncbi_id,genelist1),:);
    for listCt = 1:size(find_allpval,1)
        find_each_pval = find(find_allpval(listCt,:)<0.05); %threshold to be essential
        if ~isempty(find_each_pval)
            for ct2 = 1:numel(find_each_pval)
                found_cells{cntr}{ct2} = colt_cells{find_each_pval(ct2)};
                found_cancer{cntr}(ct2) = colt_linename(ismember(colt_line,found_cells{cntr}{ct2}));
            end
            % Check if current model is among found cell lines in COLT
            if any(strcmp(found_cells{cntr},ncinames{ct}))
                ctrNci = ctrNci + 1;
            end
            cntr = cntr + 1;
        end
    end
    % For generic metabolic model
    for innerCt = 1:size(find_allpval1,1)
        find_each_pval1 = find(find_allpval1(innerCt,:)<0.05); %threshold to be essential
        if ~isempty(find_each_pval1)
            for ct21 = 1:numel(find_each_pval1)
                found_cells1{cntr1}{ct21} = colt_cells{find_each_pval1(ct21)};
                found_cancer1{cntr1}(ct21) = colt_linename(ismember(colt_line,found_cells1{cntr1}{ct21}));
            end
            % Check if current model is among found cell lines in COLT
            if any(strcmp(found_cells1{cntr1},ncinames{ct}))
                ctrNci1 = ctrNci1 + 1;
            end
            cntr1 = cntr1 + 1;
        end
    end
    Allfound{ct}.cells = found_cells;
    Allfound{ct}.cancer = found_cancer;
    Allfound{ct}.nci = ctrNci;
    Allfound{ct}.frac1 = numel(found_cells)/numel(genelist);
    Allfound{ct}.ess1 = numel(found_cells);
    Allfound{ct}.gene1 = numel(genelist);
    Allfound{ct}.ess2 = numel(found_cells1);
    Allfound{ct}.gene2 = numel(genelist1);
    Allfound{ct}.frac2 = numel(found_cells1)/numel(genelist1); % Generic metabolic model
    Allfound{ct}.enrich = Allfound{ct}.frac1/ Allfound{ct}.frac2;
    % Compute hypergeometric test
    Allfound{ct}.p= hygecdf(numel(found_cells)-1,numel(genelist1),numel(found_cells1),numel(genelist),'upper');
end
COLTraw = Allfound;
clear Allfound
%//////////////////////////////////////////////////////////////////////
% Check with DisGeNet database
% fprintf('Checking with DisGeNet database...\n')
Allfound = ({}); 
for ct = 1:size(PotenGenes,1)
%     fprintf('Checking with DisGeNet: model:%d of %d\n',ct,size(PotenGenes,1))
    genelist = PotenGenes{ct,1}; % for each model
    genelist = unique(genelist);
    genelist1 = Recon_genes; % All unique genes in Recon1

    found_cells = geneId(ismember(geneId,genelist));
    found_cells1 = geneId(ismember(geneId,genelist1));
    Allfound{ct}.frac1 = numel(found_cells)/numel(genelist);
    Allfound{ct}.ess1 = numel(found_cells);
    Allfound{ct}.gene1 = numel(genelist);
    Allfound{ct}.ess2 = numel(found_cells1);
    Allfound{ct}.gene2 = numel(genelist1);
    Allfound{ct}.frac2 = numel(found_cells1)/numel(genelist1); % Generic metabolic model
    Allfound{ct}.enrich = Allfound{ct}.frac1/ Allfound{ct}.frac2;
    % Compute hypergeometric test
    Allfound{ct}.p= hygecdf(numel(found_cells)-1,numel(genelist1),numel(found_cells1),numel(genelist),'upper');

end
DisGeNetraw = Allfound;
clear Allfound
%//////////////////////////////////////////////////////////////////////
% Check with Luo et al. results
% fprintf('Checking with Luo et al. results...\n')
ks_p = (0);

for cnt = 1:size(PotenGenes,1)
%     fprintf('Performing KS-test for model:%d of %d\n',cnt,size(PotenGenes,1))
    found_ranks = ({});
    for ct = 1:2
        genelist = PotenGenes{cnt,ct}; % for each model
        genelist = unique(genelist);
        genesym = gene_symbol(ismember(entrez_id,genelist));
        found_rank = riger_rank(ismember(riger_gene,genesym));
        found_rank(~found_rank) = [];
        found_ranks{ct} = found_rank;
    end
    if isempty(found_ranks{1}) || isempty(found_ranks{2})
        ks_p(cnt) = 9999;
    else
        [~,ks_p(cnt)] = kstest2(found_ranks{1},found_ranks{2},'Tail','larger');
    end
end
LuoRaw = ks_p;
clear ks_p


% Analyze and display results
MetEssAnalysis(COLTraw,LuoRaw,DisGeNetraw,algorithmName)
fprintf('============================================\n')
% ====================subfunctions=========================================
function MetEssAnalysis(COLTraw,LuoRaw,DisGeNetraw,algorithmName)

COLTres = struct('frac',0,'p',1,'pstd',1,'en',0,'enstd',0);
DisGeNetres = struct('frac',0,'p',1,'pstd',1,'en',0,'enstd',0);
Luores = struct('frac',0,'p',1,'pstd',1);

c = 1; this_p = (1); this_frac = (0); 
for count = 1:numel(COLTraw)
    Allfound = COLTraw{count};
    if Allfound.p < 0.05
        this_p(c) = Allfound.p;
        this_frac(c) = Allfound.enrich;
        c=c+1;
    end
    clear Allfound
end
COLTres.frac = (c-1)/numel(COLTraw);
COLTres.p = mean(this_p);
COLTres.pstd = std(this_p);
COLTres.en = mean(this_frac);
COLTres.enstd =std(this_frac);
    
c = 1; this_p = (1); this_frac = (0); 
for count = 1:numel(DisGeNetraw)
    Allfound = DisGeNetraw{count};
    if Allfound.p < 0.05
        this_p(c) = Allfound.p;
        this_frac(c) = Allfound.enrich;
        c=c+1;
    end
    clear Allfound
end
DisGeNetres.frac = (c-1)/numel(DisGeNetraw);
DisGeNetres.p = mean(this_p);
DisGeNetres.pstd = std(this_p);
DisGeNetres.en = mean(this_frac);
DisGeNetres.enstd =std(this_frac);

this_p = (1); c = 1;
for count = 1:numel(LuoRaw)
    if LuoRaw(count) < 0.05
        this_p(c) = LuoRaw(count);
        c=c+1;
    end
end
Luores.frac = (c-1)/numel(LuoRaw);
Luores.p = mean(this_p);
Luores.pstd = std(this_p);
clear ks_p

COLTres = squeeze(cell2mat(struct2cell(COLTres)))';
DisGeNetres = squeeze(cell2mat(struct2cell(DisGeNetres)))';
Luores = squeeze(cell2mat(struct2cell(Luores)))';

fprintf('COLT-Caner:....................................\n')
disp(table(COLTres(1),COLTres(2),COLTres(4),COLTres(5)...
     ,'VariableNames',{'f','pval','En','Enstd'},'RowNames',algorithmName'))

 fprintf('DisGeNet:....................................\n')
disp(table(DisGeNetres(1),DisGeNetres(2),DisGeNetres(4),DisGeNetres(5)...
 ,'VariableNames',{'f','pval','En','Enstd'},'RowNames',algorithmName'))

 fprintf('Luo et al.:....................................\n')
 disp(table(Luores(1),Luores(2),Luores(3)...
     ,'VariableNames',{'f','pval','En_p'},'RowNames',algorithmName'))
