Contents
Function, exponential curve fit
%------------------------------------------------------------ % function [Curvefit, goodness] = myfit2(MeanInt,FrameTreat) % [Curvefit, goodness] = myfity2(MeanInt,FrameTreat) % Curve Fit % % Input % MeanInt - curve to be fitted must be exponential rising curve % FrameTreat - frame number corresponding to time of challenge/treatment % % Output % Curvefit - curve fit to data % goodness - goodness of fit least square error % % Curve Fit Equation % y(x) = alpha*[1 -exp(-beta *(x-r))] + c % Solutions: % y[ x = (1/beta) + r ] = alpha( 1 - 1/e ) + c % = alpha(0.63)+c % c = ymin % alpha = ymax - c % beta: alpha(1-1/e)+c = alpha[1-exp(-beta*(x-r))]+c % alpha(1-1/e) = alpha[1-exp(-beta*(x-r))] % (1-1/e) = 1-exp[-beta*(x-r)] % -1/e = -exp[-beta*(x-r)] % ln(1/e) = -beta*(x-r) % ln(1)-ln(e) = -beta*(x-r) % 1 = beta*(x-r) % beta = 1/(x-r) % % | % alpha+c ---> | . . . % | . % | . %0.63(alpha)+c ---> | . % | . % | . % | . % | . % c ---> | . % | % |__________________________________________ MeanInt = MeanInt(FrameTreat:FrameTreat-1+36); %3 min post treat, inclusive treat SMeanInt = smooth(MeanInt,'lowess'); %smooth the curve dSMeanInt = diff(SMeanInt); %the discrete derivative of smooth curve zerocross = diff(sign(dSMeanInt)); %the zero crossing of derivative if ~isempty(find(zerocross<=-1,1)) %get peak if zero cross fitlim = find(zerocross<=-1); %all indices of zero cross [~, findx] = max(SMeanInt(fitlim)); %max value of original smooth fx at zero corssings fitlim = fitlim(findx)+1; %index of max value zero cross else %get peak if no zero cross %[~, fitlim] = min(abs(dSMeanInt)); fitlim = 12; %1min post treat end TreatCurve = MeanInt(1+1:fitlim);%curve post treat %TreatCurve data index starts at 1 a = ones(1,length(TreatCurve)); %make array same length as TreatCurve a = a*TreatCurve(end); %make array equal to platue value TreatCurve = [TreatCurve a]; %concatenate, extend platue x = 0:(length(TreatCurve)-1); %start x-axis from zero x = x*5; %convert from frames to seconds (1frame = 5seconds) x = x+(5*FrameTreat); %shift the curve to match time stamp myc = round(min(TreatCurve)); myalpha = round(max(TreatCurve)-myc); temp = 0.63*myalpha+myc; %find value 63% temp = abs(TreatCurve-temp); %subtract value 66% of max from curve [~, tempx] = min(temp); %get index of min difference = index of 63%myalpha+myc mybeta = 1/(tempx); r = 5*FrameTreat; %create right shift known constant (seconds) s = fitoptions('Mehtod','NonlinearLeastSquares','Robust','LAR',... 'Lower',[myalpha-.01,0,myc-.01],... 'Upper',[myalpha+.01,inf,myc+.01],... 'Startpoint',[myalpha,mybeta,myc]); %set the fitting parameters f = fittype( @(alpha,beta,c,x) alpha*(1 - exp(-beta *(x-r))) +c,... 'coefficients',{'alpha','beta','c'},'options',s ); %set the curve fit function [Curvefit, goodness] = fit(x',TreatCurve',f); %fit data to an eponential rise curve figure, plot(Curvefit,'k-'),hold on %create a plot of the curve fit and original data plot(x, TreatCurve,'k.','markersize',12) %plot data points legend('Fit Curve', 'Data','Location','SouthEast') %place legend xlabel('Time(seconds)','FontSize',20) %x-axis label ylabel('Fluorescence Intensity','FontSize',20) %y-axis label Txticks = get(gca,'XTick'); %x-axis tick array Pxticks = round(length(get(gca,'XTick'))/4); %1/4 length of x-axis tick array xticks = Txticks(Pxticks); %1/4 length value of x-axis ticks tempytick = abs(x-xticks); %x-axis data - x-axis 1/4 length tick value [~, tempytick] = min(tempytick); %find position of min difference text(xticks,TreatCurve(tempytick),{['y(x) = ','\alpha',... ' * [ 1 - exp( -','\beta',' * (x-r) ) ] + c'];... ['y(x) = ',num2str(Curvefit.alpha,4),... '*[1-exp(-',num2str(Curvefit.beta,4),'*(x-',num2str(5*FrameTreat,4),... '))]+',num2str(Curvefit.c,4)];['R = ', num2str(goodness.rsquare,4)];... ['\tau = 1/\beta = ', num2str(1/Curvefit.beta,4), ' (seconds)']},... 'VerticalAlignment','Cap','FontSize',16) %place equation text %------------------------------------------------------------
Function, Outlier Test
%------------------------------------------------------------ % function [ObjectOut alphas taus] = myout3(meanInt,FrameTreat) % [ObjectOut, alphas, taus] = myfity2(MeanInt,FrameTreat) % Outlier discovery based on curve fit parameters and on outlier criteria % % Input % meanInt - [#cells,1,#frames] curve to be fitted must be exponential rising curve % FrameTreat - frame number corresponding to time of challenge/treatment % % Output % ObjectOut - numeric 1 or 0 indicating exclusion status based on outlier criteria % alphas, taus - curve fit parameters objects = size(meanInt,1); frames = size(meanInt,3); curves = reshape(meanInt,objects,frames); for k = 1:objects curve = curves(k,:); curve = curve(FrameTreat:FrameTreat-1+36); %3 min post treat, inclusive treat SMeanInt = smooth(curve,'lowess'); %smooth the curve dSMeanInt = diff(SMeanInt); %the discrete derivative of smooth curve zerocross = diff(sign(dSMeanInt)); %the zero crossing of derivative if ~isempty(find(zerocross<=-1,1)) %get peak if zero cross fitlim = find(zerocross<=-1); %all indices of zero cross [~, findx] = max(SMeanInt(fitlim)); %max value of original smooth fx at zero corssings fitlim = fitlim(findx)+1; %index of max value zero cross else %get peak if no zero cross %[~, fitlim] = min(abs(dSMeanInt)); fitlim = 12; %get 1min post treat end TreatCurve = curve(1+1:fitlim);%curve post treat %TreatCurve data index starts at 1 a = ones(1,length(TreatCurve)); %make array same length as TreatCurve a = a*TreatCurve(end); %make array equal to platue value TreatCurve = [TreatCurve a]; %concatenate, extend platue x = 0:(length(TreatCurve)-1); %start x-axis from zero x = x*5; %convert from frames to seconds (1frame = 5seconds) x = x+(5*FrameTreat); %shift the curve to match time stamp myc = round(min(TreatCurve)); myalpha = round(max(TreatCurve)-myc); temp = 0.63*myalpha+myc; %find value 63% temp = abs(TreatCurve-temp); %subtract value 66% of max from curve [~, tempx] = min(temp); %get index of min difference = index of 63%myalpha+myc mybeta = 1/(tempx); r = 5*FrameTreat; %create right shift known constant (seconds) s = fitoptions('Mehtod','NonlinearLeastSquares','Robust','LAR',... 'Lower',[myalpha-.01,0,myc-.01],... 'Upper',[myalpha+.01,inf,myc+.01],... 'Startpoint',[myalpha,mybeta,myc]); %set the fitting parameters f = fittype( @(alpha,beta,c,x) alpha*(1 - exp(-beta *(x-r))) +c,... 'coefficients',{'alpha','beta','c'},'options',s ); %set the curve fit function [Curvefit, goodness] = fit(x',TreatCurve',f); %fit data to an eponential rise curve alphas(k) = Curvefit.alpha; taus(k) = 1/Curvefit.beta; % %plot % figure, plot(Curvefit,'k-'),hold on %create a plot of the curve fit and original data % plot(x, TreatCurve,'k.','markersize',12) % legend('Fit Curve', 'Data','Location','SouthEast') % xlabel('Time(seconds)','FontSize',20) % ylabel('Fluorescence Intensity','FontSize',20) % text(x(2),0.95*(TreatCurve(2)),{['y(x) = ','\alpha',... % ' * [ 1 - exp( -','\beta',' * (x-r) ) ] + c'];... % ['y(x) = ',num2str(Curvefit.alpha,4),... % '*[1-exp(-',num2str(Curvefit.beta,4),'*(x-',num2str(5*FrameTreat,4),... % '))]+',num2str(Curvefit.c,4)];['R = ', num2str(goodness.rsquare,4)];... % ['\tau = 1/\beta = ', num2str(1/Curvefit.beta,4), ' (seconds)']},... % 'VerticalAlignment','Top','FontSize',18) end %Remove outliers. Exclude data 2SEMs away (2SEMs ~ 95% population about mean) meanalphas = mean(alphas); stealphas = std(alphas)/sqrt(length(alphas)); meantaus = mean(taus); stetaus = std(taus)/sqrt(length(taus)); for k = 1:objects if abs(alphas(k)-meanalphas) > 4*stealphas ||... abs(taus(k)-meantaus) > 4*stetaus %95percent population in norm distribution ObjectOut(k) = 1; else ObjectOut(k) = 0; end end