%% Alignment to reference figure
% Here we run all data against all yeasts... 
% Then we look at all chromosomes  Fig 2 and Suppl Fig 3
% Takes a collection of barcodes (bGen) 
% and compares it against a collection of ref barcodes (bT)
% you need the theory files already, can generate using hca tools
% this contains all for global and local alignment and plot 

%the uncommented code runs Fig 2
%% Load theory
% addpath(genpath('/Users/moacarlsson/Desktop/Yeast reference files'))
% fastas=dir('/Users/moacarlsson/Desktop/Yeast reference files/*fasta');
% addpath(genpath('C:/Users/luisma/OneDrive - Chalmers/yeastdata/ReferenceScer'))
% fastas=dir('C:/Users/luisma/OneDrive - Chalmers/yeastdata/ReferenceScer/*fasta');
%addpath(genpath('C:\Users\luisma\OneDrive - Chalmers\S.pombe genome proyect\Spombe_Genome\ref_fastas'))
%fastas=dir('C:\Users\luisma\OneDrive - Chalmers\S.pombe genome proyect\Spombe_Genome\ref_fastas\*fasta');

%for i=1:length(fastas)
%fastaFile{1,i}=strjoin({fastas(i).folder,'/',fastas(i).name},"");
%end

%quickly calculate theory for given single fasta file (fastaFile, nmbp, nmPerPx, psffac, islinear)
%[theoryStructRev_newpsf_scer,theoryStruct_newpsf_scer,bT_newpsf_scer] = prep_thry_for_local_comp(fastaFile, 0.25, 1,254,1); % I think stretch is closer to 0.2375 not 0.25

%% global alignment code only
% First we mix all data from cerevisiae and pombe
barcodeGen_mix = [barcodeGen_allmerged_scer barcodeGen_allmerged_fixed];

% we add a flag to keep track on them and see if they ennded up in the
% right place
for i=1:length(barcodeGen_mix)
    if i <= length(barcodeGen_allunmerged_scer)
    barcodeGen_mix{i}.original='scer';
    else
    barcodeGen_mix{i}.original='pombe';   
    end
end


theoryStruct=[theoryStruct_newpsf_scer
    theoryStruct_newpsf];


        import CBT.Hca.Import.import_hca_settings;
       
        [hcaSets] = import_hca_settings('hca_settings.txt');
    % change names to support previous hca
     hcaSets.bitmaskSettings=0;
     hcaSets.timestamp = datestr(clock(), 'yyyy-mm-dd_HH_MM_SS');
     hcaSets.theory.stretchFactors=0.85:0.01:1.15;
    %% bitmask repstuff
    for k=1:length(barcodeGen_mix)
     if isfield(barcodeGen_mix{1},'repmask')
     barcodeGen_mix{k}.rawBitmask= barcodeGen_mix{k}.rawBitmask & ~barcodeGen_mix{k}.repmask;
     barcodeGen_mix{k}.rawBitmask= barcodeGen_mix{k}.rawBitmask & barcodeGen_mix{k}.peakBitmask;
     end
    end
bitmask_buff=5;
barcodeGen_mix = expand_mask(barcodeGen_mix, bitmask_buff);
minsize=150;
    % size check
barcodeGen_mix_sizefiltered={};
l=0;
for i=1:length(barcodeGen_mix)

    if sum(diff(find(barcodeGen_mix{i}.rawBitmask==0))>=(minsize*max( hcaSets.theory.stretchFactors)))>=1  % This weird thing here makes sure all barcodes have a continuous bitmask thats longer that the minOverlap
    l=l+1;
    barcodeGen_mix_sizefiltered{l}=barcodeGen_mix{i};
    end

end
barcodeGen_mix=barcodeGen_mix_sizefiltered;

    %% compare theory to experiment
    barcodeGenC=barcodeGen_mix;
    import CBT.Hca.Core.Comparison.compare_distance;
    [rezMax,bestBarStretch,bestLength] = compare_distance(barcodeGenC,theoryStruct, hcaSets );
    % Extract data
    comparisonStructAll = rezMax;
    for i=1:length(comparisonStructAll)
        for j=1:length(bestBarStretch{i})
            comparisonStructAll{i}{j}.bestBarStretch = bestBarStretch{i}(j);
            comparisonStructAll{i}{j}.length = bestLength{i}(j);
        end
    end
    import CBT.Hca.Core.Comparison.combine_theory_results;
    [comparisonStruct] = combine_theory_results(theoryStruct, rezMax,bestBarStretch,bestLength);
   
%% Best chromosome extract
for n = 1:numel(theoryStruct)
 names{n}=theoryStruct(n).name;
 theoryStruct(n).rawBitmask=ones(length(theoryStruct(n).rawBarcode),1)';
 theoryStruct(n).rawBitmask(1:6)=0;
 theoryStruct(n).rawBitmask(end-6:end)=0;
 theoryStruct(n).rawBitmask=logical(theoryStruct(n).rawBitmask);
end
best_chromosomes=[];
for i = 1:numel(comparisonStruct)
      found = 0; % default if nothing matches
    for j = 1:numel(names)
        if strcmp(comparisonStruct{i}.name, names{j})
            found = j;
            break; % stop after first match
        end
    end   
    best_chromosomes(i) = found; 
end
%% hotizontal plot
    font_size=8;
cross_hca=1;
chr_to_plot=0;% 0 =all chromosomes
%chr_to_plot=1:16; % %for Scerevisiae
%chr_to_plot=17:19; %for pombe ones
makethemlines=1; % 1 for full , 2 for local
layerwidth=1;
threshold=-Inf;% add here for
unmask=0;
%whole_genome_coverage(barcodeGen, comparisonStruct, best_chromosomes, theoryStruct,layerwidth, chromosome,cross_hca,font_size,threshold,consensus_line,unmask,hcaSets)

 [f,chrconsensus,b,y] = whole_genome_coverage(barcodeGenC, comparisonStruct,best_chromosomes, theoryStruct, layerwidth,chr_to_plot,cross_hca,font_size,threshold,makethemlines,unmask,hcaSets);% other arguments are offset, chromosome (0=all), highlight (to highlight local or not ) 1/0
 set(gca,'Color','white')
 
 %% CCval histograms %we decided this was too many figs already ...unused
%  % Collect all CCvals
%  ccvals_all = cellfun(@(s) s.maxcoef(1), comparisonStruct);
%  ccdifs_all = ccvals_all - cellfun(@(s) s.maxcoef(2), comparisonStruct);
%  % Collect Og organism
%  names_all = cellfun(@(s) s.original, barcodeGen_mix, 'UniformOutput', false);
% 
%  toplot=ccvals_all;
%  % we got best_chromosomes already...
% 
%  % Let's assign color again...
%  colors_all= or(and(strcmp(names_all,'scer'),best_chromosomes<=16),and(strcmp(names_all, 'pombe'),best_chromosomes>16));
% 
% 
% figure
% tiledlayout(1,2)
% %-----ccvals
% edges = linspace(min(ccvals_all), max(ccvals_all), 30);
% nexttile
% hold on
% 
% histogram(ccvals_all(colors_all== 1), edges, ...
%     'FaceColor', '#0072B2', ...
%     'EdgeColor', 'none', ...
%     'FaceAlpha', 0.6);
% 
% histogram(ccvals_all(colors_all == 0), edges, ...
%     'FaceColor', '#E69F00', ...
%     'EdgeColor', 'none', ...
%     'FaceAlpha', 0.6);
% xlim([0,1])
% title('ccvals')
% hold off
% 
% 
% % -----ccdif
% % Common bin edges (important for fair comparison)
% edges = linspace(min(ccdifs_all), max(ccdifs_all), 30);
% 
% nexttile
% hold on
% 
% histogram(ccdifs_all(colors_all == 1), edges, ...
%     'FaceColor', '#0072B2', ...
%     'EdgeColor', 'none', ...
%     'FaceAlpha', 0.6);
% 
% histogram(ccdifs_all(colors_all == 0), edges, ...
%     'FaceColor', '#E69F00', ...
%     'EdgeColor', 'none', ...
%     'FaceAlpha', 0.6);
% xlim([0,1])
% title('ccdifs')
% hold off
% set(gca,'Color','white')
% 
% 
% %%
% % Common bin edges (important for fair comparison)
% edges = linspace(min(toplot), max(toplot), 30);
% 
% figure
% tiledlayout(2,1)
% 
% % ---- SCER ----
% nexttile
% hold on
% 
% idx_scer = strcmp(names_all, 'scer');
% 
% cc_scer = toplot(idx_scer);
% col_scer = colors_all(idx_scer);
% 
% histogram(cc_scer(col_scer == 1), edges, ...
%     'FaceColor', '#0072B2', ...
%     'EdgeColor', 'k', ...
%     'FaceAlpha', 0.6);
% 
% histogram(cc_scer(col_scer == 0), edges, ...
%     'FaceColor', '#E69F00', ...
%     'EdgeColor', 'k', ...
%     'FaceAlpha', 0.6);
% xlim([0,1])
% title('scer')
% hold off
% 
% 
% % ---- POMBE ----
% nexttile
% hold on
% 
% idx_pombe = strcmp(names_all, 'pombe');
% 
% cc_pombe = toplot(idx_pombe);
% col_pombe = colors_all(idx_pombe);
% 
% histogram(cc_pombe(col_pombe == 1), edges, ...
%     'FaceColor', '#0072B2', ...
%     'EdgeColor', 'k', ...
%     'FaceAlpha', 0.6);
% 
% histogram(cc_pombe(col_pombe == 0), edges, ...
%     'FaceColor', '#E69F00', ...
%     'EdgeColor', 'k', ...
%     'FaceAlpha', 0.6);
% xlim([0,1])
% title('pombe')
% hold off
% set(gca,'Color','white')




%% Local Reference Alignment

% From Compare consensuses

w =150;% minimum overlap length
per_stretch=5;

to_run="barcodes"; %barcodes only
bGen=barcodeGen_mix;
timestamp = datestr(clock(), 'yyyy-mm-dd_HH_MM_SS');
sF = (1-(per_stretch+8)/100):0.01:(1+per_stretch/100); 



%% Data Prep
% make bitmask logical and add the chromosome name and theory check
bT=[bT_newpsf_scer bT_newpsf_pombe];
for i=1:length(bT)
bT{i}.rawBitmask = logical(bT{i}.rawBitmask);
try
bT{i}.name=theoryStruct(i).name;
catch
end
bT{i}.theory=i;
end
% Expand mask
bitmask_buff=0;
bT = expand_mask(bT, bitmask_buff,1); % last argument adds mask (6px) on the edges




for i=1:length(bGen)
    bGen{i}.index=i;
end

if (to_run == "bargroups" || to_run == "merged" )
for i=1:length(cGenAll)

    cGenAll{i}.bargroupid=i;

end
mini_cGenAll=cGenAll(tokeep);
if isempty(mini_cGenAll)
mini_cGenAll=cGenAll;
end


if (to_run == "merged")
 [bGen_grouped,bGen_trimmed,allIdx,bGen_excluded]  = barcode_bargroup_cleaner(bGen, mini_cGenAll);
bGen=bGen_grouped;
else
bGen=mini_cGenAll;
end

end


for i = 1:numel(bGen)
    bm = bGen{i}.rawBitmask(:)';
    d = diff([0 bm 0]);
    runStarts = find(d == 1);
    runEnds   = find(d == -1) - 1;
    if isempty(runStarts)
        continue   % no 1s at all, empty barcode?
    end
    runLengths = runEnds - runStarts + 1;
    
    % Find longest run
    [~, idx] = max(runLengths);
    
    % Build new mask
    newMask = zeros(size(bm));
    newMask(runStarts(idx):runEnds(idx)) = 1;
    
    bGen{i}.rawBitmask = newMask;
end



l=0;
bGen_sizefiltered={};
numTheoryBar=0;
for i=1:length(bGen)

    if max(diff(find([0 bGen{i}.rawBitmask 0]==0)))>=(w*max(sF)) % This weird thing here makes sure all barcodes have a continuous bitmask thats longer that w 
    l=l+1;
    bGen_sizefiltered{l}=bGen{1,i};
    if isfield(bGen_sizefiltered{l},'theory')
    numTheoryBar=numTheoryBar+1;
    end
    end
end

cGen=bGen_sizefiltered;
numExpBar=length(cGen)-numTheoryBar;
% Mask peak all
% Reminder: bT has no bitmask!!


%cGen=mask_peaks(cGen,3.5);


% Warning: This could make already size checked bargroups smaller!
% Rep mask? Peak mask?
% for k=1:length(cGen)
%     if isfield(cGen{1},'repmask')
%     cGen{k}.rawBitmask= cGen{k}.rawBitmask & ~cGen{k}.repmask;
%     end
%     if isfield(cGen{1},'peakBitmask')
%     cGen{k}.rawBitmask= cGen{k}.rawBitmask & cGen{k}.peakBitmask;
%     end
% end


%% Can we do a different w for each cGen element? 
% Problem here is that theories are short!
% maybe add a check on mp_to_struct and make stuff NA if shorter than W
numExpBar=numel(bGen);
numTheoryBar=numel(bT);
template = struct( ...
    'score', [], ...
    'fullscore', [], ...
    'overlaplen', [], ...
    'pA', [], ...
    'pB', [], ...
    'or', [], ...
    'bestBarStretch', [], ...
    'lenA', [], ...
    'lenB', [], ...
    'partialScore', [], ...
    'partialLength', [], ...
    'h', [] );
minioS_mix_variable = repmat(template, length(bGen), numTheoryBar);
failedIdx = false(length(bGen),1);
tic
parfor c = 1:length(bGen)
%for c = bad'
try
tGen=[];
w=round(max(diff(find([0 bGen{c}.rawBitmask 0]==0)))/max(sF)*0.7);
tGen{1}=bGen{c};
tGen=[tGen bT];
[tempoS] = calc_overlap_mp(tGen,sF, w,timestamp);
fprintf('barcode: %d\n', c);
tempoS=tempoS(1,2:end);
if ~isempty(tempoS)
    minioS_mix_variable(c,:) = tempoS;
else
    failedIdx(c) = true;
end
catch ME
    warning('Error at c = %d', c);
    failedIdx(c) = true;
end
end
toc
bad = find(failedIdx); % if bnad is not empty, rerun uncomment line 358
clear tempoS

%% The process failed several times, this is how to look for what it missed

 missingRows = arrayfun(@(r) isempty(minioS_mix_variable(r,19).pB), ...
                       1:size(minioS_mix_variable,1));
missingIdx = find(missingRows);
%% Let's get some figures now
% first vs second first cc? %% Best match for each chromosome
% Must we filter out ccdif<0.05???
% how we get second max?
minioS=minioS_mix_variable;
unmask=0;
 
%localscore
alllocalscores= reshape([minioS(:).score], size(minioS,1),size(minioS,2)); 
alllocalscores=alllocalscores(:);
%globalscore
allglobalscores= reshape([minioS(:).fullscore], size(minioS,1),size(minioS,2));
allglobalscores=allglobalscores(:);
w = minioS(1,2).h;
alllenA = reshape([minioS(:).lenA], size(minioS,1),size(minioS,2)); 
alllenA = alllenA(:);
alllenA(alllenA==0) = nan;
alllenB = reshape([minioS(:).lenB], size(minioS,1),size(minioS,2)); 
alllenB(alllenB==0) = nan;
%partialscoresandlengths
allpartialScore= reshape([minioS(:).partialScore], size(minioS,1),size(minioS,2)); 
allpartialScore=allpartialScore(:);
allpartialLength= reshape([minioS(:).partialLength], size(minioS,1),size(minioS,2)); 
allpartialLength=allpartialLength(:);

% p-value for local
import Zeromodel.beta_ev_cdf;
pvalfun = @(x,l1,l2,nuF,w) 1-beta_ev_cdf(x,nuF*w,1,max(alphaN1*2*(max(l1,l2)-w),alphaN2*2*(l1-w+1)*(l2-w+1)),0);

k=length(alllocalscores);
minL = 50;
allpval = nan(1,k);

        for i=1:k
            allpval(i) = pvalfun(alllocalscores(i),alllenA(i),alllenB(i),alphaNu2,w);
        end
        %p-pvalue for single
        import Zeromodel.beta_cdf;
        pvalfunLeftover = @(x,l1,nuF) 1-beta_cdf(x,nuF*l1,1,0);
        allpvalLeftOver = nan(1,k);
        for i=1:k
            if allpartialLength(i) > minL
                allpvalLeftOver(i) = pvalfunLeftover(allpartialScore(i),allpartialLength(i),alphaNu1);
            end
        end
%         pval = [];
        allpvalCombined = (norminv(1-allpval)+norminv(1-allpvalLeftOver))/sqrt(2);
        allpvalCombined(isnan(allpvalCombined)) = (norminv(1-allpval(isnan(allpvalCombined))));
 

allpvalCombined(allpvalCombined==-Inf)=-8;
allpvalCombined(allpvalCombined==Inf)=8;
allscores=allpvalCombined;

% Fit a normal distribution
pd = fitdist(allscores(:), 'Normal');

% Get the value where P = 0.95 (i.e., the 95th percentile)
threshold = norminv(0.5, pd.mu, pd.sigma);
figure;
histogram(allscores, 'Normalization', 'pdf');  % 'pdf' scales it for overlay
hold on;

% Plot the fitted normal distribution
x_values = linspace(min(allscores), max(allscores), 100);
y_values = pdf(pd, x_values);
plot(x_values, y_values, 'r-', 'LineWidth', 2);

y_thresh = pdf(pd, threshold);
plot([threshold threshold], [0 y_thresh], 'k--', 'LineWidth', 1.5);

% Annotate
text(threshold, y_thresh, sprintf('  P=0.05 @ %.2f', threshold), ...
    'VerticalAlignment', 'bottom', 'Color', 'k');

xlabel('Value');
ylabel('Probability Density');
title('Histogram with Fitted Normal Distribution');
legend('Data Histogram', 'Fitted Normal PDF', '95th Percentile');

hold off;

%% Best is sorted according to local score (not stouffer)
% maybe ccdif no, but stscore or something else
ccdif_thres=0.03;
s=[];

stscores=[];
pvalLs=[];
pvalLOs=[];
Vals=[];
ccdifs=[];
best_chromosomes=[];
for r=1:height(minioS)
 [sortedVals,sortedIds,pscores,fullScores,overlaplen,partialScore,partialLength,pvalL,pvalLO,stscore]=sorted_scores(minioS(r,:),alphaN1,alphaN2,alphaNu1,alphaNu2);
 [~,maxind]=max(stscore);
 
 %how does this look?
 stscores(r)=max(stscore);
 pvalLs(r)=min(pvalL); 
 pvalLOs(r)=min(pvalLO);
Vals(r)=max(sortedVals);
ccdifs(r)=sortedVals(1)-sortedVals(2);
% WHAT TO DO HERE???
if ccdifs(r)<=ccdif_thres
best_chromosomes(r)= 0 ;% if ccdif < 0 kick them out!
else
best_chromosomes(r)=sortedIds(1); % best local score
%best_chromosomes(r)=sortedIds(maxind); %best stouffer score
end
end

 figure;
 histogram(ccdifs,10);  % 'pdf' scales it for overlay
 hold on;


font_size=8;
see_local=0; %0 no local, % 1local and global % 2 only local
chr_to_plot=0; %(0=all)
consensus_line=1; % 1 for full , 2 for local
layerwidth=1;
threshold=-Inf;%4.5;% -7.88
horizontalgap=600;
unmask=0;
cGen=[bGen bT];
%% horizontal plot all

 [f,chromosomeConsensus_bargroupsonly] = whole_genome_local_coverage(cGen, minioS,best_chromosomes,layerwidth,chr_to_plot,see_local,font_size,threshold,consensus_line,unmask,horizontalgap,alphaN1,alphaN2,alphaNu1,alphaNu2);% other arguments are offset, chromosome (0=all), highlight (to highlight local or not ) 1/0
%% hotizontal plot
    font_size=8;
cross_hca=1; 
coverageline=1;
chr_to_plot=1:16; %(0=all) scer
chr_to_plot=17:19; %(0=all) pombe
chr_to_plot=1:19; % all
makethemlines=1; 
layerwidth=1;
threshold=-Inf;
unmask=0;

%whole_genome_coverage(barcodeGen, comparisonStruct, best_chromosomes, theoryStruct,layerwidth, chromosome,cross_hca,font_size,threshold,consensus_line,unmask,hcaSets)
whole_genome_coverage(cGen, minioS,best_chromosomes, 0, layerwidth,chr_to_plot,cross_hca,font_size,threshold,makethemlines,unmask,coverageline);
ylim([-40 60]);
%% This is for Fig 3d
for t =1:length(thresholds)
 threshold=thresholds(t);
 [~,chrconsensusl,bl,yl] = whole_genome_coverage(cGen, minioS,best_chromosomes, 0, layerwidth,chr_to_plot,cross_hca,font_size,threshold,makethemlines,unmask,hcaSets);% other arguments are offset, chromosome (0=all), highlight (to highlight local or not ) 1/0
 set(gca,'Color','white')
 scer_whole_consensus=chrconsensusl{1}.forcedconsensus;
 pombe_whole_consensus=chrconsensusl{17}.forcedconsensus;
 for c = 2:16
 scer_whole_consensus=[scer_whole_consensus chrconsensusl{c}.forcedconsensus];
 end
 for p = 18:19
 pombe_whole_consensus=[pombe_whole_consensus chrconsensusl{p}.forcedconsensus];
 end
 scercov(t)=(1-sum(isnan(scer_whole_consensus))/length(scer_whole_consensus))*100;
 pombecov(t)=(1-sum(isnan(pombe_whole_consensus))/length(pombe_whole_consensus))*100;
end

%% ccdifs vs % of correct matches
% Thresholds
thresholds = 0:0.005:0.1;

% Preallocate
perc_scer  = zeros(size(thresholds));
perc_pombe = zeros(size(thresholds));

% Logical organism masks
is_scer  = strcmp(names_all, 'scer');
is_pombe = strcmp(names_all, 'pombe');

for k = 1:numel(thresholds)

    T = thresholds(k);

    % Only consider barcodes above threshold
    valid = ccdifs >= T;

    % --- SCER ---
    scer_valid = valid & is_scer;
    total_scer = sum(scer_valid);

    if total_scer > 0
        correct_scer = sum(colors_all(scer_valid) == 1);
        perc_scer(k) = 100 * correct_scer / total_scer;
    else
        perc_scer(k) = NaN;
    end

    % --- POMBE ---
    pombe_valid = valid & is_pombe;
    total_pombe = sum(pombe_valid);

    if total_pombe > 0
        correct_pombe = sum(colors_all(pombe_valid) == 1);
        perc_pombe(k) = 100 * correct_pombe / total_pombe;
    else
        perc_pombe(k) = NaN;
    end

end

figure; hold on

plot(thresholds, perc_scer,  '-o', 'LineWidth',2)
plot(thresholds, perc_pombe, '-s', 'LineWidth',2)

xlabel('ccdif threshold')
ylabel('% correctly matched')
legend('S. cerevisiae','S. pombe','Location','best')
grid on
%% betteeeer
% Thresholds
thresholds = 0.0:0.005:0.1;

% Preallocate
perc_scer      = zeros(size(thresholds));
perc_pombe     = zeros(size(thresholds));
retained_scer  = zeros(size(thresholds));
retained_pombe = zeros(size(thresholds));

% Logical masks
is_scer  = strcmp(names_all, 'scer');
is_pombe = strcmp(names_all, 'pombe');

total_scer_all  = sum(is_scer);
total_pombe_all = sum(is_pombe);

for k = 1:numel(thresholds)

    T = thresholds(k);
    valid = ccdifs > T;

    % ---- SCER ----
    scer_valid = valid & is_scer;
    total_scer = sum(scer_valid);

    retained_scer(k) = 100 * total_scer / total_scer_all;

    if total_scer > 0
        correct_scer = sum(colors_all(scer_valid) == 1);
        perc_scer(k) = 100 * correct_scer / total_scer;
    else
        perc_scer(k) = NaN;
    end

    % ---- POMBE ----
    pombe_valid = valid & is_pombe;
    total_pombe = sum(pombe_valid);

    retained_pombe(k) = 100 * total_pombe / total_pombe_all;

    if total_pombe > 0
        correct_pombe = sum(colors_all(pombe_valid) == 1);
        perc_pombe(k) = 100 * correct_pombe / total_pombe;
    else
        perc_pombe(k) = NaN;
    end

end
figure; hold on

yyaxis left
plot(thresholds, perc_scer,  '-o', 'LineWidth',2)
plot(thresholds, perc_pombe, '-s', 'LineWidth',2)
ylabel('% Correctly Matched')

yyaxis right
plot(thresholds, retained_scer,  '--o', 'LineWidth',2)
plot(thresholds, retained_pombe, '--s', 'LineWidth',2)
ylabel('% Retained')

xlabel('ccdif threshold')
legend({'Scer % correct', ...
        'Pombe % correct', ...
        'Scer % retained', ...
        'Pombe % retained'}, ...
        'Location','best')

grid on
%% with coverage
figure; hold on

plot(thresholds, perc_scer,  '-', 'LineWidth',2,'Color',[241/255 163/255 64/255])
plot(thresholds, perc_pombe, '-', 'LineWidth',2,'Color',[153/255 142/255 195/255])
ylabel('%')
ylim([0 100])

plot(thresholds, retained_scer,  '--', 'LineWidth',2,'Color',[241/255 163/255 64/255])
plot(thresholds, retained_pombe, '--', 'LineWidth',2,'Color',[153/255 142/255 195/255])

plot(thresholds, scercov,  ':', 'LineWidth',2,'Color',[241/255 163/255 64/255])
plot(thresholds, pombecov, ':', 'LineWidth',2,'Color',[153/255 142/255 195/255])

ylim([0 105])

xlabel('ccdif threshold')

legend({'Scer % correct', ...
        'Pombe % correct', ...
        'Scer % retained', ...
        'Pombe % retained', ...
        'Scer coverage', ...
        'Pombe coverage'}, ...
        'Location','best')

 %% CCval histograms
 % Collect all CCvals
 ccvals_all = Vals;
 ccdifs_all = ccdifs;
 % Collect Og organism
 names_all = cellfun(@(s) s.original, barcodeGen_mix, 'UniformOutput', false);
 
 toplot=ccdifs_all;
 % we got best_chromosomes already...

 % Let's assign color again...
 colors_all= or(and(strcmp(names_all,'scer'),best_chromosomes<=16),and(strcmp(names_all, 'pombe'),best_chromosomes>16));


figure
tiledlayout(1,2)
%-----ccvals
edges = linspace(min(ccvals_all), max(ccvals_all), 30);
nexttile
hold on

histogram(ccvals_all(colors_all== 1), edges, ...
    'FaceColor', '#0072B2', ...
    'EdgeColor', 'none', ...
    'FaceAlpha', 0.6);

histogram(ccvals_all(colors_all == 0), edges, ...
    'FaceColor', '#E69F00', ...
    'EdgeColor', 'none', ...
    'FaceAlpha', 0.6);
xlim([0,1])
title('ccvals')
hold off


% -----ccdif
% Common bin edges (important for fair comparison)
edges = linspace(min(ccdifs_all), max(ccdifs_all), 30);

nexttile
hold on

histogram(ccdifs_all(colors_all == 1), edges, ...
    'FaceColor', '#0072B2', ...
    'EdgeColor', 'none', ...
    'FaceAlpha', 0.6);

histogram(ccdifs_all(colors_all == 0), edges, ...
    'FaceColor', '#E69F00', ...
    'EdgeColor', 'none', ...
    'FaceAlpha', 0.6);
xlim([0,1])
title('ccdifs')
hold off
set(gca,'Color','white')


%%
% Common bin edges (important for fair comparison)
edges = linspace(min(toplot), max(toplot), 30);

figure
tiledlayout(2,1)

% ---- SCER ----
nexttile
hold on

idx_scer = strcmp(names_all, 'scer');

cc_scer = toplot(idx_scer);
col_scer = colors_all(idx_scer);

histogram(cc_scer(col_scer == 1), edges, ...
    'FaceColor', '#0072B2', ...
    'EdgeColor', 'k', ...
    'FaceAlpha', 0.6);

histogram(cc_scer(col_scer == 0), edges, ...
    'FaceColor', '#E69F00', ...
    'EdgeColor', 'k', ...
    'FaceAlpha', 0.6);
xlim([0,1])
title('scer')
hold off


% ---- POMBE ----
nexttile
hold on

idx_pombe = strcmp(names_all, 'pombe');

cc_pombe = toplot(idx_pombe);
col_pombe = colors_all(idx_pombe);

histogram(cc_pombe(col_pombe == 1), edges, ...
    'FaceColor', '#0072B2', ...
    'EdgeColor', 'k', ...
    'FaceAlpha', 0.6);

histogram(cc_pombe(col_pombe == 0), edges, ...
    'FaceColor', '#E69F00', ...
    'EdgeColor', 'k', ...
    'FaceAlpha', 0.6);
xlim([0,1])
title('pombe')
hold off
set(gca,'Color','white')


