function MetabolicEssentiality_sp(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

PotenGenes_main = PotenGenes;
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

% 3- DisGeNet database
load DisGeNet_curated
clear diseaseName
load DisGeNet_diseaseMapping

sp_cat = {'breast_neoplasm_Id','CNS_neoplasm_Id','colon_neoplasm_Id','leukemia_neoplasm_Id',...
    'lung_neoplasm_Id','melanoma_neoplasm_Id','ovarian_neoplasm_Id','prostate_neoplasm_Id',...
    'renal_neoplasm_Id'};
id_cat = {'brst','brst','cns','cns','co_59','co_60','le_59','le_60','lc_59','lc_60',...
    'me_59','me_60','ov','ov','pr','pr','re_59','re_60'};

    brst = [13   15    27    37    50];
    ov = [3     6    19    20    34    38    42];
    re_60 = [1     2    24    26    30    40    54    57];
    re_59 = [1     2    24    26    30    40    53    56];
    co_60 = [4    18    21    25    28    33    59];
    co_59 = [4    18    21    25    28    33    58];
    lc_60 = [5     8     9    16    23    32    41    52    55];
    lc_59 = [5     8     9    16    23    32    41    51    54];
    cns = [17    29    43    44    46    49];
    me_60 = [7    10    12 14    22    31    36    39 45    60];
    me_59 = [7    10    12 14    22    31    36    39 45   59];
    pr = [11    47];
    le_60 = [ 35    48    51    53    56    58];
    le_59 = [ 35    48    52    55    57];


    brst_colt = 1:29;
    ov_colt = 30:44;

for sp = 1:numel(sp_cat)
    fprintf('Cancer: %s\n',strrep(sp_cat{sp},'_neoplasm_Id',''))
    load DisGeNet_curated
    clear diseaseName
    load DisGeNet_diseaseMapping
    load([sp_cat{sp},'.mat'])
    eval(['neoplasia_sharedID = intersect(',sp_cat{sp},',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);
    %//////////////////////////////////////////////////////////////////////////
    COLTraw = ({}); DisGeNetraw = ({});

    load COLTCancerData
    if strcmp(sp_cat{sp},'breast_neoplasm_Id')
        colt_cells = colt_cells(brst_colt);
        coltpval = coltpval(:,brst_colt);
    end
    if strcmp(sp_cat{sp},'ovarian_neoplasm_Id')
        colt_cells = colt_cells(ov_colt);
        coltpval = coltpval(:,ov_colt);
    end

   PotenGenes = PotenGenes_main;
    if size(PotenGenes,1) < 60
        eval(['PotenGenes = PotenGenes(',id_cat{2*sp-1},',:);'])
    else
        eval(['PotenGenes = PotenGenes(',id_cat{2*sp},',:);'])
    end

    Allfound = ({});
    for ct = 1: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
    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

    % Analyze and display results
    MetEssAnalysis(COLTraw,DisGeNetraw,algorithmName,id_cat(2*sp-1))
    fprintf('====================================\n')
end
fprintf('===============DONE=========================\n')
% ====================subfunctions=========================================
function MetEssAnalysis(COLTraw,DisGeNetraw,algorithmName,id_cat)

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

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);

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

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

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