load('Figure1_Data.mat');
        
%% SYNCH to birth or G1/S
chGem=1; chArea=2;
ctrl = find(fulltracks.Treatment=="DMSO");
sns = find(fulltracks.Treatment=="SNS");
for nn=1:height(fulltracks)
    L=size(fulltracks.FullTracks{nn,1},1); 
    clear GemRiseVec
    SynchToStart{nn}=ones(L,300)*NaN;
    SynchTo_G1S{nn}=ones(L,600)*NaN;
    SynchTo_G1S_GR{nn}=ones(L,600)*NaN;
    SynchTo_G1S_GR_SL{nn}=ones(L,600)*NaN;
    SynchToStart_G1{nn}=false(L,300);
    for ii=1:L
        kk=ii;
        xdata=fulltracks.FullTracks{nn,1}(kk).Feat.Ti-fulltracks.FullTracks{nn,1}(kk).Feat.Ti(1)+1;
        ydata=fulltracks.FullTracks{nn,1}(kk).Feat.Area;
        zdata=fulltracks.FullTracks{nn,1}(kk).GrowthRate.GR;
        zzdata=fulltracks.FullTracks{nn,1}(kk).GrowthRateSL;
        xdata=xdata(4:end-4); ydata=ydata(4:end-4); 
        gemrise=round(fulltracks.FullTracks{nn,1}(kk).Gemrise);
        if length(zdata)>length(ydata) && gemrise>5 && size(ydata,1)>gemrise
            zdata=zdata(1:length(ydata));
            GemRiseVec(ii)=gemrise;
            L2=length(ydata); L3=length(zzdata);
            SynchTo_G1S_GR{nn}(ii,(302-gemrise):300)=zdata(1:(gemrise-1)); % GR from birth to 300 (G1/S)
            SynchTo_G1S_GR{nn}(ii,301:(301+L2-gemrise))=zdata(gemrise:end); % GR from 301 (G1/S) to end
            SynchTo_G1S_GR_SL{nn}(ii,(302-gemrise):300)=zzdata(1:(gemrise-1)); % GR from birth to 300 (G1/S)
            SynchTo_G1S_GR_SL{nn}(ii,301:(301+L3-gemrise))=zzdata(gemrise:end); % GR from 301 (G1/S) to end
        end
        if gemrise>5 && size(ydata,1)>gemrise
            GemRiseVec(ii)=gemrise;
            L2=length(ydata);
            SynchToStart_G1{nn}(ii,1:gemrise)=true(1,gemrise); % logical: 1 for G1, 0 for after
            SynchToStart{nn}(ii,1:L2)=ydata'; % nuclear area from birth forward
            SynchTo_G1S{nn}(ii,(302-gemrise):300)=ydata(1:(gemrise-1)); % area from birth to 300 (G1/S)
            SynchTo_G1S{nn}(ii,301:(301+L2-gemrise))=ydata(gemrise:end); % area from 301 (G1/S) to end
            SynchToStart1=SynchToStart{nn};
            SynchToStart2=SynchToStart{nn};
            SynchToStart1(~SynchToStart_G1{nn})=NaN; % make cells not in G1 into nan
            SynchToStart2(SynchToStart_G1{nn})=NaN; % make cells in G1 into nan
        end
    end  
end
%% PLOT FIGURES
% calculate average G1 length
for tt=[ctrl sns]
    g1=[];
    for kk=1:size(fulltracks.FullTracks{tt,1},1)
        g1=[g1;fulltracks.FullTracks{tt,1}(kk).G1];
    end
    if tt==ctrl
        g1_ctrl=g1;
    elseif tt==sns
        g1_sns=g1;
    end
end
g1_ctrl_CI=bsxfun(@times, std(g1_ctrl)/sqrt(length(g1_ctrl)), tinv(0.975, length(g1_ctrl)-1));
g1_sns_CI=bsxfun(@times, std(g1_sns)/sqrt(length(g1_sns)), tinv(0.975, length(g1_sns)-1));

% calculate average growth rate
for tt=[ctrl sns]
    gr=[];
    for kk=1:size(fulltracks.FullTracks{tt,1},1)
        gr=[gr;fulltracks.FullTracks{tt,1}(kk).AvgGRSL];
    end
    if tt==ctrl
        gr_ctrl=gr;
    elseif tt==sns
        gr_sns=gr;
    end
end
gr_ctrl_CI=bsxfun(@times, std(gr_ctrl)/sqrt(length(gr_ctrl)), tinv(0.975, length(gr_ctrl)-1));
gr_sns_CI=bsxfun(@times, std(gr_sns)/sqrt(length(gr_sns)), tinv(0.975, length(gr_sns)-1));

% plot bar charts of avg. G1 length and Growth Rate 
figure; b=bar([1 2],[mean(g1_ctrl) mean(g1_sns)],0.75); b.FaceColor='flat'; b.CData=[1 1 1]; hold on;
errorbar([1 2],[mean(g1_ctrl) mean(g1_sns)],[g1_ctrl_CI g1_sns_CI],'color','k','linestyle','none','linewidth',2);
scatter(normrnd(1,0.1,[length(g1_ctrl),1]),g1_ctrl,5,[0.5 0.5 0.5],'filled','markerfacealpha',0.5)
scatter(normrnd(2,0.1,[length(g1_sns),1]),g1_sns,5,[0.5 0.5 0.5],'filled','markerfacealpha',0.5)
set(gca,'xtick',1:2,'xticklabel',{'Ctrl','SNS'},'fontsize',14,'ytick',0:10:30); ylim([0 30]); ylabel('G1 Length (hr)')
[~, pval] = ttest2(g1_ctrl, g1_sns); disp(pval)

figure; b=bar([1 2],[mean(gr_ctrl) mean(gr_sns)],0.75); b.FaceColor='flat'; b.CData=[1 1 1]; hold on;
errorbar([1 2],[mean(gr_ctrl) mean(gr_sns)],[gr_ctrl_CI gr_sns_CI],'color','k','linestyle','none','linewidth',2);
scatter(normrnd(1,0.1,[length(gr_ctrl),1]),gr_ctrl,5,[0.5 0.5 0.5],'filled','markerfacealpha',0.5)
scatter(normrnd(2,0.1,[length(gr_sns),1]),gr_sns,5,[0.5 0.5 0.5],'filled','markerfacealpha',0.5)
set(gca,'xtick',[1 2],'xticklabel',{'Ctrl','SNS'},'fontsize',14,'ytick',0:4:12); ytickformat('%g'); ylim([0 12])
ylabel('Avg. Growth Rate (pix./hr)')
[~, pval] = ttest2(gr_ctrl, gr_sns); disp(pval)

% mean and CI nuclear area synch to G1/S (threshold of >10 cells)
avgSynchG1S=nan(4,size(SynchTo_G1S{1,1},2));
CISynchG1S=nan(4,size(SynchTo_G1S{1,1},2));
for tt=[ctrl sns]
    for kk=1:size(SynchTo_G1S{1,tt},2)
        if sum(~isnan(SynchTo_G1S{1,tt}(:,kk)))>20
            temp=SynchTo_G1S{1,tt}(:,kk); temp=temp(~isnan(temp));
            avgSynchG1S(tt,kk)=nanmean(temp);
            CISynchG1S(tt,kk)=bsxfun(@times, std(temp)/sqrt(length(temp)), tinv(0.975, length(temp)-1));
        else
            continue
        end
    end
end

figure; hold on
xx=find(~isnan(avgSynchG1S(ctrl,:))); y=avgSynchG1S(ctrl,:); yy=y(xx); z=CISynchG1S(ctrl,:); zz=z(xx);
smyy = smooth(xx, yy, 10)';
plot(xx, smyy,'color',[0 83 125]./255,'linewidth',2)
fill([xx xx(end:-1:1)], [smyy+zz smyy(end:-1:1)-zz(end:-1:1)], [0 83 125]./255, 'FaceAlpha',.5, 'EdgeColor','none')
xx=find(~isnan(avgSynchG1S(sns,:))); y=avgSynchG1S(sns,:); yy=y(xx); z=CISynchG1S(sns,:); zz=z(xx);
smyy = smooth(xx, yy, 10)';
plot(xx, smyy,'color',[194 87 87]./255,'linewidth',2)
fill([xx xx(end:-1:1)], [smyy+zz smyy(end:-1:1)-zz(end:-1:1)], [194 87 87]./255, 'FaceAlpha',.5, 'EdgeColor','none')
line([300 300],get(gca,'ylim'),'color','k','linewidth',1.5);axis([225 365 430 970])
xlabel('hr from G1/S');ylabel('Nuclear Area (pix.)'); 
set(gca,'xtick',[240 300 360],'xticklabel',[-15 0 15],'fontsize',14,'ytick',500:200:900); box on
line([300-mean(g1_ctrl)*4 300-mean(g1_ctrl)*4],[430 970],'linestyle','-','color',[0 83 125]./255,'linewidth',1.5);
line([300-mean(g1_sns)*4 300-mean(g1_sns)*4],[430 970],'linestyle','-','color',[194 87 87]./255,'linewidth',1.5);

% mean and CI nuclear area synch to birth (threshold of >10 cells)
avgSynchB=nan(4,size(SynchToStart{1,1},2));
CISynchB=nan(4,size(SynchToStart{1,1},2));
for tt=[ctrl sns]
    for kk=1:size(SynchToStart{1,tt},2)
        if sum(~isnan(SynchToStart{1,tt}(:,kk)))>20
            temp=SynchToStart{1,tt}(:,kk); temp=temp(~isnan(temp));
            avgSynchB(tt,kk)=nanmean(temp);
            CISynchB(tt,kk)=bsxfun(@times, std(temp)/sqrt(length(temp)), tinv(0.975, length(temp)-1));
        else
            continue
        end
    end
end

figure; hold on % stopping after G1/S
xx=find(~isnan(avgSynchB(ctrl,:))); y=avgSynchB(ctrl,:); yy=y(xx); z=CISynchB(ctrl,:); zz=z(xx);
smyy = smooth(xx, yy, 10)';
plot(xx, smyy,'color',[0 83 125]./255,'linewidth',2)
fill([xx xx(end:-1:1)], [smyy+zz smyy(end:-1:1)-zz(end:-1:1)], [0 83 125]./255, 'FaceAlpha',.5, 'EdgeColor','none')
xx=find(~isnan(avgSynchB(sns,:))); y=avgSynchB(sns,:); yy=y(xx); z=CISynchB(sns,:); zz=z(xx);
smyy = smooth(xx, yy, 10)';
plot(xx, smyy,'color',[194 87 87]./255,'linewidth',2)
fill([xx xx(end:-1:1)], [smyy+zz smyy(end:-1:1)-zz(end:-1:1)], [194 87 87]./255, 'FaceAlpha',.5, 'EdgeColor','none')
axis([0 65 450 750]); 
xlabel('Time (hr)');ylabel('Nuclear Area (pix.)'); 
set(gca,'xtick',0:20:60,'xticklabel',[0:20:60]./4,'fontsize',14,'ytick',500:100:700); box on
line([mean(g1_ctrl)*4 mean(g1_ctrl)*4],[455 930],'linestyle','-','color',[0 83 125]./255,'linewidth',1.5);
line([mean(g1_sns)*4 mean(g1_sns)*4],[455 930],'linestyle','-','color',[194 87 87]./255,'linewidth',1.5);

