clear all 
close all
clc

%% SET PARAMETERS
tf = input('Build a new standard curve. Enter the number of training files? ');
ts = input('Number of samples in training files ');
tr = input('Number of replicates per sample in each experiment? ');
f = zeros((ts-1),1);
f(1)=input('What is the lowest template concentration (in log10 copies)? ');
step = input('How many orders of magnitude are there between each standard? '); 
HThresh = input('Slope threshold? ');  % HThresh = slope threshold
z = input('Number (z) standad deviations above background for positive threshold? ');          
JPGcheck = input('Is datatype from 1 = thermal cycler, 2 = microscope TIF, 3 = microscope JPG? ');


if JPGcheck == 1
    FTchange = input('input new IPC threshold? y/n? ', 's' );
    if FTchange == 'y'
        FAMthresh = input('set IPC FAMthreshold: '); 
    else
        FAMthresh = 327.0745;
    end
    
elseif JPGcheck == 2
   FTchange = input('input new IPC threshold? y/n? ', 's' );
    if FTchange == 'y'
        FAMthresh = input('set IPC FAMthreshold: '); 
    else
        FAMthresh = 327.0745;
    end
elseif JPGcheck == 3
   FTchange = input('input new IPC threshold? y/n? ', 's' );
    if FTchange == 'y'
        FAMthresh = input('set IPC FAMthreshold: '); 
    else
        FAMthresh = 9;
    end
end
    
% READ IN TRAINING FILES

HEX = zeros(60,ts*tr*tf);
FAM = zeros(60,ts*tr*tf);

for i =1:tf
   [tfilename, tpathname] = uigetfile('*.xlsx', 'Choose training file to build standard curve');
    cd(tpathname);
    HEX(1:60,(ts*tr*(i-1)+1):(ts*tr*i))= [xlsread(tfilename, 'HEX', 'c2:af61')];
    FAM(1:60,(ts*tr*(i-1)+1):(ts*tr*i)) = [xlsread(tfilename, 'FAM', 'c2:af61')];
end

% CHECK HEX NTCs TO SET THRESHOLD FOR SUCCESSFULLY RUN SAMPLES
HEXtest = max(HEX)-min(HEX);
NTCmean = mean(HEXtest(1:ts:end));
NTCstd = std(HEXtest(1:ts:end));
posTHRESH = NTCmean + z*NTCstd;

%% TRAINING SECTION

% check for failed samples (no HEX or FAM amplification
FAMtest = max(FAM)-min(FAM);
HEXclean = zeros(60, ts*tr*tf);

HEXfail = zeros(1, ts*tr*tf);  %TESTCHECK tells indicies of negatives and invalid
for i=1:ts*tr*tf
    if HEXtest(i)>posTHRESH;
       HEXfail(i)=HEXtest(i);
       HEXclean(:,i)=HEX(:,i);
    else HEXfail(:,i)=NaN;
        HEXclean(:,i)=NaN;
    end
end

FAMfail = zeros(1, ts*tr*tf);
for i=1:ts*tr*tf
    if FAMtest(i)>JPGcheck;      %!!!!!NEED TO SET THIS VALUE APPROPRIATELY
       FAMfail(i)=FAMtest(i);
    else FAMfail(:,i)=NaN;       
    end
end

BOTHfail = isnan(HEXfail).*isnan(FAMfail);
HEXpos = ~isnan(HEXfail);
FAMpos = ~isnan(FAMfail);
valid = ~BOTHfail;

% CHANGE SAMPLES THAT DIDN'T REACH THRESHOLD TO NANs IN HEXclean
HEXclean = zeros(60, ts*tr*tf);
for i=1:ts*tr*tf    
    if HEXtest(i)>posTHRESH;
       HEXclean(:,i)=HEX(:,i);
    else HEXclean(:,i)=NaN;   
    end
end

% smooth HEXclean
HEXSMOOTH = {zeros(60,ts*tr*tf)};
for i=1:ts*tr*tf                        
    HEXCOL = HEXclean(:,i);
    HEXSMOOTH{i} = smooth(HEXCOL,5);
end
HEXSMOOTH = cell2mat(HEXSMOOTH);

% claculates rough first deriv
HDIFF = {zeros(59,ts*tr*tf)};
for i=1:ts*tr*tf                        
    HEXDIFFCOL = HEXSMOOTH(:,i);
    HDIFF{i} = diff(HEXDIFFCOL);
end
HDIFF = cell2mat(HDIFF);
% cleans matrix back to 60 rows
firstrow = zeros(1, ts*tr*tf);          
HDIFF = [firstrow; HDIFF];

% find threshold cycle
HTHRESHcycle = {zeros(1,ts*tr*tf)};
for i = 1:ts*tr*tf
    HTHRESHcycleCOL = HDIFF(:,i);
    HTHRESHcycle{i} = find(HTHRESHcycleCOL > HThresh, 1, 'first');
end
HTHRESHcycleCLEAN = cellfun('isempty', HTHRESHcycle);
% for unamplified samples, HTHRESHcycle turned to NaN.
HTHRESHcycle(HTHRESHcycleCLEAN)={nan};      
HTHRESHcycle=cell2mat(HTHRESHcycle);

% take every 6th value starting with 1, 2, etc to get a matrix with all
% threshold values corresponding to a particular input log copy #
HTHRESH0  = HTHRESHcycle(1:ts:end); 
HY = zeros(1, tr*tf*(ts-1));
for i = 1:(ts-1)
    HY((tr*tf*(i-1)+1):tr*tf*i) = HTHRESHcycle((i+1):ts:end);
end

% build an array with X values for each of the HY values

for i = 2:(ts-1)
    f(i) =f(i-1)+step;
end

x=zeros((ts-1), tr*tf);
for i=1:tr*tf                 
    x(:,i) = f;
end

xarray = zeros(1, (ts-1)*tr*tf);
for i = 1:(ts-1)
    xarray(((i-1)*tr*tf+1):(i*tr*tf)) = x(i,:);
end

% clean the vector to get rid of x spots corresponding to HY NaNs
xclean = ~isnan(HY).*xarray;
xclean(xclean == 0) = [];
xfit = [0:0.1:5];
% get rid of NaN in HY
HY(find(isnan(HY)))=[];


%% BEST FIT: EXPONENTIAL
exp_fit = fit(xclean', HY', 'exp1');
% use fit to calculate expected values
ypred = exp_fit.a*exp(exp_fit.b*xclean);
% use expected values to calculate R-squared
SSE = sum((HY-ypred).^2);
SST = sum((HY-mean(HY)).^2); 
R2_EXPONENTIAL = 1-SSE/SST

% make bestfit line for plotting later
xfit = [0:0.1:6];
yline = exp_fit.a*exp(exp_fit.b*xfit);

% Process Data and Write to Excel
HT = num2str(HThresh);
m = exp_fit.a;
b = exp_fit.b;
exceldata = [m, b, posTHRESH, HThresh, FAMthresh, f(1), step];
headers = {'m exp_fit.a','b exp_fit.b', 'posTHRESH', 'HThresh', 'FAMthresh', 'lowest concentration', 'step'};
excelfilename = strcat('STD_CURVE_thresh', HT, '.xlsx');
xlswrite(excelfilename,headers,1,'A1');
xlswrite(excelfilename,exceldata,1,'A2');


%% BEST FIT: LINEAR
LinearFit = polyfit(xclean', HY', 1);
m = LinearFit(1,1);
b = LinearFit(1,2);
% use fit to calculate expected values
LINEARpred = m*xclean+b;
% use expected values to calculate R-squared
SSE_LINEAR = sum((HY-LINEARpred).^2);
SST_LINEAR = sum((HY-mean(HY)).^2); 
R2_LINEAR = 1-SSE_LINEAR/SST_LINEAR


%% BEST FIT: 2nd ORDER POLYNOMIAL
LinearFit = polyfit(xclean', HY', 2);
a2 = LinearFit(1,1);
a1 = LinearFit(1,2);
a0 = LinearFit(1,3);
% use fit to calculate expected values
secondORDERpred = a2*xclean.^2+a1*xclean+a0;
% use expected values to calculate R-squared
SSE_SECOND = sum((HY-secondORDERpred).^2);
SST_SECOND = sum((HY-mean(HY)).^2); 
R2_SECONDORDER = 1-SSE_SECOND/SST_SECOND


%%
close all
hold on
plot(xclean,HY,'bo')
p = plot(xfit, yline, 'k-', 'LineWidth', 2);
axis([.1 f(end) 0 60])
title('log10 copies HIV vs cycle threshold');
m = exp_fit.a;
b = exp_fit.b;
ExpFit = ['y = ' num2str(m) ' * exp (' num2str(b) ' * x)'];
text(.25,5, ExpFit)