% This script plots the scatter + contour of cell size vs. cell cycle, as
% well as the heatmap for deltaK48 over the cell cycle

load('Figure4CDE_Data');
chDAPI=1; chSE=2; chGem=3; chK48=4;

%% indexing 
idx=ResultTable.condition == "DMSO";
Gem=ResultTable.NInt(idx,chGem);
DNA=ResultTable.NInt(idx,chDAPI);
SE=ResultTable.CInt(idx,chSE);
K48=ResultTable.CInt(idx,chK48);
lGem=log(abs(Gem));
lGem(isinf(lGem))=NaN;
lGem(isnan(lGem))=min(lGem);

DNA=DNA(SE<prctile(SE,98));lGem=lGem(SE<prctile(SE,98));
K48=K48(SE<prctile(SE,98));SE=SE(SE<prctile(SE,98));

%% indexing 
idx=ResultTable.condition == "CFZ";
GemCFZ=ResultTable.NInt(idx,chGem);
DNACFZ=ResultTable.NInt(idx,chDAPI);
SECFZ=ResultTable.CInt(idx,chSE);
K48CFZ=ResultTable.CInt(idx,chK48);
lGemCFZ=log(abs(GemCFZ));
lGemCFZ(isinf(lGemCFZ))=NaN;
lGemCFZ(isnan(lGemCFZ))=min(lGemCFZ);

DNACFZ=DNACFZ(SECFZ<prctile(SECFZ,98));lGemCFZ=lGemCFZ(SECFZ<prctile(SECFZ,98));
K48CFZ=K48CFZ(SECFZ<prctile(SECFZ,98));SECFZ=SECFZ(SECFZ<prctile(SECFZ,98));

%% normlization by control
K48CFZ=K48CFZ./nanmedian(K48);
K48=K48./nanmedian(K48);

SECFZ=SECFZ./nanmedian(SE);
SE=SE./nanmedian(SE);

%% normalizing DNA & Geminin
[~,~,h]=ksdensity(DNA);
[f_dna,s,~]=ksdensity(DNA,'width',0.3*h);
[mx,mxi]=max(f_dna);
DNA=DNA/s(mxi);

[~,~,h]=ksdensity(DNACFZ);
[f_dna,s,~]=ksdensity(DNACFZ,'width',0.3*h);
[mx,mxi]=max(f_dna);
DNACFZ=DNACFZ/s(mxi);

[~,~,h]=ksdensity(lGem);
[f_gem,s,h]=ksdensity(lGem,'width',h*0.6);
[pks,pksi]=findpeaks(f_gem,'MinPeakProminence',0.1);
lGem=(lGem-s(pksi(1)))./(s(pksi(2))-s(pksi(1)))+1; %change accdgly

[~,~,h]=ksdensity(lGemCFZ);
[f_gemCFZ,s,h]=ksdensity(lGemCFZ,'width',h*0.6);
[pks,pksi]=findpeaks(f_gemCFZ,'MinPeakProminence',0.1);
lGemCFZ=(lGemCFZ-s(pksi(1)))./(s(pksi(2))-s(pksi(1)))+1; %change accdgly

%% normalizing K48 
f_k48=fit(SE,K48,'poly1','robust','bisquare');
K48=K48-f_k48(SE);
K48=K48-min(K48); 

f_k48CFZ=fit(SECFZ,K48CFZ,'poly1','robust','bisquare');
K48CFZ=K48CFZ-f_k48CFZ(SECFZ);
K48CFZ=K48CFZ-min(K48CFZ); 

%% checking DNA and Gem 
ii=DNA>0.7 & DNA<2.2;
del=(DNA>1.12 & lGem<1.5); %del=DNA>1.12 & lGem<1.5;
DNA=DNA(ii&~del); lGem=lGem(ii&~del);
SE=SE(ii&~del); K48=K48(ii&~del);
% figure; scatterhist(DNA,lGem); title('DMSO')

ii=DNACFZ>0.7 & DNACFZ<2.3;
del1=(DNACFZ>1.15 & lGemCFZ<1.43); del2=(DNACFZ>1.8 & lGemCFZ<1.6);
del=del1|del2;
DNACFZ=DNACFZ(ii&~del); lGemCFZ=lGemCFZ(ii&~del);
SECFZ=SECFZ(ii&~del); K48CFZ=K48CFZ(ii&~del);
% figure; scatterhist(DNACFZ,lGemCFZ); title('DMSO+CFZ')

%% creating Ldist from DNA and Gem
L=exp(linspace(log(0.7),log(100),100)); % log(1)
D=@(n,K,L)((L).^n./(K^n+(L).^n)+1);
G=@(K,L)(log(L)/log(K)+1);
FitFunc=@(X) Adjust_DNA_Gem(DNA,lGem,X);
figure; plot(DNA,lGem,'.'); hold on; title('DMSO')
pr=fminsearch(FitFunc,[1.5 25 40]); %[2,31,40] [2,31,200]
n=pr(1); P1=pr(2); P2=pr(3);
plot(D(n,P1,L),G(P2,L),'.r')
Parametrised=[D(n,P1,L)' G(P2,L)'];
RealD=[DNA lGem];
Dst=pdist2(Parametrised,RealD);
[mn,mni]=min(Dst); ii=DNA>1.5 & lGem<1.5; mni(ii)=max(mni)+10; LL=mni;
xlim([.5 3]); ylim([0.5 2.5]) 

FitFunc=@(X) Adjust_DNA_Gem(DNACFZ,lGemCFZ,X);
figure; plot(DNACFZ,lGemCFZ,'.'); hold on; title('DMSO+CFZ')
pr=fminsearch(FitFunc,[1.4 35 40]); %[2,31,40] [2,31,200]
n=pr(1); P1=pr(2); P2=pr(3);
plot(D(n,P1,L),G(P2,L),'.r')
Parametrised=[D(n,P1,L)' G(P2,L)'];
RealD=[DNACFZ lGemCFZ];
Dst=pdist2(Parametrised,RealD);
[mn,mni]=min(Dst); ii=DNACFZ>1.5 & lGemCFZ<1.5; mni(ii)=max(mni)+10; LLCFZ=mni;
xlim([.5 3]); ylim([0.5 2.5]) 

%% heatmap of cell size over cell cycle
% figure;plot(LL,log2(SE),'o'); title('DMSO')
gridx1 = linspace(0,100,100);
gridx2 = linspace(-1.5,1.5,70);  %linspace(23.5,26.5,70); 
[x1,x2] = meshgrid(gridx1, gridx2);
x1v = x1(:); x2v = x2(:);
xi = [x1v x2v];
f=ksdensity([LL',log2(SE)],xi,'width',[2 0.12]);
f=reshape(f,[length(gridx2),length(gridx1)]);
g=@(p) (0.98-sum(f(f>p))/sum(f(:))).^2;
g0=fminsearch(g,0.0025); f(f<g0)=NaN;
% figure();clf;imagesc(log2(f)); set(gca, 'YDir','normal');title('DMSO')

% figure;plot(LLCFZ,log2(SECFZ),'o');title('DMSO+CFZ')
fCFZ=ksdensity([LLCFZ',log2(SECFZ)],xi,'width',[2 0.12]);
fCFZ=reshape(fCFZ,[length(gridx2),length(gridx1)]);
g=@(pCFZ) (0.98-sum(fCFZ(fCFZ>pCFZ))/sum(fCFZ(:))).^2;
g0=fminsearch(g,0.0025);
fCFZ(fCFZ<g0)=NaN;
% figure();clf;imagesc(log2(fCFZ)); set(gca, 'YDir','normal');title('DMSO+CFZ')
%% PLOT scatter with contour of log(SE) vs. Ldist -- FIgure 4D
figure; plot(LL,log2(SE),'o','markerfacecolor',[0.7 0.7 0.7],'markeredgecolor','non','markersize',3); hold on
contour(x1,x2,log2(f),'k','linewidth',1.5); axis([0 100 -1.5 1.5]);%axis([0 100 23.3 26.7]);
xlabel('Cell cycle stage (a.u.)'); ylabel('log_2(Cell size) (a.u.)');
set(gca,'xtick',0:20:100,'ytick',-1.5:.5:1.5,'fontsize',13); title('DMSO')

%% heatmap of Ldist vs. K48
f_k48=fit([LL' log2(SE)],K48,'lowess','span',0.2);
f_k48_d=f_k48(x1,x2);
f_k48_d(isnan(f))=NaN;
% figure; pcolor(x1,x2,f_k48_d)
% shading interp; colormap jet; colorbar; title('DMSO')
 
f_k48CFZ=fit([LLCFZ' log2(SECFZ)],K48CFZ,'lowess','span',0.2);
f_k48_d_CFZ=f_k48CFZ(x1,x2);
f_k48_d_CFZ(isnan(fCFZ))=NaN;
% figure; pcolor(x1,x2,f_k48_d_CFZ)
% shading interp; colormap jet; colorbar; title('DMSO+CFZ')

% save DMSO variables
f_DMSO=f;
f_k48_d_DMSO=f_k48_d;
% figure; % sanity check
% pcolor(x1,x2,f_k48_d_DMSO); shading interp; colormap jet; colorbar; caxis([.7 1.05]);
% hold on; contour(x1,x2,log2(f_DMSO),'k','linewidth',1.5); title('DMSO')

% subtract DMSO from DMSO+CFZ -- Figure 4E
f_k48_delta=f_k48_d_CFZ-f_k48_d;
figure; pcolor(x1,x2,f_k48_delta); shading interp; colormap jet; colorbar; caxis([2.2 2.73]);
hold on; contour(x1,x2,log(f_DMSO),'k','linewidth',1.5); 
axis([0 100 -1.5 1.5]); xlabel('Cell Cycle Stage'); ylabel('log_2(Cell Size) (a.u.)') %axis([0 100 23.5 26.5]);
set(gca,'fontsize',13,'xtick',0:20:100,'ytick',-1.5:.5:1.5);

%%
function [Dst]=Adjust_DNA_Gem(DNA,lGem,X)

n=X(1);
P1=X(2);
P2=X(3);
% P2=100;
L=exp(linspace(log(1),log(100),100));
D=@(n,K,L)((L).^n./(K^n+(L).^n)+1);
G=@(K,L)(log(L)/log(K)+1);
% % figure(7);clf;
% % plot(DNA,lGem,'.')
% % hold on
% % n=3; %5
% % P1=30; %20
% % P2=60; %30
% % plot(D(n,P1,L),G(P2,L),'or')
% % CERYL: see if you can find values for P1, P2 and n that optimize the fit
% %%
Parametrised=[D(n,P1,L)' G(P2,L)'];
RealD=[DNA lGem];
Dst=pdist2(Parametrised,RealD);
Dst=min(Dst);
Dst=Dst(Dst<0.2 & (DNA<1.85)');
Dst=mean(Dst);
% Dst=X(1)+X(2)+X(3);
end