
% This code use the array 'trackmate', once cell segmentation is performed, to extract the cell cycle trajectories
% and determine the different parameters of interest, e.g the cell cycle
% duration ,the initial cell dry mass and the final cell drymass.
% This code is in two parts, the first part deals with the detection of the dividing cells, and the second part with
% the analysis of the cell cycle trajectories.

%%%%%%%%%%%%%%%% Detection of divison looking for new cell appearing close to a max. thickness signal

% This first part of the code use the dataset 'trackmate' to detect cell divisions
% The output is the array 'divisions' featuring the positions of the
% divisions in the 'trackmate' dataset, the generation number, and the
% track number of the possible daughter cell. In the 'trackmate array',
% colum 2 'LABEL' becomes the 'generation number' from 0 to inf



bVisu=false ;  % boolean for visual inspection

detection_threshold_thickness=5; %thickness threshold of mitosis detection [µm]


detection_threshold_time1=200;   %time threshold in minute, maximum mitotic duration
detection_threshold_time2=60;    %time threshold in minute, duration of mother track after division detection
detection_threshold_distance=50;   %distance max. of the apparition of the 'new' daughter cell in µm
plot_length=64;plot_length_min=-2;

trackmate(:,2)=0;
track_pos=trackmate(:,4);track_pos=unique(track_pos);% track number
start_pos=find((trackmate(2:end,4)-trackmate(1:end-1,4))>0)+1; % line number of starting points of the tracks
if (not(isempty(start_pos)))
    start_trackmate=trackmate(start_pos,:);
end
divisions=struct('N',[],'daughters',[],'pos',0,'posini',0,'frame',0,'generation',0,'detection',0,'detection_pos',0,'detection_thickness',0);
divisions_count=0;
t_divisions_count=0;
for k=track_pos'
    if (mod(k,1000)==0)
        display(strcat(num2str(k),'\',num2str(size(track_pos,1))));
    end
    divisions_count=divisions_count+1;
    pos=find(trackmate(:,4)==k); % analysis of track number k
    pos_initial=pos(1);
    track_thickness=trackmate(pos,15);
    track_phase=trackmate(pos,8);
    track_CDM=trackmate(pos,12);
    track_x=trackmate(pos,6);
    track_y=trackmate(pos,7);
    track_frame=trackmate(pos,10);
    track_daughter=(trackmate(pos(1),2)>0);  % tag on track corresponding to a daughter cell already found
    track_CDM_diff=diff(track_CDM);
    pos=find(track_CDM_diff>500);
    if ~isempty(pos)
        if (pos>1) & (pos<length(track_thickness)-2)
            track_CDM(pos+1)=(track_CDM(pos+2)+track_CDM(pos))/2;
            track_thickness(pos+1)=(track_thickness(pos+2)+track_thickness(pos))/2;
        end
    end
    window=[ones(60/rate,1)' -1*ones(60/rate,1)'];
    track_CDM_filter=-conv(track_CDM,window')/length(window);
    track_CDM_filter=track_CDM_filter(1:length(track_CDM));
    track_CDM_filter=track_CDM_filter*2./circshift(track_CDM,[60/rate 0])*1000;
    track_thickness_filter=-conv(track_thickness,window')/length(window);
    track_thickness_filter=track_thickness_filter(1:length(track_CDM))*10;

    % division detection
    results=double(zeros(length(track_frame),6));
    for m=1:length(track_frame)
        % calculation of distances in time and length from track(m) to starting pooints
        % of all tracks
        d=sqrt((trackmate(start_pos,6)-track_x(m)).^2+(trackmate(start_pos,7)-track_y(m)).^2);  %distance in µm
        t1=trackmate(start_pos,10)-track_frame(m);  %time in numbers of frame
        t2=track_frame(end)-trackmate(start_pos,10);  %time in numbers of frame
        r=[start_trackmate(:,4) d t1 t2]; % for every point of the trajectories, calculation of the distance to all trajectories starting points
        % r(:,1)=N°track, r(:,2)=distance r(:,3)=delta frame r(:,4)=delta frame end
        % first criteria of division detection 
        pos=find((r(:,2)<detection_threshold_distance)&(r(:,3)*rate<detection_threshold_time1)&(r(:,3)>0)&(r(:,4)*rate>detection_threshold_time2));
        results(m,2)=track_thickness(m);
        
        % if a new daughter cell is present in the close surrounding
        if not(isempty(pos))
            if (size(pos,1)>1)
                [mmin,ppos]=min(r(pos,2)+r(pos,3)*rate*1);
                pos=pos(ppos); %we only look at the closet new daughter cell
            end
            results(m,1)=size(pos,1); % number of possible daughter tracks detections
            % results(m,2)=track_thickness(m); 
            results(m,3)=r(pos(1),1); % number of possible daughter tracks 
            results(m,4)=r(pos(1),2); % distance
            results(m,5)=r(pos(1),3); % delta frame
            results(m,6)=(results(m,3)>0)&(results(m,2)>detection_threshold_thickness);     %boolean of a possible division detection 
            results(m,7)=size(pos,1);
        end
    end
    % first criteria of division detection 
    results_division=results(:,3).*results(:,6);  % list of possible daughter cells, when the cell thickness is above the threshold in the cell track of interest
    results_division=unique(results_division);results_division(results_division==0)=[];
    if isempty(results_division)
        trackmate(pos_initial:pos_initial+length(track_frame)-1,2)=trackmate(pos_initial,2);
    else
    % if divisions detections are found, the exact positions in the tracks are extracted    
        for m=1:size(results_division,1)
            pos=find(results(:,3)==results_division(m,1)); % for every possible daughter cells 
            [mmax,ppos]=max(results(pos,2));               % we look at the maxium thickness of the cell track 
            results_division(m,2)=pos(ppos);               % position in the cell track of interest
            results_division(m,3)=mmax;                    % maximum thickness found
        end
        divisions(k).N=k;
        divisions(k).detection=results_division(:,1)';
        divisions(k).detection_pos=results_division(:,2)';
        divisions(k).detection_thickness=results_division(:,3)';
        % remove division before first 6 hours when it is a daughter cell
        pos=find(results_division(:,2)<6*60/rate);
        if not(isempty(pos))&(track_daughter)
            results_division(pos,:)=[];
        end
        % remove division too close by 6 hours
        if (size(results_division,1)>1)
            flag=true;
            while(flag)
                pos=find(abs((results_division(2:end,2)-results_division(1:end-1,2)))<360/rate);
                if not(isempty(pos))
                    if (results_division(pos+1,3)>results_division(pos,3))
                        results_division(pos(1),:)=[];
                    else
                        results_division(pos(1)+1,:)=[];
                    end
                else
                    flag=false;
                end
                if not(size(results_division,1)>1)
                    flag=false;
                end
            end
        end
        if (isempty(results_division))
            % generation calculation
            trackmate(pos_initial:pos_initial+length(track_frame)-1,2)=trackmate(pos_initial,2);
        else
            t_divisions_count=t_divisions_count+size(results_division,1);
            for m=1:size(results_division,1)
                pos=find(trackmate(start_pos,4)==results_division(m));  
                pos=start_pos(pos);
                % generation calculation
                if (m==1)
                    trackmate(pos_initial:pos_initial+results_division(1,2)-1,2)=trackmate(pos_initial,2);
                else
                    trackmate(pos_initial+results_division(m-1,2):pos_initial+results_division(m,2)-1,2)=trackmate(pos_initial+results_division(m-1,2)-1,2)+1;
                end
                % the generation number is given as the one of the mother cell +1 
                trackmate(pos,2)=trackmate(pos_initial+results_division(m,2)-1,2)+1;
            end
            trackmate(pos_initial+results_division(end,2):pos_initial+length(track_frame)-1,2)=trackmate(pos_initial+results_division(end,2)-1,2)+1;
            divisions(k).daughters=results_division(:,1)';
            divisions(k).generation=trackmate(pos_initial,2);
            divisions(k).pos=results_division(:,2)'+pos_initial-1;
            divisions(k).posini=pos_initial;
            divisions(k).frame=trackmate(pos_initial,10)+results_division(:,2)';
            if(divisions(k).generation>0) % if it is daughter cell, the beginning of the track is given as division
                divisions(k).pos=[pos_initial divisions(k).pos];
            end
        end
    end
    if (bVisu)&(~isempty(results_division))
        figure(1);clf;subplot(2,2,1);plot(results(:,2),'Color','k');ylim([0 10]);title('thickness');
        subplot(2,2,2);plot(results(:,4),'Color','k');ylim([0 30]);title('neighbour distance');
        subplot(2,2,3);plot(results(:,3),'Color','k');ylim([0 100]);title('neighbour track number');
        subplot(2,2,4);plot(results(:,6),'Color','k');ylim([0 2]);title('boolean');
        text(5,1.5,strcat('k=',num2str(k)));
        text(5,1.2,strcat(num2str(results_division)));hold off;drawnow;
        set(1,'units','normalized','position',[0.01 0.01 0.45 0.8]);set(1,'Color','w');drawnow;
        
        figure(4);clf;
        set(4,'units','normalized','position',[0.47 0.01 0.35 0.9]);set(4,'Color','w');
        subplot(4,1,1); plot(track_CDM,'Color','k');ylim([0 1250]);xlim([0 800]);hold on;
        set(gca,'XTick',[0,72,144,216,288,360,432,504,576]);
        if ~isempty(results_division)
            pos=divisions(k).pos-divisions(k).posini;
            track_abcisse=[1:length(track_CDM)];
            track_red=track_abcisse*0;track_red(max(pos,1))=2000;
            plot(track_red,'Color','r');
        end
        xl=get(gca, 'XTickLabel');xl=num2str(0.1*round((str2double(xl(:,:))*rate/60)*10));set(gca, 'XTickLabel',xl);set(gca,'FontSize',8);
        xlabel('Duration [hours]','FontSize',8); ylabel('Cell dry mass [pg]','FontSize',8);
        subplot(4,1,2);plot(track_thickness,'Color','k');ylim([0 10]);xlim([0 800]);hold on;
         if ~isempty(results_division)
              plot(track_red,'Color','r');
         end
        set(gca,'XTick',[0,72,144,216,288,360,432,504,576]);
        xl=get(gca, 'XTickLabel');xl=num2str(0.1*round((str2double(xl(:,:))*rate/60)*10));set(gca, 'XTickLabel',xl);set(gca,'FontSize',8);
        xlabel('Duration [hours]','FontSize',8); ylabel('Thickness [µm]','FontSize',8);
        subplot(4,1,3); plot(track_CDM,'Color','k');ylim([0 1250]);xlim([0 1000]);hold on;
        plot(track_CDM_filter,'Color','r');xlim([0 800]);hold on;ylim([0 1250])
        set(gca,'FontSize',8);
       
    end
end
divisions(max(trackmate(:,4))).N=max(trackmate(:,4));
display(' ');
display(strcat('divisions_count=',num2str( t_divisions_count )));

% Extraction of the cell cycle trajectories
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This second part of the code use the array 'divisions' and the dataset 'trackmate' to
% extract the cell cycle trajectories. For each trajectory a linear fit is
% performed. Based on the R2 and the slope of the cell mass trajectory a
% criteria is applied to remove 'bad' cell cycle tracks.
% The ouput is the structure lineage, the arrays 'intercept' (initial dry
% mass (pg]), 'last' (en dry mass [pg]) and 'slope' (growth rate [pg/h])
%
% The logfit function can be found in the Mathworks file exchange.

lineage=[];
lineage=struct('N0',0,'N',0,'N1',0,'N2',0,'type',0,'xi',0,'yi',0,'ti',0,'xf',0,'yf',0,'tf',0,'D',0,'Vi',0,'Vf',0,'Phase',0,'RawPhase',0,'VolumePhase',0,'Track',0);
lineage_count=0;
divisions_count=0;
for k=1:size(divisions,2)
    if (mod(k,1000)==0)
        display(strcat(num2str(k),'\',num2str(size(divisions,2))));
    end
    divisions_count=divisions_count+1;
    if not(isempty(divisions(divisions_count).N))

        if (size(divisions(divisions_count).pos,2)>1)
            pos=divisions(divisions_count).pos;
            for n=1:size(divisions(divisions_count).pos,2)-1
                track=trackmate(pos(n):pos(n+1),:);
                if(size(track,1)>0)
                    lineage_count=lineage_count+1;
                    lineage(lineage_count).N0=divisions(divisions_count).generation;
                    lineage(lineage_count).N=divisions(divisions_count).N;
                    lineage(lineage_count).N1=0;
                    lineage(lineage_count).N2=0;
                    lineage(lineage_count).type=1;
                    lineage(lineage_count).x=track(:,6);
                    lineage(lineage_count).y=track(:,7);
                    lineage(lineage_count).xi=track(1,6);
                    lineage(lineage_count).yi=track(1,7);
                    lineage(lineage_count).xf=track(end,6);
                    lineage(lineage_count).yf=track(end,7);
                    lineage(lineage_count).ti=track(1,10)*rate/60;
                    lineage(lineage_count).tf=track(end,10)*rate/60;
                    lineage(lineage_count).D=size(track,1)*rate/60;
                    lineage(lineage_count).Phase=track(:,8);
                    lineage(lineage_count).Area=track(:,11);
                    lineage(lineage_count).DryMassPhase=track(:,12);
                    lineage(lineage_count).VolumePhase=track(:,13);
                end
            end
      
        end
    end
end

if (size(lineage,2)>1)
    tobesaved=zeros(1,1);
    n_lineage=1;
    for k=1:size(lineage,2)
        if (length(lineage(k).DryMassPhase)*rate/60>4)
            tobesaved(1:length(lineage(k).DryMassPhase),n_lineage)=lineage(k).DryMassPhase;
            n_lineage=n_lineage+1;
        end
    end
    n_lineage=n_lineage-1;
    if (n_lineage>0)
        b=[0:rate:(size(tobesaved,1)-1)*rate];b=b';
        tobesaved=[b tobesaved];
        xlswrite(strcat(folder_in,'drymass.xlsx'),tobesaved);
        display(strcat(folder_in,'drymass.xlsx    N=',num2str(n_lineage)))
    end
end

filename=strcat(folder_in,strcat('lineage.mat'));
save(filename,'lineage','-v7.3');

%%%%%%%%%%%% Cell growth_rate fit analysis %%%%%%%%%%%%%%%%%%%%

bVisu=false;
if (size(lineage,2)>1)
    display('Cell cycles growth_rate fit analysis ...');
    %%% Analyse growth_rate
    time=transpose([1:rate:10000]);
    count=0;t_initial=zeros(length(lineage),1);cycle_length=zeros(length(lineage),1);distance=zeros(length(lineage),1);
    int_phase=zeros(length(lineage),1);R=zeros(length(lineage),1);slope=zeros(length(lineage),1);intercept=zeros(length(lineage),1);last=zeros(length(lineage),1);ratio=zeros(length(lineage),1);
    for k=1:length(lineage)
        if (mod(k,500)==0)
            display(strcat(num2str(k),'\',num2str(size(lineage,2))));
        end
        t_initial(k,1)=lineage(k).ti;
        if (lineage(k).type>0)
            cycle_length(k,1)=lineage(k).tf-lineage(k).ti;
        end
        distance(k,1)=sqrt((lineage(k).xf-lineage(k).xi)^2+(lineage(k).yf-lineage(k).yi)^2);
        int_phase(1:size(lineage(k).DryMassPhase,1),k+1)=transpose(lineage(k).DryMassPhase);
        % slope analyse
        if(length(lineage(k).DryMassPhase(:,1))>12)
            [slope(k,1), intercept(k,1), R(k,1)] = logfit(time(6:size(lineage(k).DryMassPhase,1)-5),lineage(k).DryMassPhase(6:end-5,1),'linear'); close
            last(k,1) = intercept(k,1)+size(time(6:size(lineage(k).DryMassPhase,1)-5),1)*slope(k,1)*rate;
            slope(k,1)=slope(k,1)*60;
            ratio(k,1)=last(k,1)./intercept(k,1);
        else
            slope(k,1)=0; intercept(k,1)=0; R(k,1)=0;
        end
        if (bVisu)
            fig1=figure(1);
            clf;plot(time(1:size(lineage(k).DryMassPhase,1)-5)/60,lineage(k).DryMassPhase(1:end-5,1));hold on;ylim([0 1000]);
            plot(time(1:size(lineage(k).DryMassPhase,1))/60,intercept(k,1)+time(1:size(lineage(k).DryMassPhase,1))*slope(k,1)/60);
            text(1.5,850,strcat('track=',num2str(round(lineage(k).N))));
            text(1.5,800,strcat('R^2=',num2str(round(R(k,1)))));
            text(1.5,750,strcat('growth rate=',num2str(  round(10*slope(k,1))*0.1 ),' pg/h'));
            text(1.5,700,strcat('division ratio=',num2str(  round(10*ratio(k,1))*0.1  )));
            text(1.5,650,strcat('cycle length=',num2str(  round(10*cycle_length(k,1))*0.1 ),' hours'));
            xlabel('Duration [hours]','FontSize',12);ylabel('Cell dry mass [pg]','FontSize',12);
            set(1,'Color','w');drawnow;pause(0.5);
        end
    end

    % The cell cycle trajectories are filtered out 
    selection=(cycle_length(:,1)>5)& (slope(:,1)<25) & (R(:,1)<10000);
    selection_count=sum(selection);
    display(strcat('selection_count=',num2str( sum(selection) )));
    display(strcat('overall_efficiency=',num2str(0.1*round(selection_count/size(lineage,2)*1000)),'%'))
    display(strcat('Cell cycle duration=',num2str(0.1*round(mean(cycle_length)*10)),'H+/-',num2str(0.1*round(std(cycle_length)*10))));
    ratio(selection==0)=[];last(selection==0)=[];intercept(selection==0)=[];slope(selection==0)=[];cycle_length(selection==0)=[];t_initial(selection==0)=[];R(selection==0)=[];
    lineage(selection==0)=[];
end

%%%%%%%%%%%% Cell cycle trajectories measurements %%%%%%%%%%%%%%%%%%%%

if (size(lineage,2)>1)
    for k=1:length(lineage)
        CDM=lineage(k).DryMassPhase;
        CDM2=CDM;
        pos=find( (diff(CDM)/rate*10>25) | (diff(CDM)/rate*10<-25)) ;
        CDM2(pos)=[];
        
        if (length(CDM2)>30)
            s_CDM2=CDM2(1:15);s_CDM2=sort(s_CDM2);
            f_CDM2=CDM2(end-15:end);f_CDM2=sort(f_CDM2,'descend');
            cycle_length(k,1)=lineage(k).tf-lineage(k).ti;  % cell cycle duration of track number k
            intercept(k)=mean(s_CDM2(1:uint8(50/rate)));    % initial cell dry mass of track number k
            last(k)=mean(f_CDM2(1:uint8(50/rate)));         % final cell dry mass of track number k
            ratio(k)=last(k)/intercept(k);                  % ratio between final and initial cell dry mass of track number k
            slope(k)=(last(k)-intercept(k))/cycle_length(k);
        end
    end
end


