% estimateEffectSize
clear -except instruct

% % % Prompt user to load file containing combinedSummitStruct
baseDir = "C:\";

showPValues = 0;

% % Select data and analysis options
ctrlChoice = 'Generic';
[gfile,gdir] = uigetfile('*.xlsx',"Select input xlsx");
genoXlsx = [gdir,gfile];
[ifile,idir] = uigetfile('*.mat',"Select input struct");
inStruct = [idir,ifile];
outDir=gdir;
figWd = 10;
numBins = 5; % used to generate histograms of line categories later on
numSubPlots = 8;
showGenos = 0;
plotPVals = 0;
showCats = 1;
ymax = 1.5;
sortByES = 1; %sort by effect size, 1 = true, 0 = false

target = split(genoXlsx,' - ');
target = split(target(2),'.');
target = target(1);

load(inStruct{1});

genoTable = readtable(genoXlsx);

try
    mkdir(outDir);
catch
end

%ctrlChoice = 'Generic';  % 'Best' = best control, 'Cs' = Cs WF for everyone, 'Generic' = controls held constant across perturbation type
numGenos2Display = 0; % set as 0 to plot all genotypes, otherwise display # of genos indicated here
figHt = 4;  % figure height, in inches
minWd = figHt;

if ~numGenos2Display
    numGenos2Display = height(genoTable);
end

mkrSz = 10; % marker size for plotting mean

% % % Settings for bootstrapping w/ replacement
iterN = 1000;  % number of iterations for each E_C pair (# of t-tests performed)

[effectArray,pArray,p2Array] = deal(NaN(iterN,height(table)));
genoLabel = genoTable{:,'E_geno'};

% % Save categories for each line tested in catIds
catNames = genoTable.Properties.VariableNames(6:end-2);
catMarks = cell(1,numel(catNames));
for i=6:numel(catNames)+5
    catMarks{i-5} = genoTable{:,genoTable.Properties.VariableNames{i}};
end

nTable = NaN(height(genoTable),2);

for row=1:height(genoTable)
    
    % Grab geno names for E and C
    e_geno = genoTable{row,'E_geno'};
    
    switch ctrlChoice
        case 'Best'
            c_geno = genoTable{row,'Best_C'};
        case 'Cs'
            c_geno = {'Cs WF'}; % use Cs WF as control for everyone
        case 'Generic'
            c_geno = genoTable{row,'Generic_C'}; % uses generic control for each perturbation (e.g. mutants > w1118, TNT-E > Cs/TNT-E etc)
    end
        
    % Find corresponding indices in combinedSummit struct
    e_idx = findGeno(combinedSummit.geno,e_geno{1});
    c_idx = findGeno(combinedSummit.geno,c_geno{1});
    nTable(row,1) =  combinedSummit.numCad{e_idx};
    nTable(row,2) =  combinedSummit.numCad{c_idx};
    
    for j = 1:iterN
        % Bootstrap dtata
        e_bs = datasample(combinedSummit.sumInt{e_idx},combinedSummit.numCad{e_idx});
        c_bs = datasample(combinedSummit.sumInt{c_idx},combinedSummit.numCad{c_idx});
        
        % Perform 2-sided ttest on bootstrapped data
        [h,p,ci,stats] = ttest2(e_bs,c_bs);
        
        % Perform 2-sided ttest on original data (NOT bootstrapped)
        [h2,p2,ci2,stats2] = ttest2(combinedSummit.sumInt{e_idx},combinedSummit.sumInt{c_idx});
        
        % Save pvalues in pArray
        pArray(j,row) = p; % bootstrapped p value
        p2Array(j,row) = p2; % unbootstrapped p-value
        
        % Calculate effect size relative to best control
        e_mean = nanmean(e_bs);
        c_mean = nanmean(c_bs);
        effect = (e_mean-c_mean)/(c_mean);
        effectArray(j,row) = effect;
    end    
end

%% Generate violin plot of effect sizes of each geno on summiting

% % Find mean of t value for each geno and determine order sorted by this
% mean
means = nanmean(effectArray);
pmean = p2Array(1,:); %all values in array are the same, reflect two-sided t-test between indicated groups 

if sortByES
    [mVal,mIdx] = sort(means);
else
    [mVal,mIdx] = find(means);
end

fh = figure;
subplot(numSubPlots,1,1:5);
hold on;

% Calculate edges for each x-tile bin and plot dashed lines
yRange = -1.5:0.1:3.5;
edges = ones(1,numBins+1);
edges(length(edges)) = height(genoTable);
for i=1:numBins
    if i~=numBins
        edges(i+1) = round(height(genoTable)*(1/numBins)*i);
    end
    % %Plot dotted lines showing x-tiles
    plot(repmat(edges(i+1),length(yRange)),yRange,'--','Color',[0.7 0.7 0.7]);
end

% Plot horizontal line at 0
plot(0:length(means)+1,zeros(1,length(means)+2),'-','Color',[0.7 0.7 0.7]);

% Sort genotypes as desired (by listing in XLSX doc or by effect size
% (toggled w/ sortByES)
meansSort = means(mIdx);
pSort = pmean(mIdx);
disp(mIdx)
genoSort = genoLabel(mIdx);

% Save numerical results to text file, for incorporating into results
outAnalysis = strcat(baseDir,'analysis_',num2str(dataSet),'.txt');
fID = fopen(outAnalysis,'wt');
fprintf(fID,'Geno\tMeanES\tpVal\n');
for m =1:numel(meansSort)
    fprintf(fID,strcat(genoSort{m},'\t',num2str(meansSort(m)),'\t',num2str(pSort(m)),'\n'));
end
fclose(fID);

v = violinPlot(effectArray(:,mIdx),[0 0 0],fh);

ylabel('Effect size');

h = gca;

h.Clipping = 'off';
set(gca,'xtick',[],'xticklabel',[],'FontSize',12);
[h.YAxis.TickLength, h.XAxis.TickLength] = deal([0 0]);

xlim([0 numGenos2Display+1]);

ylim([-1.5 ymax]);

title(target);
xtickangle(90);

% % Plot p-values for comparisons on right y axis
if plotPVals
    yyaxis right;
    h2 = gca;

    pCutoff = 0.05;
    pOrdered = pmean(mIdx);
    pSig = find(pOrdered<pCutoff);
    pNonSig = find(pOrdered>=pCutoff);
    
    scatter(pSig,-1*log10(pOrdered(pSig)),'*','MarkerEdgeColor',[0 0 0]);
    disp(pSig);
    disp(pOrdered(pSig));
    scatter(pNonSig,-1*log10(pOrdered(pNonSig)),'*','MarkerEdgeColor',[0.7 0.7 0.7]);
    
    plot(0:length(pmean)+1,repmat(-1*log10(pCutoff),1,length(pmean)+2),':k');
    plot(0:length(pmean)+1,repmat(-1*log10(0.01),1,length(pmean)+2),':k');
    plot(0:length(pmean)+1,repmat(-1*log10(0.001),1,length(pmean)+2),':k');
    [h2.YAxis.TickLength, h2.XAxis.TickLength] = deal([0 0]);
    ylim([0 20]);
    ylabel('-log10(pvalue)');
    h2.YColor='black';
end

if showPValues
    pOrdered = pmean(mIdx);
end

if showGenos && ~showCats
    set(gca,'xtick',1:1:height(genoTable),'xticklabel',genoSort,'FontSize',6);
    xtickangle(90);
    xlabel('Genotype');
end

if showCats
    subplot(numSubPlots,1,6:8);
    hold on;
    
    % Plot dashed lines to indicate cutoff points between chosen percentiles (#
    % percentiles = numBins)
    yRange = 0:0.1:numel(catMarks);
    for i=1:numBins
        plot(repmat(edges(i+1),length(yRange)),yRange,'--','Color',[0.7 0.7 0.7]);
    end
    
    for i=1:numGenos2Display
        if mod(i,2)==0
            area([i-0.5 i+0.5],[numel(catMarks)+1 numel(catMarks)+1],'EdgeColor','none','FaceColor',[0.7 0.7 0.7],'FaceAlpha',0.1)
        end
    end
    
    % Specify color scheme
    color = hsv(numel(catNames)); % generate rainbow!
    
    % Initialize array to save which indices correspond to each category, so I
    % can use this info later on instead of having to rerun this loop
    saveCatHits = nan(height(genoTable),numel(catMarks));
    
    colorCt = 1;
    order = numel(catMarks);
    medCatRank = nan(1,numel(catMarks)); % initialize array to store average rank value for each category
    for i=1:numel(catMarks)
        catSort = catMarks{i}(mIdx);
        catHits = find(catSort>0);     % find genos that match this category (indicated in .xslx by 1)
        medCatRank(i) = median(catHits); % find mean rank for this category and store
        for j=1:length(catHits)
            if (numGenos2Display && catHits(j)<=numGenos2Display) || ~numGenos2Display
                plot([catHits(j)-0.25 catHits(j)+0.25],[order order],'Color',color(colorCt,:),'LineWidth',5)
            end 
        end
        if colorCt==length(color)
            colorCt=1;
        end
        colorCt = colorCt+1;
        order = order -1;
        saveCatHits(1:length(catHits),i) = catHits;
    end
    
    % Sort categories by median category rank (calculated in loop above)
    % Not currently using this for anything
    [b, catOrder] = sort(medCatRank);
    
    set(gca,'xtick',[],'ytick',1:1:numel(catMarks),'yticklabel',flip(catNames),'TickLabelInterpreter','None','FontSize',12);
    h2 = gca;
    [h2.YAxis.TickLength, h2.XAxis.TickLength] = deal([0 0]);
    ylim([0.5 numel(catMarks)+0.5]);
    xlim([0 numGenos2Display+1]);
    
    if showGenos
        set(gca,'xtick',1:1:height(genoTable),'xticklabel',genoLabel(mIdx),'FontSize',6);
        xtickangle(90);
        xlabel('Genotype');
    end
end

% Resize and save figure
[ts, ds] = Stamps;
set(gcf,'PaperUnits','inches','PaperPosition',[0 0 figWd figHt + (2/5)*figHt]);

imgOutBase = strcat(outDir,ds,'_',ts,'_',ctrlChoice,'_control_ES_iter=',num2str(iterN),'_','bins=',num2str(numBins),'_',target);

print(strcat(imgOutBase,'.png'),'-dpng');
print(strcat(imgOutBase,'.pdf'),'-dpdf');
close();

% --------------------------------------------------------------

% Initialize variable to hold counts of hits in each percentile, specified
% by numBins
[catHistoCounts, catHistoCountsNorm] = deal(nan(numel(catMarks),numBins));
    
% For each category, count hits in each percentile and save 
for i=1:numel(catMarks)
    % normalized counts per bin (divided by total number of lines with this
    % annotation)
    catHistoCountsNorm(i,:) = histcounts((saveCatHits(:,i)),edges)/sum(~isnan(saveCatHits(:,i)));
    
    %raw counts per bin
    catHistoCounts(i,:) = histcounts(saveCatHits(:,i),edges);

end

percentileHistos = transpose(catHistoCountsNorm);

% Plot one histogram for each category

figure;
hold on;
xrange = 0.5:0.1:numBins+0.5;

for j=1:numel(catNames)
    subplot(2,ceil(numel(catNames)/2),j);
    hold on;
    bar(catHistoCountsNorm(j,:),'FaceColor',color(j,:),'EdgeColor','none');
    ylim([0 0.8]);
    xlim([0.5 numBins+0.5]);
    plotText = strcat("N=",num2str(sum(~isnan(saveCatHits(:,j)))));
    text(numBins-1,0.7,plotText);
    title(catNames(j),'Interpreter','none');
    plot(xrange,repmat(1/numBins,1,length(xrange)),':k');
end

%% Find intersection between each category and PI
if dataSet == 1
    
    [overlapHist,overlapHistNorm] = deal(NaN(numel(catMarks)-1,numBins));
    for j=1:9
        % Find genotypes that have both annotation for idx j and PI (idx 10)
        hits = intersect(saveCatHits(:,j),saveCatHits(:,10));
        
        % Make histogram of where these indices are in terms of effect size
        overlapHist(j,:) = histcounts(hits,edges);
        overlapHistNorm(j,:) = overlapHist(j,:)./catHistoCounts(j,:);
    end
end

set(gcf,'PaperUnits','inches','PaperPosition',[0 0 figWd 0.5*figWd]);
print(strcat(imgOutBase,'_cat_histo.png'),'-dpng');
print(strcat(imgOutBase,'_cat_histo.pdf'),'-dpdf');
close();
