%========================
%Data Input
%========================
clear all;
tic
cd('path to where input file can be found');
%Enter your input filename here (including file extension) - 
filename = 'InputFilename.txt';
[TimeData, LabelData, LabelNames] = platereaderread3(filename);
%Platereaderread3 saves the input data into
%InputFilename.mat which is quicker to load that reading the original data
%file.

% load BC-latency-2.mat

toc
disp('data loaded')

%%
%========================
%set smoothing parameters
%========================
%specify a smoothing method for absorbance data and fluorescence data.
%Options are - '','ma',fb','co'.  See platereadersmooth.m for more info.
Amethod = '';  
Fmethod = '';
%========================
%set data parameters
%========================
timepoints = [1:end]; %What timepoints to consider, all by default.
replicates = 3; %the number of replicates of each sample.  For example, 
%if this is 3, then sample one is the first three wells, sample two 
%is the next three wells etc.
Absbackgroundsamples = [4]; %the sample that should be used as an absorbance background
GFPbackgroundsamples = [4]; %the sample that should be used as a fluorescence background
experimentsamples = [7:4:31,8:4:32]; % experimental samples

names = ['sample 1      ';...%name of sample 1, each row needs to be same length
         'sample 2      ';...%name of sample 2
         'sample 3      ';...%name of sample 3
         'GFP background';...
         'Abs background';...
         ];
wellnames = cellstr(names);

samples = [experimentsamples,GFPbackgroundsamples,Absbackgroundsamples];
Absconversion = 'ce'; %specify what conversion to use on the absorbance data
%Options are:
%'ab' for no conversion
%'OD' if a conversion from absorbance to OD exists
%'ce' if a conversion from absorbance to cell density exists
GFPconversion = 'abs'; %specify what conversion to use on the fluorescence data
%Options are:
%'rel' for no conversion
%'abs' if a conversion from relative to absolute units exists.
%========================
%set display parameters
%========================
plotunsmooth = 0; %1 to plot unsmoothed data on the plots
plotsmooth = 1; %1 to plot smooted data on the plots
%Define some figure handles
raw = 1; %define some figure handles
synthesis = 2;

%based on the unit conversion, set the axis labels for later
if Absconversion == 'Ab'
    cellunits = 'Absorbance @ 600nm';
elseif Absconversion == 'OD'
    cellunits = 'Optical Density (600nm)';
elseif Absconversion == 'ce'
    cellunits = 'Cells per well';
end
if GFPconversion == 'rel'
    GFPunits = 'GFP counts (A.U.)';
    GFPperODtitle = 'GFP per OD (counts/OD)';
elseif GFPconversion == 'abs'
    GFPunits = 'GFP (molecules)';
    GFPperODtitle = 'GFP molecules per cell';
    GFPsynthesistotal = 'GFP synthesis rate (apparent molecules/well.s)';
    GFPsynthesispercell = 'GFP synthesis rate (apparent molecules/CFU.s)';
end

%========================
%Reshape data
%========================
%Start the data processing 

Time = TimeData(timepoints);
rawAbs = LabelData(timepoints,:,1);
rawGFP = LabelData(timepoints,:,2);
%reshape time series so that the third dimension represents the different replicates
orderedAbs = reshape(rawAbs,[size(rawAbs,1),...
    replicates,size(rawAbs,2)/replicates]);
orderedGFP = reshape(rawGFP,[size(rawGFP,1),...
    replicates,size(rawGFP,2)/replicates]);
%reorder the dimensions of the array so that the third dimension is the
%number of replicates.
orderedAbs = permute(orderedAbs,[1,3,2]);
orderedGFP = permute(orderedGFP,[1,3,2]);
uAbs = orderedAbs(:,samples,:);
uGFP = orderedGFP(:,samples,:);
toc %check how much time this took
disp('data sorted')
   

%========================
%Time Series Smoothing
%========================
%Put the data to be smoothed into a cell array so the function can handle
%arbitrary numbers of labels.
if ~isempty(Amethod)
    timeseries{1} = uAbs;
    [smoothtimeseries] = platereadersmooth(timeseries,Amethod);
    %now get the smoothed data back out into individual arrays.
    fAbs = smoothtimeseries{1};
else
    fAbs = uAbs;
end
if ~isempty(Fmethod)
    timeseries{1} = uGFP;
    [smoothtimeseries] = platereadersmooth(timeseries,Fmethod);
    %now get the smoothed data back out into individual arrays.
    fGFP = smoothtimeseries{1};
else
    fGFP = uGFP;
end
toc
disp('data smoothed')

%========================
%Background subtraction
%========================
%This implementation assumes the OD of your samples is similar to that of the
%control at each time point.

if ~isempty(Absbackgroundsamples)
    correcteduAbs = Absbackgroundsubtract(uAbs);
    correctedfAbs = Absbackgroundsubtract(fAbs);
else
    correcteduAbs = uAbs;
    correctedfAbs = fAbs;
end

if ~isempty(GFPbackgroundsamples)
    [correctedfGFP]=FPbackgroundsubtract(fGFP);
    [correcteduGFP]=FPbackgroundsubtract(uGFP);
else
    correctedfGFP = fGFP;
    correcteduGFP = uGFP;    
end

%========================
%Unit conversion
%========================
if Absconversion == 'Ab'
    ucells = correcteduAbs;
    fcells = correctedfAbs;
elseif Absconversion == 'OD'
    ucells = 3.113.*correcteduAbs-0.0158;
    fcells = 3.113.*correctedfAbs-0.0158;
elseif Absconversion == 'ce'
    uOD = 3.113.*correcteduAbs-0.0158;%convert to OD
    fOD = 3.113.*correctedfAbs-0.0158;
    ucells = uOD*5E8;%convert to cells/ml
    fcells = fOD*5E8;
    ucells = ucells*0.2;%convert to cells per well.
    fcells = fcells*0.2;
end

if GFPconversion == 'rel'
    
elseif GFPconversion == 'abs'
    correctedfGFP=1.16E-6*(correctedfGFP)+9.95E-4;%convert to nmoles
    correctedfGFP=correctedfGFP*6.022E14;%convert to molecules
    correcteduGFP=1.16E-6*(correcteduGFP)+9.95E-4;%convert to nmoles
    correcteduGFP=correcteduGFP*6.022E14;%convert to molecules
end
    
%===========================

%Calculate the dilution rate, assuming exponential growth.  The function
%returns the means and the confidence intervals
[fgammad, Mfgammad, CIfgammad, Mfmu, CIfmu] = AbsProcess(Time,fcells);
[ugammad, Mugammad, CIugammad, Mumu, CIumu] = AbsProcess(Time,ucells);

%===========================
%GFP Synthesis Rate Calculations

[MdfGFPdt, CIdfGFPdt, Mfsynthesis, CIfsynthesis, Mfaccumulation, CIfaccumulation] = FPBulkProcess(Time,fcells,correctedfGFP);
[MduGFPdt, CIduGFPdt, Musynthesis, CIusynthesis, Muaccumulation, CIuaccumulation] = FPBulkProcess(Time,ucells,correcteduGFP);

%Get stats for each of the raw/smooth time series replicates
[Mucells,CIucells] = platereaderstats(ucells);
[Mfcells,CIfcells] = platereaderstats(fcells);

[MuGFP,CIuGFP] = platereaderstats(uGFP);
[MfGFP,CIfGFP] = platereaderstats(fGFP);
toc
disp('data processed')

%========================
%Data plotting
%========================

%how many columns should I plot?
figure1columns = size(experimentsamples+1,2);
figure2columns = size(experimentsamples,2);

figure(raw);

subplot(3,1,1)
    hold off;
    if plotunsmooth
        semilogy(Time,Mucells(:,1:figure1columns),'s','MarkerSize',5)
        hold on;
     
    end
    
    if plotsmooth
      semilogy(Time,Mfcells(:,1:figure1columns),'LineWidth',1.5); 
    end
    
    title('Growth Curve')
    xlim([0,inf]);
subplot(3,1,2)
    hold off;
    if plotunsmooth
        errorbar(repmat(Time,[1,figure1columns]),MuGFP(:,1:figure1columns),CIuGFP(:,1:figure1columns),'-.k')
        hold on;
    end
    
    if plotsmooth
        errorbar(repmat(Time,[1,figure1columns]),MfGFP(:,1:figure1columns),CIfGFP(:,1:figure1columns))
    end
    
    title('Fluorescence Signal')
    ylabel(GFPunits);
    xlim([0,inf]);
subplot(3,1,3)
    hold off; 
    if plotunsmooth
        errorbar(repmat(Time(1:end-1),[1,figure2columns]),MduGFPdt(:,1:figure2columns),CIduGFPdt(:,1:figure2columns),'o','MarkerSize',5)
        hold on;
%         ylim([-200, 500]);
    end
    if plotsmooth
            errorbar(repmat(Time(1:end-1),[1,figure2columns]),MdfGFPdt(:,1:figure2columns),CIdfGFPdt(:,1:figure2columns),'LineWidth',0.5)
    end
    title('Bulk GFP synthesis rate');
    xlabel('Time (s)');
    ylabel(GFPsynthesistotal);
    xlim([0,inf]);

figure(synthesis)
    hold off;
    if plotunsmooth
        errorbar(repmat(Time(1:end-1),[1,figure2columns]),Musynthesis(:,1:figure2columns),CIusynthesis(:,1:figure2columns),'-k')
        hold on;
        
    end
    if plotsmooth
        errorbar(repmat(Time(1:end-1),[1,figure2columns]),Mfsynthesis(:,1:figure2columns),CIfsynthesis(:,1:figure2columns))
        hold on;
    end
    title('GFP synthesis rate per cell','FontSize',18,'FontName','Helvetica');
%     ylim([-200, 700]);
    xlim([0,inf]);
    set(gca,'FontSize',14,'FontName','Helvetica');
    xlabel('Time (s)','FontSize',18,'FontName','Helvetica');
    ylabel(GFPsynthesispercell,'FontSize',18,'FontName','Helvetica')
    legend(names(experimentsamples,:));