%% generation times
out = out_filtered;
treatment = 82;
frameTime = 3;
totalDivisions = zeros(1000,1);
totalDivDurations = zeros(1000,1);
totalcells = zeros(1000,1)
for jj = 1:length(out)
    for kk = 1:length(out(jj).Lines)
        indices = find(out(jj).Lines(kk).divDurations>7 & out(jj).Lines(kk).divDurations<30); %remove segmentation-division error outliers
        if numel(indices)>0
            frames = out(jj).Lines(kk).divTimes(indices);
            indices(frames==0) = [];
            frames(frames==0) = [];
            totalDivisions(frames) = totalDivisions(frames) + 1;
            totalDivDurations(frames) = totalDivDurations(frames) + out(jj).Lines(kk).divDurations(indices)';
            totalcells(frames) = totalcells(frames) + 1;
        end
    end
end
framesAnalysis = 1:400;
t = (framesAnalysis - treatment)*frameTime;
y = totalDivDurations(framesAnalysis)./totalDivisions(framesAnalysis)*frameTime;
figure;
plot(t,y,'g')


%% mismatch rate curve
out= out_filtered
treatment = 66;
frameTime = 3;
totalCells = zeros(1000,1);
totalFoci = zeros(1000,1);
nFoci = nan(10000,1);
ll = 1;
for jj = 1:length(out)
    for kk = 1:length(out(jj).Lines)
        indices = find(out(jj).Lines(kk).YFPAvg>0);
        frames = out(jj).Lines(kk).frameNo(indices);
        totalCells(frames) = totalCells(frames) + 1;
	%remove single off frames when foci are visible for multiple frames:
	filterFociNo = out(jj).Lines(kk).fociNo;
	for ii = 1:length(out(jj).Lines(kk).fociNo)-2
    	if out(jj).Lines(kk).fociNo(ii) > 0 && out(jj).Lines(kk).fociNo(ii+1) == 0 && out(jj).Lines(kk).fociNo(ii+2) > 0
        	filterFociNo(ii+1) = out(jj).Lines(kk).fociNo(ii);
    	end
        end
        fociIndices = find(diff(filterFociNo)>0)+1;
        if ~isempty(fociIndices)
            totalFoci(out(jj).Lines(kk).frameNo(fociIndices)) = totalFoci(out(jj).Lines(kk).frameNo(fociIndices)) + 1;
        end
    end
end
framesAnalysis =1:1000;
t = (framesAnalysis-treatment)*frameTime;
nFoci(isnan(nFoci)) = [];
y = (totalFoci(framesAnalysis)./totalCells(framesAnalysis))/frameTime;
figure;
plot(t,y,'c')


%% CFP average response curves
out = out_filtered;
treatment = 82;
frameTime=3;
totalCells = zeros(1000,1);
totalCFP = zeros(1000,1);
for jj = 1:length(out)
    for kk = 1:length(out(jj).Lines)
        indices = find(out(jj).Lines(kk).CFPAvg>0);
        frames = out(jj).Lines(kk).frameNo(indices);
        totalCells(frames) = totalCells(frames) + 1;
        totalCFP(frames) = totalCFP(frames) + out(jj).Lines(kk).CFPAvg(indices)';
    end
end
framesAnalysis = 1:330;
%correcting for fluroescent protein maturation:
maturationTime = 6.4;
lambda=log(2)/(maturationTime/frameTime);
maturation = lambda*exp(-lambda*[0:19]);
maturation = maturation/sum(maturation);
t = (framesAnalysis-treatment)*frameTime;
y = totalCFP(framesAnalysis)./totalCells(framesAnalysis);
deconvolutedCFP = deconv(y(framesAnalysis),maturation);
figure;
plot(t(1:end-19),deconvolutedCFP);


%% survival time distribution
out = out_filtered;
treatment = 0;
frameTime=3;
survival = [];
totalcells = zeros(1000,1);
for jj = 1:length(out)
    for kk = 1:length(out(jj).Lines)
        indices = find(out(jj).Lines(kk).mCherryAvg>0);
        frames = out(jj).Lines(kk).frameNo(indices);
        if max(frames) > treatment
            survival = [survival, max(frames) - treatment];
        end
        totalcells = totalcells + 1;
    end
end
[y,x,flo,fup] = ecdf(survival);
t = x*frameTime;
figure;
stairs(t,1-y,'c')

%% elongation rate 
out = out_filtered;
treatment = 82;
elongationRateMean = [];
framesAnalysis = 1:1000;
pixel = 0.13; %pixel length in micrometer
frameTime = 3;
t = (framesAnalysis-treatment)*frameTime;
elongationRateSum = zeros(1,length(t));
cellCount = zeros(1,length(t));

for jj = 1:length(out)
    for kk = 1:length(out(jj).Lines)
        frames = out(jj).Lines(kk).frameNo;
        for ii = 2:numel(frames)-1
            cellDiv = find(out(jj).Lines(kk).divTimes == frames(ii));
            cellDiv2 = find(out(jj).Lines(kk).divTimes == frames(ii-1));
            if isempty(cellDiv) && isempty(cellDiv2) && frames(ii)>0
                elongationRate = (((out(jj).Lines(kk).MajorAxisLength(1,ii) -...
                    out(jj).Lines(kk).MajorAxisLength(1,ii-1))*pixel)/frameTime);
                elongationRateSum(frames(ii)) = (elongationRateSum(frames(ii)) + elongationRate);
                cellCount(frames(ii)) = cellCount(frames(ii)) + 1;
            end
        end
    end
end

figure;
plot(t,elongationRateSum./cellCount)

%% CFP expression rate plot
out = out_filtered;
frameTime = 3;
treatment = 52;
framesAnalysis = 1:500;
cellCount = zeros(1,length(framesAnalysis));
CFPSum = zeros(1,length(framesAnalysis));
dCFPdtASum = zeros(1,length(framesAnalysis));
dCFPdtACellCount = zeros(1,length(framesAnalysis));
expressionRate = [];
for jj = 1:length(out)
    for kk = 1:length(out(jj).Lines)
        indices = find(out(jj).Lines(kk).frameNo>0);
        frames = out(jj).Lines(kk).frameNo(indices);
        if frames(end)> 70
        cellCount(frames) = cellCount(frames) + 1;
        CFPSum(frames) = CFPSum(frames) + out(jj).Lines(kk).CFPAvg(indices);    
        for ii = 2:numel(indices)-1 
            cellDiv = find(out(jj).Lines(kk).divTimes == frames(ii), 1);
            cellDiv2 = find(out(jj).Lines(kk).divTimes == frames(ii-1),1);
            if isempty(cellDiv) && isempty(cellDiv2) && frames(ii)>0
            dCFPdtA = (((out(jj).Lines(kk).CFPTot(1,indices(ii)) -...
                    out(jj).Lines(kk).CFPTot(1,indices(ii)-1)))/frameTime)/...
                    out(jj).Lines(kk).Area(1,indices(ii));
            dCFPdtASum(frames(ii)) = (dCFPdtASum(frames(ii)) + dCFPdtA);
            dCFPdtACellCount(frames(ii)) = dCFPdtACellCount(frames(ii)) + 1;  
            end
        end
        end
    end
end

CFPMean = CFPSum./cellCount;
dCFPdtAMean = dCFPdtASum./dCFPdtACellCount;
smoothingWindow = 3;
dCFPdtAMovMean = movmean(dCFPdtAMean,smoothingWindow);
expressionRate = dCFPdtAMovMean;

%correcting for fluroescent protein maturation:
maturationTime = 6.4;
lambda=log(2)/(maturationTime/frameTime);
maturation = lambda*exp(-lambda*[0:19]);
maturation = maturation/sum(maturation);

deconvolutedRate = deconv(expressionRate(framesAnalysis),maturation);
deconvolutedMean = deconv(CFPMean(framesAnalysis),maturation);
t = (framesAnalysis-treatment)*frameTime;

figure
yyaxis left
y1 = deconvolutedRate;
plot(t(1:end-19),deconvolutedRate)
hold on
yyaxis right
y2 = deconvolutedMean;
plot(t(1:end-19),y2)
hold off

%% cell length
out = out_filtered;
pixel = 0.13;
treatment = 82;
frameTime=3;
totalCells = zeros(1000,1);
totalCFP = zeros(1000,1);
totalLength = zeros(1000,1);
ll = 1;
for jj = 1:length(out)
    for kk = 1:length(out(jj).Lines)
        indices = find(out(jj).Lines(kk).CFPAvg>0);
        frames = out(jj).Lines(kk).frameNo(indices);
        totalCells(frames) = totalCells(frames) + 1;
        totalLength(frames) = totalLength(frames) + out(jj).Lines(kk).MajorAxisLength(indices)'*pixel;
    end
end
framesAnalysis = 1:500;
x = (framesAnalysis-treatment)*frameTime;
figure;
hold all
plot(x,totalLength(framesAnalysis)./totalCells(framesAnalysis),'c')
title('Cell length')
hold off

%% diffusion coefficient histogram

[filename,pathname] = uigetfile('*.tracks', 'Tracks data:')
tracks = importdata([pathname filename]);
pixel = 0.096; % length per pixel
dT = 0.01548; % time per frame
DhistMinSteps = 4; % minimum number of steps for a track to be analyzed
rangeD = -0.2:0.1:8 ; % D range for the histogram
nMolecules = max(tracks(:,4));
MSD = zeros(nMolecules,1);
kk = 1;

for ii = 1:nMolecules
    
    xx = find(tracks(:,4)==ii);
    
    if numel(xx) > DhistMinSteps 
        % sum all squared displacement in the track
        for jj = 1:numel(xx)-1      
            MSD(kk) = MSD(kk) + ((tracks(xx(jj+1),1) - tracks(xx(jj),1))^2 +...
                (tracks(xx(jj+1),2) - tracks(xx(jj),2))^2);
        end
        MSD(kk) = MSD(kk)/jj; % mean square displacement
        kk = kk + 1;
    end
    ii
end

MSD(kk:end) = []; % delete unused rows
MSD = MSD * pixel^2; % convert from pixel to length units
% calculate D from MSD and correct for localization noise
D = MSD/(4*dT);

binSpacing = rangeD(2)-rangeD(1);
histDCount = histc(D,rangeD + binSpacing/2);
figure;
bar(rangeD+binSpacing,histDCount./(sum(histDCount)*binSpacing));
xlim([min(rangeD), max(rangeD)]);
xlabel('diffusion coefficient [um^2/s]')
ylabel('probability density')