Skip to content

GitLab

Commit 81f96ac8 authored by drobot's avatar drobot
Browse files

Initial commit

parents
Hide whitespace changes
Inline Side-by-side
.DS_Store 0 → 100644
View file @ 81f96ac8
File added
FRAP/.DS_Store 0 → 100644
View file @ 81f96ac8
File added
FRAP/BuildAndFit.m 0 → 100644
View file @ 81f96ac8
%%%%%%%%%%%%%%Build And Measure
clear all; close all; clc
AnzExp=1; % You can choose between monoexponential (1) or biexponential (2) fitting
timeMax=0; % Maximal time for analysis (0 - all time points were analysed)
%% get Names and indices of the files in the current folder (produced by Fiji using the Frap_in_Fiji.imj plugin)
endung='txt';
Names2=['*.' endung];
FileNames=dir(Names2);
[~,index] = sortrows({FileNames.name}.'); FileNames = FileNames(index);
for i=1:size(FileNames)
sampleorder{i,1}=FileNames(i).name;
end
idx=find(sampleorder{1}=='_');
idx=idx(1);
AllNames=sampleorder{1}(idx+1:end);
i=1;
while sampleorder{i}(1)=='I'
AllNamesidx=sampleorder{i}(idx+1:end);
Num = regexp(AllNamesidx,'\d');
Names(i)=evalin('base',sprintf(sampleorder{i}((idx)+Num)));
i=i+1;
end
i=i-1;
idx=find(sampleorder{i}=='_');
AllNames=sampleorder{i}(idx+1:end);
Num = regexp(AllNamesidx,'\d');
AllNames(end-7:end)=[];
AllNames(Num)=[];
clear i endung FileNames Sortierung Num AllNamesidx sampleorder Names2 idx index
%% User input for defining the type FRAP experiment
display('Which type of FRAP')
FRAPtype = input('spot in droplet (1), bulk phase (2), whole droplet (3): ');
%% Read and normalise data
for i=1:length(Names)
assignin ('base',['Image' num2str(Names(i))], dlmread(sprintf(['timeSeries_' AllNames num2str(Names(i)) '_Out.txt']),'\t',1,0));
if timeMax==0
tend=length(dlmread(sprintf(['timeSeries_' AllNames num2str(Names(i)) '_Out.txt']),'\t',1,0));
else
assignin ('base',['time'], evalin('base',sprintf(['Image' num2str(Names(i)) '(1:end,1)'])));
tend=find(time>timeMax);
if isempty(tend)
tend=length(dlmread(sprintf(['timeSeries_' AllNames num2str(Names(i)) '_Out.txt']),'\t',1,0));
else
tend=tend(1);
end
end
assignin ('base',['Raw' num2str(Names(i))], evalin('base',sprintf(['Image' num2str(Names(i)) '(1:tend,2:end)'])));
timetrace=evalin('base',sprintf(['Raw' num2str(Names(i)) '(1:tend,1)']));
for j=2:size(timetrace,1)
timetraceDiff(j)=timetrace(j)-timetrace(j-1);
end
t0(i)=find(timetraceDiff==min(timetraceDiff));
%%%Normalization
if FRAPtype==3 % whole droplet
if ~isnan(evalin('base',sprintf(['Raw' num2str(Names(i)) '(1,3)'])))
assignin ('base',['RawBG' num2str(Names(i))], evalin('base',sprintf(['Raw' num2str(Names(i)) '-Raw' num2str(Names(i)) '(:,3)' ])));
else
assignin ('base',['RawBG' num2str(Names(i))], evalin('base',sprintf(['Raw' num2str(Names(i))])));
end
assignin ('base',['NormRef' num2str(Names(i))], evalin('base',sprintf(['RawBG' num2str(Names(i)) './RawBG' num2str(Names(i)) '(:,2)' ])));
assignin ('base',['NormDrop' num2str(Names(i))], evalin('base',sprintf(['NormRef' num2str(Names(i))])));
elseif FRAPtype==2 % bulk phase
assignin ('base',['RawBG' num2str(Names(i))], evalin('base',sprintf(['Raw' num2str(Names(i))])));
assignin ('base',['NormRef' num2str(Names(i))], evalin('base',sprintf(['RawBG' num2str(Names(i))])));
assignin ('base',['NormDrop' num2str(Names(i))], evalin('base',sprintf(['NormRef' num2str(Names(i))])));
elseif FRAPtype==1 % spot in droplet ==> double normalization
assignin ('base',['RawBG' num2str(Names(i))], evalin('base',sprintf(['Raw' num2str(Names(i)) '-Raw' num2str(Names(i)) '(:,4)' ])));
assignin ('base',['NormRef' num2str(Names(i))], evalin('base',sprintf(['RawBG' num2str(Names(i)) './RawBG' num2str(Names(i)) '(:,3)' ])));
assignin ('base',['NormDrop' num2str(Names(i))], evalin('base',sprintf(['NormRef' num2str(Names(i)) './ NormRef' num2str(Names(i)) '(:,2)' ])));
end
assignin ('base',['timetrace' num2str(Names(i))], evalin('base',sprintf(['NormDrop' num2str(Names(i)) '(1:tend,1)/mean( NormDrop' num2str(Names(i)) '(1:t0(i)-1,1))' ])));
assignin ('base',['Diam' num2str(Names(i))], dlmread(sprintf(['Measurements_' AllNames num2str(Names(i)) '_Out.txt']),'\t',[1 8 1 8]));
assignin ('base',['time' num2str(Names(i))], evalin('base',sprintf(['Image' num2str(Names(i)) '(1:tend,1)'])));
end
%% call Fit
for i=1:length(Names)
Data=evalin('base',sprintf(['timetrace' num2str(Names(i))]));
timeuse=evalin('base',sprintf(['time' num2str(Names(i))]));
[Fit,Parameter,Results]=FitFRAP(timeuse,Data,AnzExp);
assignin ('base',['Fit' num2str(Names(i))], evalin('base',sprintf(['Fit'])));
assignin ('base',['Parameter' num2str(Names(i))], evalin('base',sprintf(['Parameter'])));
assignin ('base',['Results' num2str(Names(i))], evalin('base',sprintf(['Results'])));
end
%% call plotting
FRAPplot
%%
clear timeMax
FRAP/FRAPplot.m 0 → 100644
View file @ 81f96ac8
% % % % % AverageFrap
showparam=1; % 1 for displaying parameters in the 'Command window'
if FRAPtype==2 || FRAPtype==1
doDiffusion=1; %calculating diffusion coefficients
else
doDiffusion=0; %no calculation of diffusion coefficients
end
sample='FAMsubstrate_PP_'; %% defines variable names for diffusion coefficients
close all
set(0,'DefaultTextInterpreter','tex') %%%for plotting subscripts
%% defining colors
cmap= [0 0.4470 0.7410 0.3
0.8500 0.3250 0.0980 0.3
0.9290 0.6940 0.1250 0.3
0.4940 0.1840 0.5560 0.3
0.4660 0.6740 0.1880 0.3
0.3010 0.7450 0.9330 0.3
0.6350 0.0780 0.1840 0.3];
%% getting the time of the bleached frame
alltimesbleach=[];
for i=1:length(Names)
assignin ('base', 't0',evalin('base',sprintf(['Results' num2str(Names(i)) '.t0(1)' ]) ))
assignin ('base', 'alltimesbleach',[alltimesbleach t0]);
end
idx=find(alltimesbleach==max(alltimesbleach));
toff=alltimesbleach(idx(1))-alltimesbleach;
%% plotting raw data
figure(1)
hold on
mintimeuse=[];
maxtimeuse=[];
for i=1:length(Names)
assignin ('base','timeuse', evalin('base',sprintf(['time' num2str(Names(i))]))-alltimesbleach(i));
mintimeuse=min([mintimeuse timeuse']);
maxtimeuse=max([maxtimeuse timeuse']);
plot(evalin('base',sprintf(['timeuse'])),evalin('base',sprintf(['timetrace' num2str(Names(i))])),'.','color',[0.5 0.5 0.5],'linewidth',1);
end
axis([mintimeuse maxtimeuse 0 1.1]);
tau1=[];
A1=[];
c=[];
if AnzExp==2
tau2=[];
A2=[];
end
%% extracting parameters from different Results
for i=1:length(Names)
assignin('base',['c'],[c evalin('base',sprintf('Results%d.c', [Names(i)]))]);
assignin('base',['tau1'],[tau1 evalin('base',sprintf('Results%d.tau1', [Names(i)]))]);
assignin('base',['A1'],[A1 evalin('base',sprintf('Results%d.A1', [Names(i)]))]);
if AnzExp==2
assignin('base',['tau2'],[tau2 evalin('base',sprintf('Results%d.tau2', [Names(i)]))]);
assignin('base',['A2'],[A2 evalin('base',sprintf('Results%d.A2', [Names(i)]))]);
end
end
cMean=mean(c(1,:));
cStd=std(c(1,:));
A1Mean=mean(A1(1,:));
A1Std=std(A1(1,:));
if AnzExp==2
A2Mean=mean(A2(1,:));
A2Std=std(A2(1,:));
tau2Mean=mean(tau2(1,:));
tau2Std=std(tau2(1,:));
end
tau1Mean=mean(tau1(1,:));
tau1Std=std(tau1(1,:));
y0=0;
%% calculating the Diffusion coefficient
if doDiffusion==1
sample1=['Diff_' sample];
sample2=['Diff_' sample 'All'];
assignin ('base', sprintf(sample2),[]);
for i=1:length(Names)
Data=evalin('base',sprintf('timetrace%d', [Names(i)]));
Fit=evalin('base',sprintf('Fit%d', [Names(i)]))';
assignin ('base', sprintf([sample1 num2str(Names(i))] ),0.88.*(evalin('base',sprintf('Diam%d', [Names(i)]))/2).^2/(4*evalin('base',sprintf('Results%d.tau1(1)', [Names(i)]))*log(2)));
assignin ('base', sprintf([sample2] ), [evalin('base',sprintf([sample2])) evalin('base',sprintf([sample1 num2str(Names(i))] ))]);
end
assignin ('base', sprintf([sample2 '_Mean']),mean(evalin('base',sprintf([sample2] ))));
assignin ('base', sprintf([sample2 '_STD']),std(evalin('base',sprintf([sample2] ))));
end
%% calculating the Fits
clear testFit FitAll
tvAll=mintimeuse:0.01:maxtimeuse;
if AnzExp==1
for i=1:length(tvAll)
if tvAll(i)<0
Exp1=0;
FitAll(i)=(1-y0)-Exp1;
else
Exp1=A1Mean*exp(-(tvAll(i))/tau1Mean);
FitAll(i)=(1-cMean-y0)-Exp1;
end
end
elseif AnzExp==2
for i=1:length(time)
if tvAll(i)<0
Exp1=0;
FitAll(i)=(1-y0)-Exp1;
else
Exp1=A1Mean*exp(-(tvAll(i))/tau1Mean)+A2Mean*exp(-(tvAll(i))/tau2Mean);
FitAll(i)=(1-cMean-y0)-Exp1;
end
end
end
clear Exp1 Exp2
%% Averaging fits for plotting
for i=1:length(Names)
Data=evalin('base',sprintf('timetrace%d', [Names(i)]));
Fit=evalin('base',sprintf('Fit%d', [Names(i)]))';
assignin ('base', sprintf('Error%d', [Names(i)]),Data-Fit);
end
%% plotting average fits
plot(tvAll,FitAll,'linewidth',1,'color',cmap(1,1:3))
xlabel('Time (s)');
ylabel('I_{norm}');
% inserting the residual plot
axes('Position',[.6 .45 .30 .3])
box on
minError=[];
maxError=[];
for i=1:length(Names)
assignin ('base','timeuse', evalin('base',sprintf(['time' num2str(Names(i))]))-alltimesbleach(i));
plot(timeuse, evalin('base',sprintf(['Error' num2str(Names(i))])),'.','color',[0.5 0.5 0.5]);hold on
minError=min([minError evalin('base',sprintf(['Error' num2str(Names(i))]))' ]);
maxError=max([maxError evalin('base',sprintf(['Error' num2str(Names(i))]))' ]);
end
axis([mintimeuse maxtimeuse roundsd(minError*1.1,1) roundsd(maxError*1.1,1)])
yticks([roundsd(minError,1) 0 roundsd(maxError,1)])
set(gca,'fontsize',7)
hold off
%% display parameters
if showparam==1
clc
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('mean parameters +/- standard deviation')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
if doDiffusion==1
disp([sprintf(['Diffusion coefficient = ' num2str(roundsd(mean(evalin('base',sprintf([sample2] ))),3)) ' +/- ' num2str(roundsd(std(evalin('base',sprintf([sample2] ))),2)) ' ']) char(956) 'm' char(178) '/s'])
else
disp([char(964) sprintf([' = ' num2str(roundsd(tau1Mean,2)) ' +/- ' num2str(roundsd(tau1Std,2)) ' s'])])
end
if AnzExp==2
disp([char(964) sprintf(['2 = ' num2str(roundsd(tau2Mean,2)) ' +/- ' num2str(roundsd(tau2Std,2)) ' s'])])
disp(sprintf(['Bleached Fraction = ' num2str(roundsd((A1Mean+A2Mean+cMean)*100,2)) ' +/- ' num2str(roundsd((A1Std+A2Mean+cStd)*100,2)) ' %%']))
else
disp(sprintf(['Bleached Fraction = ' num2str(roundsd((A1Mean+cMean)*100,2)) ' +/- ' num2str(roundsd((A1Std+cStd)*100,2)) ' %%']))
end
disp(sprintf(['Immobile Fraction = ' num2str(roundsd(cMean*100,2)) ' +/- ' num2str(roundsd(cStd*100,2)) ' %%']))
end
%%
clear alltimesbleach tvAll idx maxError minError i k mintimeuse maxtimeuse idx Data Fit j Parameter Results sample sample1 sample2 tend
clear timetrace timetraceDiff timeuse cmap FitAll showparam doDiffusion
\ No newline at end of file
FRAP/FitFRAP.m 0 → 100644
View file @ 81f96ac8
function [Fit,Parameters,Results]=FitFRAP(time,timetrace,AnzExp,t0)
close all
%%%%%%%%%%%%%%%%%%%%%Unhindered Diffusion should be monoexponetial
%%%%%%%%%%%%%%%%%%%%But interaction with Polymers could yield
%%%%%%%%%%%%%%%%%%%%%biexponential
%%%%time = vector with timesteps
%%%%timetrace = vector with normalized data
%%%%AnzExp = number of exponentials
%%%%optimize an initial offset from 1 (1 for optimize)
doy0=0;
%%%%show fit evolution (1 for show)
doshow=0;
%%%%weighting of the curve (the value gives the number of points
%%%%after the bleached frame to get a higher wighting (5 x))
weight=0;
%% find t0
for i=2:size(timetrace,1)
timetraceDiff(i)=timetrace(i)-timetrace(i-1);
end
dot0=find(timetraceDiff==min(timetraceDiff));
%% set fitting options
set(0,'defaulttextinterpreter','latex')
clear J R ci se h p st chi2red beta CovB MSE parametererror dmono datymono datymono xmono x0mono Spectra A1Fakt A2Fakt A3Fakt A4Fakt
opts = optimoptions(@lsqcurvefit,'Algorithm','levenberg-marquardt',...
'display','off',...
'ScaleProblem','jacobian',...
'TolFun',10^-6,...
'TolX',10^-6);
%Constrains (had to be empty for levenberg marquardt)
lb=[];
ub=[];
%%%Test if 'Statistics Toolbox' is available
%%%(needed for error estimation with nlparci)
testToolbox=license('checkout','Statistics_Toolbox');
if testToolbox==1
nltest=1;
else
nltest=0;
end
%% Initial guesses
%%%%%%Automatical estimation of final offset
c2guess=1-mean(timetrace(end-5:end));
%%%%%%Automatical estimation of tau1
timetraceRec=timetrace(dot0:end);
idx=find(timetraceRec>timetraceRec(1)+(timetraceRec(end)-timetraceRec(1))/1.5);
if any(idx>0)
tau1guess=time(idx(1));
else
tau1guess=time(round(length(time)/3));
end
%%%%%%Automatical estimation of the amplitudes
if AnzExp==1
A1guess=1-timetrace(dot0)+c2guess;
else
A1guess=(1-timetrace(dot0)+c2guess)/2;
end
if AnzExp==1
A2guess=1;
else
A2guess=(1-timetrace(dot0)+c2guess)/2;
end
%%%%%%Guessing tau2
tau2guess=6*tau1guess;
y0guess=0;
t0guess=time(dot0);
%StartValues=[tau1 A1 c2 A2 tau2 y0 t0]
StartValues= [tau1guess A1guess c2guess A2guess tau2guess y0guess t0guess];
if dot0~=1
StartValues(7)=[];
end
if doy0~=1
StartValues(6)=[];
end
if AnzExp==1
StartValues(4:5)=[];
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Parameter fitting
if weight>0
ToFit=zeros(size(timetrace));
else
ToFit=timetrace;
end
StdError=[];
if AnzExp==1
FaltFourExpGlob=@(x,datx) FaltFourExpFuncGlob1y0(x,datx,AnzExp,timetrace,doshow,doy0,dot0,weight);
[Param,resnorm,residual,exitflag,output,lambda,J]=lsqcurvefit(FaltFourExpGlob,StartValues,time,ToFit,lb,ub,opts);
if nltest==1
[ci] = nlparci(Param,residual,'jacobian',J);
StdError=Param-ci(:,1)';
StdError=roundsd(StdError,2)';
end
elseif AnzExp==2
FaltFourExpGlob=@(x,datx) FaltFourExpFuncGlob2y0(x,datx,AnzExp,timetrace,doshow,doy0,dot0,weight);
[Param,resnorm,residual,exitflag,output,lambda,J]=lsqcurvefit(FaltFourExpGlob,StartValues,time,ToFit,lb,ub,opts);
if nltest==1
[ci] = nlparci(Param,residual,'jacobian',J);
StdError=Param-ci(:,1)';
StdError=roundsd(StdError,2)';
end
end
%% calculating the fit using optimized parameters
Parameters=abs(Param);
if AnzExp==1
x=Param;
tau1=x(1);
A1=x(2);
c2=x(3);
if doy0==1
y0=x(4);
else
y0=doy0;
end
if dot0==1
t0=x(5);
else
t0=time(dot0);
end
for i=1:length(time)
if time(i)<t0
Exp1=0;
Fit(i)=(1-y0)-Exp1;
else
Exp1=A1*exp(-(time(i)-t0)/tau1);
Fit(i)=(1-c2-y0)-Exp1;
end
end
elseif AnzExp==2
x=Param;
tau1=x(1);
A1=x(2);
c2=x(3);
A2=x(4);
tau2=x(5);
if doy0==1
y0=x(6);
else
y0=doy0;
end
if dot0==1
t0=x(7);
else
t0=time(dot0);
end
for i=1:length(time)
if time(i)<t0
Exp1=0;
Fit(i)=(1-y0)-Exp1;
else
Exp1=A1*exp(-(time(i)-t0)/tau1)+A2*exp(-(time(i)-t0)/tau2);
Fit(i)=(1-c2-y0)-Exp1;
end
end
end
%% putting the standard error in the second row of the parameters
if any(StdError>0)
Parameters=[Parameters;StdError'];
end
%% Generating a results table
Results=table([0; 0],[0; 0],[0; 0],[0; 0],[0; 0],[0; 0],[0; 0]);
Results.Properties.VariableNames = {'tau1' 'A1' 'c' 'tau2' 'A2' 'y0' 't0'};
Results.Properties.RowNames = {'Value' 'Std'};
Results.tau1=Parameters(:,1);
Results.A1=Parameters(:,2);
Results.c=Parameters(:,3);
if AnzExp==1
if doy0==1
Results.y0=Parameters(:,4);
else
Results.y0=[doy0;0];
end
if dot0==1
Results.t0=Parameters(:,5);
else
Results.t0=[time(dot0);0];
end
elseif AnzExp==2
Results.A2=Parameters(:,5);
Results.tau2=Parameters(:,4);
if doy0==1
Results.y0=Parameters(:,6);
else
Results.y0=[doy0;0];
end
if dot0==1
Results.t0=Parameters(:,7);
else
Results.t0=[time(dot0);0];
end
end
end
function F=FaltFourExpFuncGlob1y0(x,datx,AnzExp,timetrace,doshow,doy0,dot0,weight)
tau=x(1);
A1=x(2);
c2=x(3);
if doy0==1
y0=x(4);
else
y0=doy0;
end
if dot0==1
t0=x(5);
else
t0=datx(dot0);
end
for i=1:length(datx)
if datx(i)<t0
Exp1=0;
F(i)=(1-y0)-Exp1;
else
Exp1=A1*exp(-(datx(i)-t0)/tau);
F(i)=(1-c2-y0)-Exp1;
end
end
F=F';
% doshow
if doshow==1
plot(datx,timetrace,'o','markers',12);hold on
plot(datx,F,'linewidth',2)
drawnow; hold off
end
if weight>0
F=timetrace-F;
F(dot0:dot0+weight)=F(dot0:dot0+weight)*5;
end
end
function F=FaltFourExpFuncGlob2y0(x,datx,AnzExp,timetrace,doshow,doy0,dot0,weight)
tau1=x(1);
A1=x(2);
c2=x(3);
A2=x(4);
tau2=x(5);
if doy0==1
y0=x(6);
else
y0=doy0;
end
if dot0==1
t0=x(7);
else
t0=datx(dot0);
end
for i=1:length(datx)
if datx(i)<t0
Exp1=0;
F(i)=(1-y0)-Exp1;
else
Exp1=A1*exp(-(datx(i)-t0)/tau1)+A2*exp(-(datx(i)-t0)/tau2);
F(i)=(1-c2-y0)-Exp1;
end
end
F=F';