load('Figure4AB_Data.mat') % same as Figure3BCD_Data.mat

%% Reorganize data by treatments
Tt = [1 8 24]; 
cellcycle = ["All cells" "G1" "G1S" "S" "G2+M"];

conditions = table; 
conditions.Drug = ["DMSO"; "SNS"];
NN = size(conditions,1); Nt = length(Tt);
conditions.deltaIF = cell(NN,Nt); % conditions by timepoints
conditions.deltaIF_mean = nan(NN,length(cellcycle),Nt);
conditions.deltaIF_CI = nan(NN,length(cellcycle),Nt);
conditions.cellsize = cell(NN,Nt);
conditions.cellsize_mean = nan(NN,length(cellcycle),Nt);
conditions.cellsize_CI = nan(NN,length(cellcycle),Nt);
conditions.cellnum = cell(NN,Nt);
for nt = 1:Nt
    for nn = 1:NN
        temp = statTab(statTab.Drug==conditions.Drug(nn) & statTab.DrugTime==Tt(nt) & statTab.AbAdded==1,:);
        idx = temp.CFZ==1;
        % calculate delta_K48 by substracting CFZ+ & CFZ-
        nrep = sum(idx);
        temp1IF = temp.MedianIF(idx,1:end-1); temp2IF = temp.MedianIF(~idx,1:end-1); 
        temp1size = temp.MedianSize(idx,1:end-1); temp2size = temp.MedianSize(~idx,1:end-1);
%         conditions.deltaIF{nn, nt} = temp1IF-temp2IF;
%         conditions.cellsize{nn, nt} = temp1size; %(temp1size+temp2size)/2;
        for cc = 1:length(cellcycle)
            temp1c = repmat(temp1IF(:,cc), nrep, 1);
            temp2c = repmat(temp2IF(:,cc), 1, nrep); temp2c = temp2c'; temp2c = temp2c(:);
            conditions.deltaIF{nn, nt}(:,cc) = temp1c-temp2c;
            temp1c = repmat(temp1size(:,cc), nrep, 1);
            temp2c = repmat(temp2size(:,cc), 1, nrep); temp2c = temp2c'; temp2c = temp2c(:);
            conditions.cellsize{nn, nt}(:,cc) = (temp1c+temp2c)/2;
        end
        conditions.deltaIF_mean(nn,:,nt) = mean(conditions.deltaIF{nn, nt}); 
        conditions.deltaIF_CI(nn,:, nt) = tinv(.9725, nrep-1)*std(conditions.deltaIF{nn, nt})/sqrt(nrep);%95% confidence interval range
        
        conditions.cellsize_mean(nn,:,nt) = mean(conditions.cellsize{nn, nt});
        nrep = size(conditions.cellsize{nn, nt},1);
        conditions.cellsize_CI(nn,:, nt) = tinv(.9725, nrep-1)*std(conditions.cellsize{nn, nt})/sqrt(nrep); %95% confidence interval range
        
        conditions.cellnum{nn, nt} = temp.CellNum(:,1); 
    end
    normFactorIF = mean(conditions.deltaIF{1, nt}(:,2)); % G1 of controls
    normFactorCellSize = mean(conditions.cellsize{1, nt}(:,2)); % G1 of controls
    
    conditions.deltaIF_mean(:,:, nt) = conditions.deltaIF_mean(:,:, nt)/normFactorIF;
    conditions.deltaIF_CI(:,:, nt) = conditions.deltaIF_CI(:,:, nt)/normFactorIF;
    conditions.cellsize_mean(:,:, nt) = conditions.cellsize_mean(:,:, nt)/normFactorCellSize;
    conditions.cellsize_CI(:,:, nt) = conditions.cellsize_CI(:,:, nt)/normFactorCellSize;
    for nn = 1:NN
        conditions.deltaIF{nn, nt} = conditions.deltaIF{nn, nt}/normFactorIF;
        conditions.cellsize{nn, nt} = conditions.cellsize{nn, nt}/normFactorCellSize;
    end
end

clear ans idx nn nt nrep temp temp1 temp2 normFactorCellSize normFactorIF temp1c temp2c  temp1IF temp2IF temp1size temp2size

%% over cell cycle
figure; hold on
for cc = 2:length(cellcycle)
    temp = conditions.cellsize{1,3}(:,cc);
    scatter(cc+randn(length(temp),1)*0.05, temp, 3, 'b', 'filled', 'MarkerFaceAlpha', .5)
    temp = conditions.cellsize{2,3}(:,cc);
    scatter(cc+randn(length(temp),1)*0.05, temp, 3, 'r', 'filled', 'MarkerFaceAlpha', .5)
end
errorbar(2:length(cellcycle), conditions.cellsize_mean(1,2:5,3), conditions.cellsize_CI(1,2:5,3), 'bs')
errorbar(2:length(cellcycle), conditions.cellsize_mean(2,2:5,3), conditions.cellsize_CI(2,2:5,3), 'rs')
fitobj = fit([2 4 5]', conditions.cellsize_mean(1,[2 4 5],3)','poly1');
plot([1.8 5.2], fitobj.p1*[1.8 5.2]+fitobj.p2,'k')
fitobj = fit([2 4 5]', conditions.cellsize_mean(2,[2 4 5],3)','poly1');
plot([1.8 5.2], fitobj.p1*[1.8 5.2]+fitobj.p2,'k')
title("24h"); xlim([1.5 5.5]); ylim([0.9 1.9]);
xticks(2:length(cellcycle)); xticklabels(cellcycle(2:end)); 
ylabel('Cell Size (a.u.)');


figure; hold on
for cc = 2:length(cellcycle)
    temp = conditions.deltaIF{1,3}(:,cc);
    scatter(cc+randn(length(temp),1)*0.05, temp, 3, 'b', 'filled', 'MarkerFaceAlpha', .5)
    temp = conditions.deltaIF{2,3}(:,cc);
    scatter(cc+randn(length(temp),1)*0.05, temp, 3, 'r', 'filled', 'MarkerFaceAlpha', .5)
end
errorbar(2:length(cellcycle), conditions.deltaIF_mean(1,2:5,3), conditions.deltaIF_CI(1,2:5,3), 'bs')
errorbar(2:length(cellcycle), conditions.deltaIF_mean(2,2:5,3), conditions.deltaIF_CI(2,2:5,3), 'rs')
fitobj = fit([2 4 5]', conditions.deltaIF_mean(1,[2 4 5],3)','poly1');
plot([1.8 5.2], fitobj.p1*[1.8 5.2]+fitobj.p2,'k')
fitobj = fit([2 4 5]', conditions.deltaIF_mean(2,[2 4 5],3)','poly1');
plot([1.8 5.2], fitobj.p1*[1.8 5.2]+fitobj.p2,'k')
title("24h"); xlim([1.5 5.5]); ylim([0.8 2.3]);
xticks(2:length(cellcycle)); xticklabels(cellcycle(2:end)); 
ylabel('delta-K48 polyUb');
