% Analyze differences between three clusters on a sample basis using
% results generated with post-index NMPC data with 3 clusters with
% mine_cdi_model
% Michael A, Henson, 5/20/2020
% Modified to add errorbars to plots
% Michael A, Henson, 8/24/2020

% Find samples indices and patients belonging to the three clusters
[maxClVal,iClMax] = max(recurNum./nSampCl);
idxClMax = idxHigh(idxSamp(1:nSampCl(iClMax),iClMax));
patMax = unique(text16S(2,idxClMax+1),'stable');
[minClVal,iClMin] = min(recurNum./nSampCl);
idxClMin = idxHigh(idxSamp(1:nSampCl(iClMin),iClMin));
patMin = unique(text16S(2,idxClMin+1),'stable');
iClMed = setdiff([1:numCl],[iClMax iClMin]);
idxClMed = idxHigh(idxSamp(1:nSampCl(iClMed),iClMed));
patMed = unique(text16S(2,idxClMed+1),'stable');

% Find significantly different NMPCs between two clusters of interest for
% post-index samples
% iClMax, idxClMax = high recurrence, high Enterobactericeae cluster
% iClMed, idxClMed = medium recurrence, high Enterococcus cluster
% iClMin, idxClMin = medium recurrence, high Bacteroides cluster
iCl1 = iClMax;
iCl2 = iClMed;
idxCl1 = idxClMax;
idxCl2 = idxClMed;
for i = 1:size(dataMp,1)
    numSet1 = dataMp(i,idxCl1);
    numSet2 = dataMp(i,idxCl2);
    pMet(i) = ranksum(numSet1,numSet2);
end

% Correct p-values for false positive discovery due to multiple hypothesis
% testing
fdrMet = mafdr(pMet,'BHFDR',true);
idxMet = find(fdrMet < 0.05);
recSig = [];
for i = idxMet
    numSet1 = dataMp(i,idxCl1);
    numSet2 = dataMp(i,idxCl2);
    if fdrMet(i) < 0.05 && (mean(numSet1) > 50 || mean(numSet2) > 50) && (mean(numSet1)/mean(numSet2) < 0.5 || mean(numSet1)/mean(numSet2) > 2)
        recSig = [recSig; [i fdrMet(i)]];
    end
end
nSig = length(recSig);

% Plot significantly different NMPCs between two clusters of interest for
% post-index samples as Figure S5, S6 or S7 in paper
figure(21)
sigIdx = recSig(:,1);
tempLabels = erase(textMp(sigIdx+1,1),"EX_");
tempLabels = erase(tempLabels,"[fe]");
xLabels5 = strrep(tempLabels,'_','-');
dataSet1 = mean(dataMp(sigIdx,idxCl1)')';
dataSet2 = mean(dataMp(sigIdx,idxCl2)')';
nPlt = round(nSig/2);
subplot(2,1,1)
bar([dataSet1(1:nPlt) dataSet2(1:nPlt)]);
set(gca, 'XTick', 1:nPlt, 'XTickLabel', xLabels5(1:nPlt));
xtickangle(45)
ylabel('Maximal production (mmol/day)')
legend("High Bacteroides cluster, " + nSampCl(iCl1) + " samples",...
    "High Enterococcus cluster, " + nSampCl(iCl2) + " samples")
set(gca,'FontSize',12)
subplot(2,1,2)
bar([dataSet1(nPlt+1:end) dataSet2(nPlt+1:end)]);
set(gca, 'XTick', 1:length(xLabels5(nPlt+1:end)), 'XTickLabel', xLabels5(nPlt+1:end));
xtickangle(45)
ylabel('Maximal production (mmol/day)')
set(gca,'FontSize',12)

% Specify metabolites to be analyzed in more detail: BA = bile acids
baMets = {'cholate'; 'C02528'; '12dhchol'; '3dhcdchol'; '3dhchol'; '7dhcdchol'; '7ocholate'; 'HC02191';...
    'HC02194'; 'adchac'; 'alchac'; 'dchac'; 'icdchol'; 'isochol'; 'uchol'};
metNames = erase(textMp(:,1),"EX_");
metNames = erase(metNames,"[fe]");
nBA = length(baMets);
for i = 1:nBA
    idxBA(i) = find(strcmp(metNames,baMets{i}))-1;
end

% Calculate mean and standard error of the mean of each BA metabolite for
% each cluster
dataMax = mean(dataMp(idxBA,idxClMax)')';
semMax = std(dataMp(idxBA,idxClMax)')'/sqrt(nSampCl(iClMax));
dataMed = mean(dataMp(idxBA,idxClMed)')';
semMed = std(dataMp(idxBA,idxClMed)')'/sqrt(nSampCl(iClMed));
dataMin = mean(dataMp(idxBA,idxClMin)')';
semMin = std(dataMp(idxBA,idxClMin)')'/sqrt(nSampCl(iClMin));
meanData = [dataMax dataMed dataMin];
semData = [semMax semMed semMin];

% Plot mean and standard deviation of each BA metabolite for each cluster
% as Figure 3 in paper
figure(22)
subplot(2,1,1)
b = bar(meanData,'grouped');
hold on
ngroups = size(meanData,1);
nbars = size(meanData,2);
groupwidth = min(0.8, nbars/(nbars + 1.5));
for i = 1:nbars
    x = (1:ngroups) - groupwidth/2 + (2*i-1) * groupwidth / (2*nbars);
    errorbar(x, meanData(:,i), semData(:,i), 'k', 'linestyle', 'none');
end
hold off
ylim([0 400])
tempLabels = erase(textMp(idxBA+1,1),"EX_");
tempLabels = erase(tempLabels,"[fe]");
xLabels6 = strrep(tempLabels,'_','-');
set(gca, 'XTick', 1:nBA, 'XTickLabel', xLabels6);
xtickangle(45)
ylabel('Maximal production (mmol/day)')
legend("High Enterobactericeae cluster, " + nSampCl(iClMax) + " samples",...
    "High Bacteriodes cluster, " + nSampCl(iClMed) + " samples",...
    "High Enterococcus cluster, " + nSampCl(iClMin) + " samples")
set(gca,'FontSize',12)

% Specify metabolites to be analyzed in more detail: AAA = aromatic amino
% acids
aaaMets = {'2hyoxplac'; '34dhpha'; '4hoxpacd'; '4hphac'; 'Lkynr'; 'acac'; 'chor'; 'dhpppn'; 'ind3ppa';...
    'indole'; 'indprp'; 'pac'; 'pcresol'; 'phe_L'; 'phenol'; 'phppa'; 'phpyr'; 'pppn'; 'trp_L'; 'trypta';...
    'tym'; 'tyr_L'};
nAAA = length(aaaMets);
for i = 1:nAAA
    idxAAA(i) = find(strcmp(metNames,aaaMets{i}))-1;
end

% Calculate mean and standard deviation of each AAA metabolite for each
% cluster
dataMax1 = mean(dataMp(idxAAA,idxClMax)')';
semMax1 = std(dataMp(idxAAA,idxClMax)')'/sqrt(nSampCl(iClMax));
dataMed1 = mean(dataMp(idxAAA,idxClMed)')';
semMed1 = std(dataMp(idxAAA,idxClMed)')'/sqrt(nSampCl(iClMed));
dataMin1 = mean(dataMp(idxAAA,idxClMin)')';
semMin1 = std(dataMp(idxAAA,idxClMin)')'/sqrt(nSampCl(iClMin));
meanData1 = [dataMax1 dataMed1 dataMin1];
semData1 = [semMax1 semMed1 semMin1];

% Plot mean and standard deviation of each BA metabolite for each cluster
figure(22)
subplot(2,1,2)
b1 = bar(meanData1,'grouped');
hold on
ngroups1 = size(meanData1,1);
nbars1 = size(meanData1,2);
groupwidth1 = min(0.8, nbars1/(nbars1 + 1.5));
for i = 1:nbars1
    x1 = (1:ngroups1) - groupwidth1/2 + (2*i-1) * groupwidth1 / (2*nbars1);
    errorbar(x1, meanData1(:,i), semData1(:,i), 'k', 'linestyle', 'none');
end
hold off
ylim([0 400])
tempLabels = erase(textMp(idxAAA+1,1),"EX_");
tempLabels = erase(tempLabels,"[fe]");
xLabels6 = strrep(tempLabels,'_','-');
set(gca, 'XTick', 1:nAAA, 'XTickLabel', xLabels6);
xtickangle(45)
ylabel('Maximal production (mmol/day)')
set(gca,'FontSize',12)