Commit e5f18bca authored by Antonio Politi's avatar Antonio Politi
Browse files

matlab code

parent 4c332351
%%
pari = MO.par;
pari.kappa = 1;
[parf, norm] = MO.fitmodel([10, 20], {'N', 'tauD1'}, MO.par, data, 2, 10000)
[parf, norm] = MO.fitmodel([parf.N, parf.tauD1, parf.tauT, parf.thetaT], {'N', 'tauD1', 'tauT', 'thetaT'}, parf, data, 2, 10000)
pari.kappa = 2;
[parf, norm] = MO.fitmodel([10, 20], {'N', 'tauD1'}, pari, data, 2, 10000)
[parf, norm] = MO.fitmodel([parf.N, parf.tauD1, parf.tauT, parf.thetaT], {'N', 'tauD1', 'tauT', 'thetaT'}, parf, data, 2, 10000)
%%
% absFcsmodel abstract class for fitting FCS data to using the autocorrelation curve that comes from
% Fluctualize M. Wachsmuth
%
% EMBL, Antonio Politi, November 2016
%
classdef absFcsmodel < handle
properties (Access = public, Abstract)
% these are the default parameters for a protein
%par = struct('N', 10, 'thetaT', 0.2, 'tauT', 100, 'f1', 1, 'tauD1', 500, 'alpha1', 1,
%'tauD2', 5000, 'alpha2', 1, 'kappa', 5.5); % parameter of model lb = struct('N', 0.0001,
%'thetaT', 0.001,'tauT', 0.1, 'f1', 1, 'tauD1', 100, 'alpha1', 0.5, 'tauD2', 500, 'alpha2', 0.5,
%'kappa', 1); hb = struct('N', 10000, 'thetaT', 1, 'tauT', 1000, 'f1', 1, 'tauD1',
%50000,'alpha1', 1.2, 'tauD2', 500000, 'alpha2', 2, 'kappa', 20);
par; % parameter values a struct with name of parameter and value. This needs to be defined in the derived class
lb; % low boundary of parameter values
hb; % high boundary of parameter values
end
methods (Abstract)
Gcomp(obj); %model to compute the ACF given the parameters
end
methods
function MO = absFcsmodel(parnames, parvalues, lbvalues, hbvalues)
% ABSFCSMODEL constructor
% MO = absFcsmodel() - just create class MO = absFcsmodel(parnames, parvalues) - create class and
% assign parvalues to parameters
% named in parnames.
% MO = absFcsmodel(parnames, parvalues) - create class and assign parvalues and lower boundary to
% parameters
% named in parnames.
% MO = absFcsmodel(parnames, parvalues)- create class and assign parvalues, lower and higher
% boundary to parameters
% named in parnames.
% parnames: cell array containing name of paramers parvalues: vector with values of parameters
% lbvalues: vector with values of lower bounds hbvalues: vector with values of higher bounds
if nargin < 3 % should this not be 2? TODO
return
end
% update default parameter
for iname = 1:size(parnames, 1)
name = parnames{iname};
assert(isfield(MO.par, name), [name ' Parameter does not exist!'])
setfield(MO.par,name, parvalues(iname));
if nargin < 4
return
end
setfield(MO.lb, name, lbvalues(iname));
if nargin < 5
return
end
setfield(MO.hb, name, hbvalues(iname));
end
end
function [N, tauD1, tauD2] = initialGuess(MO, data, tauStart, tauEnd)
% INITIALGUESS given the AC data returns an initial guess forN, tauD1 of first component, tauD2 of
% second component
% [N, tauD1, tauD2] = MO.initialGuess(data, tauStart, tauEnd)
% data - is [timept, AC] array tauStart - is value in us of time where to start the analysis
% tauEnd - is value in us of time where to end analysis
% Method from Wachsmuth, M., et al. (2015). Nature Biotechnology, 33(4), 384?389.
% https://doi.org/10.1038/nbt.3146
ist= find(data(:,1) >= tauStart); % index of start time
ien = find(data(:,1) >= tauEnd); % index of end time
ist = ist(1);
ien = ien(1);
N = 1/mean(data(ist:ist+5,2)); % get initial guess from average of AC first 5 time pts
N = N + (1-2*rand)*N/8; % randomize initial guess
tauDm = N*trapz(data(ist:ien,1),data(ist:ien,2))/4.5; % compute integral of equation to obtain a gues of average time
% another way of computing tauDm would be by using a linear approximation of AC
tauDm = tauDm + (1-2*rand)*tauDm/8; % randomize tauDm
% assign value of tauD1 and tauD2 so that tauD1 ~ 1/10 of tauD2. 1/tauDm = f1/tauD1 + (1-f1)/tauD2
tauD1 = (9*MO.par.f1+1)/10*tauDm;
tauD2 = (9*MO.par.f1+1)*tauDm;
end
function Gdiff = Gdiffcomp(MO, par, data)
% GDIFFCOMP computes weighted difference of ACF data and model
% par - full parameter struct of model data - Ntx3 matrix of data [time, ACF, ACF_std]
Gdiff = (MO.Gcomp(par, data(:,1)) - data(:,2))./data(:,3);
end
function Gdiff = Gdiffcomp_unw(MO, par, data)
% GDIFFCOMP_UNW computes unweighted difference of ACF data and model
% par - full parameter struct of model data - Ntx3 matrix of data [time, ACF, ACF_std]
Gdiff = (MO.Gcomp(par, data(:,1)) - data(:,2));
end
function G = Gcompfit(MO, parf, names, par, tau)
% GCOMPFIT computes model ACF using parameter specified for fit
% parf - parameter vector names - name of parameter to be changed par - default parameter
% struct tau - vector of ACF times
for i = 1:length(names)
par = setfield(par, names{i}, parf(i));
end
G = MO.Gcomp(par, tau);
end
function Gdiff = Gdiffcompfit(MO, parf, names, par, data)
% GDIFFCOMPFIT computes difference model to ACF using parameter specified for fit
% parf - parameter vector names - name of parameter to be changed par - default parameter
% struct tau - vector of ACF times
for i = 1:length(names)
par = setfield(par, names{i}, parf(i));
end
Gdiff = MO.Gdiffcomp(par, data);
end
function Gdiff = Gdiffcompfit_unw(MO, parf, names, par, data)
% GDIFFCOMPFIT_UNW computes difference model to ACF using parameter specified for fit. Unweighted
% parf - parameter vector names - name of parameter to be changed par - default parameter
% struct tau - vector of ACF times
for i = 1:length(names)
par = setfield(par, names{i}, parf(i));
end
Gdiff = MO.Gdiffcomp_unw(par, data);
end
function parS = parvecToparstruct(MO, parV)
% PARVECTOPARSTRUCT return a struct of type MO.par with values given by parV
% MO.parvecToparstruct(parV) parV - parameter vector of length MO.par
parS = MO.par;
parnames = fieldnames(MO.par);
for i=1:length(parnames)
parS = setfield(parS, parnames{i}, parV(i));
end
end
function plotTrace(MO, par, data)
% PLOTTRACE plot model and exp data in a graph
% par - full parameter struct data - Ntx3 matrix of data [time, ACF, ACF_std]
G = MO.Gcomp(par, data(:,1));
semilogx(data(:,1), G,'-', data(:,1), data(:,2),'o');
end
function [par, norm] = fitmodel_unw(MO, parf, names, par, data, tauStart, tauEnd)
% FITMODEL_UNW fit model to data using unweighted difference
% [par, norm] = MO.fitmodel_unw(parf, names, par, data) - fit in default time interval (set by the
% data) [par, norm] = MO.fitmodel_unw(parf, names, par, data, tauStart) - fit in time intervael
% [tauStart, tauMax] [par, norm] = MO.fitmodel_unw(parf, names, par, data, tauStart, tauEnd) - fit
% in time intervael [tauStart, tauEnd]
% INPUT:
% parf - vector of parameter values for the fit names - cell array with names of parameters.
% par - struct of parameters data - Ntx3 matrix of data [time, ACF, ACF_std] tauStart -
% start time in us for fitting domain tauEnd - end time in us for fitting domain
% OUTPUT:
% par - struct containing fitted parameters norm - sum of squared difference at fitting
% minima
if nargin <6
tauStart = data(1,1);
end
if nargin < 7
tauEnd = data(end,1);
end
ist= find(data(:,1) >= tauStart);
ien = find(data(:,1) >= tauEnd);
ist = ist(1);
ien = ien(1);
lb = [];
hb = [];
for i = 1:length(names)
lb = [lb;getfield(MO.lb,names{i})];
hb = [hb;getfield(MO.hb,names{i})];
end
options = optimset;
[parf, norm] = lsqnonlin(@MO.Gdiffcompfit_unw, parf, lb, hb, options, names, par, data(ist:ien,:));
for i = 1:length(names)
par = setfield(par, names{i}, parf(i));
end
end
function [par, norm] = fitmodel(MO, parf, names, par, data, tauStart, tauEnd)
% FITMODEL_UNW fit model to data using std weighted difference
% [par, norm] = MO.fitmodel(parf, names, par, data) - fit in default time interval (set by the
% data) [par, norm] = MO.fitmodel(parf, names, par, data, tauStart) - fit in time intervael
% [tauStart, tauMax] [par, norm] = MO.fitmodel(parf, names, par, data, tauStart, tauEnd) - fit in
% time intervael [tauStart, tauEnd]
% INPUT:
% parf - vector of parameter values for the fit names - cell array with names of parameters.
% par - struct of parameters data - Ntx3 matrix of data [time, ACF, ACF_std] tauStart -
% start time in us for fitting domain tauEnd - end time in us for fitting domain
% OUTPUT:
% par - struct containing fitted parameters norm - sum of squared difference at fitting
% minima
if nargin <6
tauStart = data(1,1);
end
if nargin < 7
tauEnd = data(end,1);
end
ist= find(data(:,1) >= tauStart);
ien = find(data(:,1) >= tauEnd);
ist = ist(1);
ien = ien(1);
lb = [];
hb = [];
for i = 1:length(names)
lb = [lb;getfield(MO.lb,names{i})];
hb = [hb;getfield(MO.hb,names{i})];
end
options = optimset;
[parf, norm] = lsqnonlin(@MO.Gdiffcompfit, parf, lb, hb, options, names, par, data(ist:ien,:));
for i = 1:length(names)
par = setfield(par, names{i}, parf(i));
end
end
function AC = readZeissfcs(MO, fname)
% READZEISSFCS read AC data that has been generated in Zeiss ZEN software *.fcs files and returns a
% cell array
% MO.readZeissfcs(fname) - read data in fname. This is the .fcs file OUTPUT:
% OUTPUT:
% AC - Cell array containing the AC for each measurment: {time, AC1_1, AC2_2}, {time, AC1_2 ... time is in us. AC has
% been converted to 0
if nargin < 1
fname = 'W:\Shotaro\Live-cell-imaging\HeLa-FCS\160701-HeLa-Nup107\160701_gfpNUP107z26z31\Calibration\Alexa\488nm.fcs';
end
fileID = fopen(fname);
% read whole file at once
C = textscan(fileID, '%s', 'delimiter','\n');
% find AC values
ACm = find(strncmp(C{1}, 'CorrelationArray =',18 ));
for i=1:length(ACm)-1
dim = str2num(C{1}{ACm(i)}(19:end));
AC{i} = cell2mat(cellfun(@str2num, C{1}(ACm(i)+1:ACm(i)+dim(1)), 'un', 0));
% convert to us and AC base 0
AC{i}(:,1) = AC{i}(:,1)*1e6;
AC{i}(:,2:end) = AC{i}(:,2:end)-1;
%define default std column
if size(AC{i},2) == 3
AC{i}(:,4) = AC{i}(:,3);
AC{i}(:,5) = 1;
AC{i}(:,3) = 1;
end
if size(AC{i},2) == 2
AC{i}(:,3) = 1;
end
end
end
function AC = readFAcor(MO, fname)
% READFACOR read AC data that has been generated with Fluctualizer
% MO.readFAcor(fname) -
% fname: name of .cor file. If fname is a cell array than all files are read at once
% output:
% AC: cell array containing the AC as - time AC_Ch1 std_ACh1 ACh2 std_ACh2 XC_Ch1Ch2 std_XC_Ch1Ch2
if nargin < 1
fname = {'W:\Shotaro\Live-cell-imaging\HeLa-FCS\160701-HeLa-Nup107\160701_gfpNUP107z26z31\Calibration\Alexa\488nm_R1_P1_K1_Ch2.zen.kappa55.cor'};
end
if ~iscell(fname)
fname = {fname};
end
for i=1:length(fname)
AC{i} = dlmread(fname{i}, '\t', 1, 0);
end
end
function AC = readFAxml(MO, fname)
% READFAXML read AC, parameters, annotations etc from xml file generated per measurement by FA
end
function writeFAxml(MA, fname, parstruct)
% WRITEFAXML write fits parameters, annotations, and settings xml file generated per measurement by FA
end
end
end
function [average_traces] = average_autocorrelation_traces(traces)
% INPUT:
% traces = autocrrelation traces as a function of time in the
% matrix. Time increases along the rows. All should have the
% same time.
%NOTE: The 1/N fitting parameter is meaningless here since the traces have been rescaled
fittedG0 = traces{3}(1,:); %1/N Scaling parameter from the fits to reduce noise
traces_rescaled_f = traces{3}*((1./fittedG0)'); % Averaging the rescaled data, fits and the errors
traces_rescaled = traces{2}*((1./fittedG0)');
traces_rescaled_er = sum(traces{4},2);
average_traces{1} = traces{1}(:,1);
average_traces{2} = traces_rescaled./size(traces{2},2);
average_traces{3} = traces_rescaled_f./size(traces{3},2);
average_traces{4} = traces_rescaled_er./size(traces{4},2);
end
\ No newline at end of file
function [] = chi_square()
%Calculates the Chi Square value of the the fittings
end
\ No newline at end of file
classdef dyeFcsmodel < absFcsmodel
properties (Access = public)
% these are the default parameters for a dye
par = struct('N', 10, 'thetaT', 0.2, 'tauT', 7, 'f1', 1, 'tauD1', 20, 'alpha1', 1, 'tauD2', 5000, 'alpha2', 1, 'kappa', 5.5);
lb = struct('N', 0.0001, 'thetaT', 0.001,'tauT', 0.1, 'f1', 1, 'tauD1', 1, 'alpha1', 0.5, 'tauD2', 500, 'alpha2', 0.5, 'kappa', 1);
hb = struct('N', 10000, 'thetaT', 1, 'tauT', 15, 'f1', 1, 'tauD1', 50000,'alpha1', 1.2, 'tauD2', 500000, 'alpha2', 2, 'kappa', 20);
end
methods
function MO = dyeFcsmodel(parnames, parvalues, lbvalues, hbvalues)
par = struct('N', 10, 'thetaT', 0.2, 'tauT', 7, 'f1', 1, 'tauD1', 20, 'alpha1', 1, 'tauD2', 5000, 'alpha2', 1, 'kappa', 5.5);
lb = struct('N', 0.0001, 'thetaT', 0.05,'tauT', 0.1, 'f1', 1, 'tauD1', 1, 'alpha1', 0.5, 'tauD2', 500, 'alpha2', 0.5, 'kappa', 1);
hb = struct('N', 10000, 'thetaT', 1, 'tauT', 15, 'f1', 1, 'tauD1', 50000,'alpha1', 1.2, 'tauD2', 500000, 'alpha2', 2, 'kappa', 20);
if nargin < 1
parnames = fieldnames(par);
end
if nargin < 2
parvalues = struct2array(par);
end
if nargin < 3
lbvalues = struct2array(lb);
end
if nargin < 4
hbvalues = struct2array(hb);
end
MO@absFcsmodel(parnames, parvalues, lbvalues, hbvalues)
end
function G = Gcomp(MO, par, tau)
% GCOMP returns ACF for dye like blinking
% par - struct of parameter values
% tau - time points to evaluate function
triplet = (1 - par.thetaT + par.thetaT*exp(-tau/par.tauT))/(1-par.thetaT);
comp1 = ((1+(tau/par.tauD1).^par.alpha1).^-1).*((1+1/par.kappa^2*(tau/par.tauD1).^par.alpha1).^-0.5);
comp2 = ((1+(tau/par.tauD2).^par.alpha2).^-1).*((1+1/par.kappa^2*(tau/par.tauD2).^par.alpha2).^-0.5);
G = 1/par.N*triplet.*(par.f1*comp1 + (1-par.f1)*comp2);
end
function vol = computeVolume(MO, par, dyetype)
% COMPUTEVOLUME Compute focal volume return value in picoliter
% par - vector or struct of parameters
if dyetype == 1
D = 464; %Alexe488 at 37 degrees
end
if dyetype == 2
D = 521; %Alexe568 at 37 degrees
end
if isstruct(par)
vol = pi^(3/2).*((2*sqrt(D*par(:,5))/1000).^3).*par(:,9);
else
vol = pi^(3/2).*((2*sqrt(D*par.tauD1)/1000).^3).*par.kappa;
end
end
function w0 = computeRadius(MO, par, dyetype)
% COMPUTERADIUS Compute focal radius ans return value in um
% par - vector or struct of parameters
if dyetype == 1
D = 464; %Alexe488 at 37 degrees
end
if dyetype == 2
D = 521; %Alexe568 at 37 degrees
end
if isstruct(par)
w0 = 2*sqrt(D*par.tauD1)/1000;
else
w0 = 2*sqrt(D*par(:,5))/1000;
end
end
%% fit thetaTtauT for a fixed kappa value
function [outPar, norm] = fitThetaTtauT(MO, data, kappa, weight)
if nargin < 3
kappa = 5;
end
if nargin < 4
weight = 1;
end
%first find optimal tauT
pari = MO.par;
pari.kappa = kappa;
%find initial value for tauD1 and N
for itry = 1:5
[N, tauD1, tauD2] = MO.initialGuess(data, 2, 10000);
if weight
[parf{itry}, normT(itry)] = MO.fitmodel([N, tauD1], {'N', 'tauD1'}, pari, data, 2, 10000);
else
[parf{itry}, normT(itry)] = MO.fitmodel_unw([N, tauD1], {'N', 'tauD1'}, pari, data, 2, 10000);
end
end
[val, pos] = min(normT);
parf1 = parf{pos};
%fit thetaT and tauT
for itry = 1:5
if weight
[parf{itry}, normT(itry)] = MO.fitmodel([parf1.N, parf1.tauD1, pari.thetaT+(1-2*rand)*pari.thetaT/8, pari.tauT+(1-2*rand)*pari.tauT/8], {'N', 'tauD1', 'thetaT', 'tauT'}, parf1, data, 2, 10000);
else
[parf{itry}, normT(itry)] = MO.fitmodel_unw([parf1.N, parf1.tauD1, pari.thetaT+(1-2*rand)*pari.thetaT/8, pari.tauT+(1-2*rand)*pari.tauT/8], {'N', 'tauD1', 'thetaT', 'tauT'}, parf1, data, 2, 10000);
end
end
[val, pos] = min(normT);
parf2 = parf{pos};
outPar(1,:) = struct2array(parf2);
norm(1,1) = normT(pos);
end
%% fit thetaTtauT for a fixed kappa value for an array of files
function [thetaT, tauT] = fitThetaTtauTArray(MO, dataC, kappa, weight)
if ~iscell(dataC)
dataC = {dataC};
warning('fitThetaTtauTArray: data is not a cell array');
end
outParCorr = [];
outParCorr_thetaT = [];
for idata = 1:length(dataC)
data = dataC{idata}(:,[1 2 3]);
[thetTfit(idata,:), normCorr_thetaT(:,idata)] = MO.fitThetaTtauT(data, kappa, weight);
end
tauT = mean(thetTfit(:,3));
thetaT = mean(thetTfit(:,2));
end
%%
function [outPar, norm] = fitKappaVal(MO, data, kappaVal, weight)
if nargin < 4
weight = 1;
end
for i = 1:length(kappaVal)
pari = MO.par;
pari.kappa = kappaVal(i);
clear('parf')
for itry = 1:5
[N, tauD1, tauD2] = MO.initialGuess(data, 2, 10000);
if weight
[parf{itry}, normT(itry)] = MO.fitmodel([N, tauD1], {'N', 'tauD1'}, pari, data, 2, 10000);
else
[parf{itry}, normT(itry)] = MO.fitmodel_unw([N, tauD1], {'N', 'tauD1'}, pari, data, 2, 10000);
end
end
[val, pos] = min(normT);
parf1 = parf{pos};
if weight
[parf, norm(i,1)] = MO.fitmodel([parf1.N, parf1.tauD1], {'N', 'tauD1'}, parf1, data, 2, 10000);
else
[parf, norm(i,1)] = MO.fitmodel_unw([parf1.N, parf1.tauD1], {'N', 'tauD1'}, parf1, data, 2, 10000);
end
outPar(i,:) = struct2array(parf);
end
end
function [outPar, norm] = fitKappaValThetaT(MO, data, kappaVal, weight)
if nargin < 4
weight = 1;
end
%fit data varying kappaVal
for i = 1:length(kappaVal)
pari = MO.par;
pari.kappa = kappaVal(i);
for itry = 1:5
[N, tauD1, tauD2] = MO.initialGuess(data, 2, 10000);
if weight
[parf{itry}, normT(itry)] = MO.fitmodel([N, tauD1], {'N', 'tauD1'}, pari, data, 2, 10000);
else
[parf{itry}, normT(itry)] = MO.fitmodel_unw([N, tauD1], {'N', 'tauD1'}, pari, data, 2, 10000);
end
end
[val, pos] = min(normT);
parf1 = parf{pos};
%just fit thetaT
for itry = 1:5
if weight
[parf{itry}, normT(itry)] = MO.fitmodel([parf1.N, parf1.tauD1, pari.thetaT+(1-2*rand)*pari.thetaT/10 ], {'N', 'tauD1', 'thetaT'}, parf1, data, 2, 10000);
else
[parf{itry}, normT(itry)] = MO.fitmodel_unw([parf1.N, parf1.tauD1, pari.thetaT+(1-2*rand)*pari.thetaT/10 ], {'N', 'tauD1', 'thetaT'}, parf1, data, 2, 10000);
end
end
[val, pos] = min(normT);
parf2 = parf{pos};
norm(i,1) = normT(pos);
outPar(i,:) = struct2array(parf2);
end
end
end
end
function [G] = dye_diffusion_triplet(par,tau)
%Global Fitting for dye
no_traces = (size(par,2)-5)./4;
for i = 1:no_traces
%Defaults
N = par(6+4*(i-1));
Theta_T = par(7+4*(i-1));
Tau_T = par(8+4*(i-1));
f1 = par(9+4*(i-1));
Tau_D1 = par(1);
Alpha1 = par(2);
Kappa = par(5);
Tau_D2 = par(3);
Alpha2 = par(4);
G(:,i) = (1- Theta_T.*(1-exp(-tau./Tau_T))).*(f1.*((1+(tau./Tau_D1).^Alpha1).^-1)...
.*((1+((tau./Tau_D1).^Alpha1)./(Kappa^2)).^-0.5) + (1-f1)...
.*((1+(tau./Tau_D2).^Alpha2).^-1).*((1+((tau./Tau_D2).^Alpha2)./(Kappa^2)).^-0.5))./(N.*(1-Theta_T));
end
\ No newline at end of file
function [G] = dye_diffusion_triplet2(par,tau)
%Global Fitting for dye
no_traces = (size(par,2)-5)./3;
for i = 1:no_traces
%Defaults
N = par(6);
Theta_T = par(7+3*(i-1));
Tau_T = par(8+3*(i-1));
f1 = par(9+3*(i-1));
Tau_D1 = par(1);
Alpha1 = par(2);
Kappa = par(5);
Tau_D2 = par(3);
Alpha2 = par(4);
G(:,i) = (1- Theta_T.*(1-exp(-tau./Tau_T))).*(f1.*((1+(tau./Tau_D1).^Alpha1).^-1)...
.*((1+((tau./Tau_D1).^Alpha1)./(Kappa^2)).^-0.5) + (1-f1)...
.*((1+(tau./Tau_D2).^Alpha2).^-1).*((1+((tau./Tau_D2).^Alpha2)./(Kappa^2)).^-0.5))./(N.*(1-Theta_T));
end
\ No newline at end of file
function [fit_par,resnorm,par_table] = fcs_fitting(traces,method,Kappa,file_list)
%Defaults
N = 10;
Theta_T = 0.2;
Tau_T = 100;
f1 = 1;
Tau_D1 = 500;
Alpha1 = 1;
%Kappa = 5.5;
Tau_D2 = 5000;
Alpha2 = 1;
Tau_gfp = 271;