% Script to mine NMPC predictions and metadata for for CDI FMT samples
% Michael A. Henson, 2/16/2020
% Modified to add errorbars to plots
% Michael A, Henson, 8/25/2020

clearvars

% Read unnormalized sample data and find samples with "high" unnormalized
% abundance
raw16Ssum = xlsread('fmt_abd_nmpc.xlsx',4,'C96:AP96');
idxHigh = find(ismember(raw16Ssum > 0.85, 1));

% Read normalized sample data
[data16S,text16S,all16S] = xlsread('fmt_abd_nmpc.xlsx',4,'A99:AP138');
abd16S = xlsread('fmt_abd_nmpc.xlsx',4,'C100:AP138');
X = abd16S(:,idxHigh);
numSamp = size(X,2);

% Read NMPC data
[dataMp,textMp,allMp] = xlsread('fmt_abd_nmpc.xlsx',5,'A1:AO413');
Mp = dataMp(:,idxHigh);

% Cluster NMPC data
numCl = 2;
idxSamp = zeros(size(Mp,2),numCl);
iCl = kmeans(Mp',numCl,'Replicates',1000);
for i = 1:numCl
    idx = find(ismember(iCl == i, 1));
    idxSamp(1:length(idx),i) = idx;
end
nSampCl = sum(idxSamp ~= 0);
nSampCl(numCl+1) = sum(nSampCl(1:numCl));

% Average abundance data over each cluster
abdCl = zeros(numCl,size(X,1));
for i = 1:numCl
    idx = nonzeros(idxSamp(:,i));
    if length(idx) == 1
        abdCl(i,:) = X(:,idx)';
    else
        abdCl(i,:) = mean(X(:,idx)');
    end
end
abdCl = abdCl';

% Find taxa for which average abundance in at least one cluster is above a
% threshold
threshAbd = 0.05;
idxAbd = [];
for i = 1:numCl
    idxAbd = union(idxAbd,find(ismember(abdCl(:,i) > threshAbd, 1)));
end

% Generate Figure 5 in paper

% Plot selected taxa abundances for each cluster
figure(55)
subplot(2,2,1)
bar(abdCl(idxAbd,:))
xLabels1 = erase(text16S(2:end,1),"pan");
set(gca,'XTick',1:length(idxAbd),'XTickLabel',xLabels1(idxAbd));
xtickangle(25)
ylabel('Average abundance')
legend("Cluster 1, " +  nSampCl(1) + " samples", "Cluster 2, " +  nSampCl(2) + " samples")

% Find type of each clustered sample
sampTypeLabel = text16S(1,3:end);
sampType = zeros(numSamp,numCl);
for i = 1:numCl
    idx = find(contains(sampTypeLabel(idxSamp(1:nSampCl(i),i)),'pre'));
    sampType(idx,i) = 1;
    idx = find(contains(sampTypeLabel(idxSamp(1:nSampCl(i),i)),'don'));
    sampType(idx,i) = 2;
        idx = find(contains(sampTypeLabel(idxSamp(1:nSampCl(i),i)),'post'));
    sampType(idx,i) = 3;
end

% Calculate number of pre, donor and post samples in each cluster
for i = 1:numCl
    for j = 1:3
        sampTypeNum(i,j) = sum(sampType(1:nSampCl(i),i) == j);
    end
end
for j = 1:3
    sampTypeNum(numCl+1,j) = sum(sampTypeNum(1:numCl,j));
end

% Plot number of pre, donor and post samples in each cluster
figure(55)
subplot(2,2,4)
bar([sampTypeNum])
xLabels2 = { 'Cluster 1'; 'Cluster 2'; 'Total'};
set(gca,'XTick',1:length(xLabels2),'XTickLabel',xLabels2);
ylabel('Number of samples')
legend('Pre-FMT', 'Donor', 'Post-FMT')

% Test clustered samples for significant differences between 3 sample types
for i = 1:numCl+1
    for j = i:numCl+1
        for k = 1:3
            occurData = [sampTypeNum(i,k) nSampCl(i)-sampTypeNum(i,k); sampTypeNum(j,k) nSampCl(j)-sampTypeNum(j,k)];
            [h,p,stats] = fishertest(occurData);
            if p < 0.05
                display(['Clusters ' num2str(i), ',' num2str(j) ' significant for sample type ' num2str(k) ' with ' num2str(p)])
            end
        end
    end
end

% Perform PCA on NMPC data
[coeff,score,latent,tsquared,explained,mu] = pca(Mp');
nSamp = size(Mp,2);

% Find pre-FMT, donor and post-FMT samples
idxPre = find(contains(sampTypeLabel,'pre'));
numSampPre = length(idxPre);
idxDon = find(contains(sampTypeLabel,'don'));
numSampDon = length(idxDon);
idxPos = find(contains(sampTypeLabel,'post'));
numSampPos = length(idxPos);

% Calculate PCA components pre-FMT, donor and post-FMT samples
xPca = zeros(nSamp,numCl+3);
yPca = zeros(nSamp,numCl+3);
xPca(1:numSampPre,1) = score(idxPre,1);
yPca(1:numSampPre,1) = score(idxPre,2);
xPca(1:numSampDon,2) = score(idxDon,1);
yPca(1:numSampDon,2) = score(idxDon,2);
xPca(1:numSampPos,3) = score(idxPos,1);
yPca(1:numSampPos,3) = score(idxPos,2);
for i = 1:numCl
    xPca(1:nSampCl(i),i+3) = score(idxSamp(1:nSampCl(i),i),1);
    yPca(1:nSampCl(i),i+3) = score(idxSamp(1:nSampCl(i),i),2);
end

% Plot PCA results
figure(55)
subplot(2,2,3)
scatter(xPca(1:numSampPre,1),yPca(1:numSampPre,1),175,'bs')
hold on
scatter(xPca(1:numSampDon,2),yPca(1:numSampDon,2),125,'ro')
scatter(xPca(1:numSampPos,3),yPca(1:numSampPos,3),175,'md')
text(xPca(1:nSampCl(1),4),yPca(1:nSampCl(1),4),'1','HorizontalAl','center','color','k')
text(xPca(1:nSampCl(2),5),yPca(1:nSampCl(2),5),'2','HorizontalAl','center','color','k')
hold off
xlabel("PC 1 (" + explained(1)+ "% variance)");
ylabel("PC 2 (" + explained(2)+ "% variance)")
legend('Pre-FMT','Donor','Post-FMT')

% Calculate geometric means of all samples using only sufficiently positive
% abundances
spAbd = abd16S';
numSamp = size(spAbd,1);
numTaxa = size(spAbd,2);
clr = zeros(numSamp,numTaxa);
for i = 1:numSamp
    xt1 = spAbd(i,:);
    xt2 = xt1(xt1>1e-3);
    nDim(i) = length(xt2);
    gmean(i) = prod(xt2)^(1/nDim(i));
    clr(i,:) = spAbd(i,:)/gmean(i);
end

% Find cluster with most pre-FMT samples
if sum(iCl(idxPre) == 1)/length(idxPre) > 0.5
    iClMax = 1;
    iClMin = 2;
else
    iClMax = 2;
    iClMin = 1;
end

% Calculate rho correlation values for all taxa pairs for all samples in
% high pre-FMT cluster
idxSet = [1:nSamp];
rho = zeros(numTaxa,numTaxa);
for i = 1:numTaxa
    for j = 1:numTaxa
        data1 = clr(idxSet,i);
        data2 = clr(idxSet,j);
        rho(i,j) = 1-var(data1-data2)/(var(data1)+var(data2));
    end
end

% Plot rho values values above a threshold for taxa of interest
figure(55)
subplot(2,2,2)
iSp = 2;
idxP = find(abs(rho(iSp,:)) > 0.05);
[rhos,idx] = sort(rho(iSp,idxP));
bar(rhos(1:end-1))
set(gca,'XTick',1:length(idxP)-1,'XTickLabel',xLabels1(idxP(idx(1:end-1))));
xtickangle(25)
ylabel("Proportionality with " + xLabels1(idxP(idx(end))));

% Find samples in clusters of interest
preFrac = sampTypeNum(:,1)./sum(sampTypeNum')';
[maxClVal,iClMax] = max(preFrac);
idxClMax = idxSamp(1:nSampCl(iClMax),iClMax);
[minClVal,iClMin] = min(preFrac);
idxClMin = idxSamp(1:nSampCl(iClMin),iClMin);

% Find significantly different NMPCs between chosen sample sets
for i = 1:size(Mp,1)
    numSet1 = dataMp(i,idxClMax);
    numSet2 = dataMp(i,idxClMin);
    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);
metSig = [];
for i = idxMet
    numSet1 = dataMp(i,idxClMax);
    numSet2 = dataMp(i,idxClMin);
    if fdrMet(i) < 0.05 && (mean(numSet1) > 50 || mean(numSet2) > 50) && (mean(numSet1)/mean(numSet2) < 0.5 || mean(numSet1)/mean(numSet2) > 2)
        metSig = [metSig; [i fdrMet(i)]];
    end
end
nSig = length(metSig);

% Generate Figure S11 in paper

% Plot significantly different NMPCs between chosen sample sets
figure(56)
sigIdx = metSig(:,1);
tempLabels = erase(textMp(sigIdx+2,1),"EX_");
tempLabels = erase(tempLabels,"[fe]");
xLabels3 = strrep(tempLabels,'_','-');
dataMax = mean(dataMp(sigIdx,idxClMax)')';
dataMin = mean(dataMp(sigIdx,idxClMin)')';
nPlt = round(nSig/2);
subplot(2,1,1)
bar([dataMax(1:nPlt) dataMin(1:nPlt)]);
set(gca, 'XTick', 1:nPlt, 'XTickLabel', xLabels3(1:nPlt));
xtickangle(45)
ylabel('Maximal production (mmol/day)')
legend("High pre-FMT cluster, " + nSampCl(iClMax) + " samples",...
    "Low pre-FMT cluster, " + nSampCl(iClMin) + " samples");
subplot(2,1,2)
bar([dataMax(nPlt+1:end) dataMin(nPlt+1:end)]);
set(gca, 'XTick', 1:length(xLabels3(nPlt+1:end)), 'XTickLabel', xLabels3(nPlt+1:end));
xtickangle(45)
ylabel('Maximal production (mmol/day)')

% 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}))-2;
end

% 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}))-2;
end

% Calculate mean and standard error of each BA metabolite for each cluster
dataMax1 = mean(dataMp(idxBA,idxClMax)')';
semMax1 = std(dataMp(idxBA,idxClMax)')'/sqrt(nSampCl(iClMax));
dataMin1 = mean(dataMp(idxBA,idxClMin)')';
semMin1 = std(dataMp(idxBA,idxClMin)')'/sqrt(nSampCl(iClMin));
meanData1 = [dataMax1 dataMin1];
semData1 = [semMax1 semMin1];

% Generate Figure S6 in paper

% Calculate mean and standard error of each BA metabolite for each cluster
figure(57)
subplot(2,1,1)
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 350])
tempLabels = erase(textMp(idxBA+2,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 pre-FMT cluster, " + nSampCl(iClMax) + " samples",...
    "Low pre-FMT cluster, " + nSampCl(iClMin) + " samples")
set(gca,'FontSize',12)

% Calculate mean and standard error of each AAA metabolite for each cluster
dataMax2 = mean(dataMp(idxAAA,idxClMax)')';
semMax2 = std(dataMp(idxAAA,idxClMax)')'/sqrt(nSampCl(iClMax));
dataMin2 = mean(dataMp(idxAAA,idxClMin)')';
semMin2 = std(dataMp(idxAAA,idxClMin)')'/sqrt(nSampCl(iClMin));
meanData2 = [dataMax1 dataMin1];
semData2 = [semMax1 semMin1];

% Calculate mean and standard error of each AAA metabolite for each cluster
subplot(2,1,2)
b2 = bar(meanData2,'grouped');
hold on
ngroups2 = size(meanData2,1);
nbars2 = size(meanData2,2);
groupwidth2 = min(0.8, nbars2/(nbars2 + 1.5));
for i = 1:nbars2
    x2 = (1:ngroups2) - groupwidth2/2 + (2*i-1) * groupwidth2 / (2*nbars2);
    errorbar(x2, meanData2(:,i), semData2(:,i), 'k', 'linestyle', 'none');
end
hold off
ylim([0 350])
tempLabels = erase(textMp(idxAAA+2,1),"EX_");
tempLabels = erase(tempLabels,"[fe]");
xLabels7 = strrep(tempLabels,'_','-');
set(gca, 'XTick', 1:nAAA, 'XTickLabel', xLabels7);
xtickangle(45)
ylabel('Maximal production (mmol/day)')
set(gca,'FontSize',12)