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

matlab files moved to FCSAnalyze

parent 86ad3d09
%%
% 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
tauBoundary;
model;
end
properties (Access = public)
weightChi2 = 1; % 1, use weighted Chi2; 0 use unweighted Chi2
maxfittry = 5; % number of repetitions for fit
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 [idxs] = getDataIndexes(MO, data, tauBoundary)
% GETDATAINDEXES get indexes within tauBoundary
% data - is [tau, corr] array
% tauBoundary - is [], [tauStart], [tauStart, tauEnd]
% idxs - are indexes of where to start (base 1) and
% where to end. idxs = [1, 0] use whole data set. Compatible
% with FA notation
tauStart = data(1,1);
tauEnd = data(end,1);
if nargin == 3
if ~isempty(tauBoundary)
tauStart = tauBoundary(1);
if length(tauBoundary) == 2
tauEnd = tauBoundary(2);
end
end
end
ist= find(data(:,1) >= tauStart); % index of start time
ien = find(data(:,1) >= tauEnd); % index of end time
idxs(1) = ist(1);
idxs(2) = size(data,1) - ien(1);
end
function [N, tauD1, tauD2] = initialGuess(MO, data, tauBoundary)
% 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 [tau, corr] 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
if nargin < 3
tauBoundary = MO.tauBoundary;
end
idxboundary = MO.getDataIndexes(data, tauBoundary);
data = data(idxboundary(1):end-idxboundary(2),:);
% only positive values for N and tauDm
idx = find(data(:,2)>0);
data = data(idx, :);
N = 1/mean(data(1:5,2));
% compute integral of equation to obtain a gues of average time
% only positive values
tauDm = N*trapz(data(:,1), data(:,2))/4.5;
%randomize
N = N + (1-2*rand)*N/8;
% 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, weight)
% GDIFFCOMP computes weighted difference of ACF data and model
% par - full parameter struct of model data - Ntx3 matrix of data [time, ACF, ACF_std]
if nargin < 4
weight = MO.weightChi2;
end
if weight
Gdiff = (MO.Gcomp(par, data(:,1)) - data(:,2))./data(:,3);
else
Gdiff = (MO.Gcomp(par, data(:,1)) - data(:,2));
end
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, weight)
% 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
if nargin < 6
weight = MO.weightChi2;
end
Gdiff = MO.Gdiffcomp(par, data, weight);
end
function [par, norm] = getbestfit(MO, parCA,normV)
% GETBESTFIT returns parameter and norm of best fit from
% several fits
% parCA - a cell array
% normV - a vector containing the norm
[val, pos] = min(normV);
norm = normV(pos);
par = parCA{pos};
end
function [R2, SSres] = compR2(MO, data, par, tauBoundary, weight)
%%
% COMPR2 computes R2 given the parameters
% Note that we are fitting a non-linear model so this has not
% much meaning for comparing models! It is just used to asses
% curves that are kind of flat.
% R2adj is used in FA
if nargin < 4
tauBoundary = MO.tauBoundary;
end
if nargin < 5
weight = MO.weightChi2;
end
idxboundary = MO.getDataIndexes(data, tauBoundary);
data = data(idxboundary(1):end-idxboundary(2),:);
if isstruct(par)
SSres = sum(MO.Gdiffcomp(par, data, weight).^2);
else
for i=1:size(par,1)
parS = MO.parvecToparstruct(par(i,:));
SSres(i,1) = sum(MO.Gdiffcomp(parS, data, weight).^2);
end
end
if weight
SStot = sum((data(:,2)./data(:,3) - 1./data(:,3)*mean(data(:,2))).^2);
else
SStot = sum((data(:,2) - mean(data(:,2))).^2);
end
R2 = 1- SSres./SStot;
% R2adj is corrected for the number of parameters that have
% been fitted
% R2adj_2c = 1 - (1-R2_2c)*(Nrpts2c-1)/(Nrpts2c - nrpara2c - 1) (https://en.wikipedia.org/wiki/Coefficient_of_determination#Adjusted_R2)
% We don't pass the Npar fitted to the function here
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(MO, parf, names, par, data, tauBoundary)
% FITMODEL fit model to data depending weightChi2 this use weighted difference or not
% [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, tauBoundary) - fit in time intervael
% [tauBoundary(1), tauMax] or [tauBoundary(1) tauBoundary(2)]
% if tauBoundary has length 2
% 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
% norm - sum of squared residuals at minima
if nargin < 6
tauBoundary = MO.tauBoundary;
end
idxboundary = MO.getDataIndexes(data, tauBoundary);
data = data(idxboundary(1):end-idxboundary(2),:);
lb = [];
hb = [];
for i = 1:length(names)
lb = [lb;getfield(MO.lb,names{i})];
hb = [hb;getfield(MO.hb,names{i})];
end
options = optimset('display','off');
[parf, norm] = lsqnonlin(@MO.Gdiffcompfit, parf, lb, hb, options, names, par, data);
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
% reads all entries without differentiating between Ch1 and Ch2.
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');
fclose(fileID);
%%
% find repetion position and channel
RPC = [];
positionIdx = find(strncmp(C{1}, 'Position = ', 11));
repetitionIdx = find(strncmp(C{1}, 'Repetition = ', 13));
channelIdx = find(strncmp(C{1}, 'Channel = ', 10));
for i=1:length( positionIdx)-1
tmp = C{1}(positionIdx(i));
[tokp] = regexp(tmp{:},'(\d+)', 'tokens');
tmp = C{1}(repetitionIdx(i));
[tokr] = regexp(tmp{:},'(\d+)', 'tokens');
tmp = C{1}(channelIdx(i));
[tokc] = regexp(tmp{:},'(\d+)', 'tokens');
% invert channel for compability with FA this assumes APD!!!
% should check for ChS1 and ChS2
if str2num(char(tokc{:}))==2
ch = 1;
elseif str2num(char(tokc{:}))==1
ch = 2;
end
RPC = [RPC; str2num(char(tokr{:}))+1 str2num(char(tokp{:}))+1 ch ];
end
% index for measurements is based on repetition or measured points
idxR = 2;
if max(RPC(:,1))>1
idxR = 1;
end
if max(RPC(:,2))>1
idxR = 2;
end
%%
%488_561nm_52f9463241b00a79480ea084e4155cbc_R1_P1_K1_Ch2.raw
ACm = find(strncmp(C{1}, 'CorrelationArray =',18 )); %this reads all entries.
AC = cell(1,max(max(RPC(:,1:2))));
%last entry is fitting etc.
for i=1:length(ACm)-1
dim = str2num(C{1}{ACm(i)}(19:end));
ACl = cell2mat(cellfun(@str2num, C{1}(ACm(i)+1:ACm(i)+dim(1)), 'un', 0));
idx = RPC(i, idxR);
% convert to us and AC base 0
if isempty(AC{idx})
AC{idx} = zeros(size(ACl,1), 7);
end
AC{idx}(:,1) = ACl(:,1)*1e6;
AC{idx}(:,RPC(i,3)*2:RPC(i,3)*2+1) = [ACl(:,2)-1 ones(size(ACl,1),1)];
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);
% add columns as required to have all 6 entries Ch1, Ch2, XC
if size(AC{i},2) < 7
AC{i}(:,size(AC{i},2):7) = 0;
end
end
end
function AC = readFAxml(MO, fnames)
% READFAXML read AC, parameters, annotations etc from xml file generated per measurement by FA
if nargin < 2
%example file
path = fileparts(fileparts(mfilename('fullpath')));
fnames = {fullfile(path, 'example_data', '488_561nm_R1_P1_K1_Ch2-Ch1.zen.1c.xml')}
end
for i=1:length(fnames)
AC{i} = xmlFA.readcor(fnames{i});
end
end
end
end
classdef dyeFcsmodel < absFcsmodel
% DYEFCSMODEL class to fit auto correlation function coming from dye measurement
% a one component model specific for dye like blinking can be fitted
% There are specific fitting workflows implemented.
% fit is performed by
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', 0, '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);
diffCoeff = [464, 521]; %Alexe488 at 37 degrees, Alexe568 at 37 degrees Malte Wachsmuth SPIM FCS measurement
tauBoundary = [2, 10000]; %boundary of where to fit the data in us
model = 1;
end
methods
function MO = dyeFcsmodel(parnames, parvalues, lbvalues, hbvalues)
% DYEFCSMODEL constructor
% MO = dyeFcsmodel() - just create class with default parameter values and boundary values
% MO = dyeFcsmodel(parnames, parvalues) - create class and assign parvalues to parameters named parnames
% MO = dyeFcsmodel(parnames, parvalues, lbvalues)- create class and assign parvalues, lower boundary to parameter named in parnames.
% MO = dyeFcsmodel(parnames, parvalues, lbvalues,hbvalues) - create class and assign parvalues, lower and higher boundary to parameter named in parnames.
% parnames - cell array containing name of parameteres to be changed
% parvalues - vector with values of parameters same length than parnames
% lbvalues - vector with values of low boundary parameters same length than parnames
% hbvalues - vector with values of high boundary parameters same length than parnames
% we cannot call MO.par in the constructor before call of super class this is the reason for the repetitions of the par values
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
if ~isstruct(par)
par = MO.parvecToparstruct(par);
end
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 femtoliter
% par - vector or struct of parameters
if isstruct(par)
vol = pi^(3/2).*MO.computeRadius(par, dyetype).^3.*par.kappa;
else
vol = pi^(3/2).*MO.computeRadius(par, dyetype).^3.*par(:,9);
end
end
function vol = computeConcentration(MO, par, dyetype)
% COMPUTECONCENTRATION Compute concentration and returns value
% in nM
% par - vector or struct of parameters
% Na = 0.6022140 (=
% 6.022140*10^23*10^-15(femtoliter)*10^-9(to get nM))
Na = 0.6022140;
if isstruct(par)
vol = par.N./MO.computeVolume(par, dyetype)/Na;
else
vol = par(:,1)./MO.computeVolume(par, dyetype)/Na;
end
end
function w0 = computeRadius(MO, par, dyetype)
% COMPUTERADIUS Compute focal radius and return value in um
% par - vector or struct of parameters
if isstruct(par)
w0 = 2*sqrt(MO.diffCoeff(dyetype)*par.tauD1)/1000;
else
w0 = 2*sqrt(MO.diffCoeff(dyetype)*par(:,5))/1000;
end
end
function [outPar, norm] = fitThetaTtauT(MO, data, kappa)
% FITTHETATTAUT fit N, tauD1, thetaT and tauT for a fixed kappa value
% INPUT:
% data - correlation data [tau, corr]
% kappa - kappa value will be kept fixed. Typically at 5 (this is more stable than allowing for verying thetaT and kappa simultaneously)
% OUTPUT:
% outPar - 'N', 'tauD1', 'thetaT', 'tauT'
% norm - norm at optimum
if nargin < 3
kappa = 5;
end
%first find optimal tauT
pari = MO.par;
pari.kappa = kappa;
%find initial value for tauD1 and N repeat 5 times
for itry = 1:MO.maxfittry
[N, tauD1, tauD2] = MO.initialGuess(data, MO.tauBoundary);
[parf{itry}, normT(itry)] = MO.fitmodel([N, tauD1], {'N', 'tauD1'}, pari, data, MO.tauBoundary);
end
[val, pos] = min(normT);
parf1 = parf{pos};
%fit thetaT and tauT for a fixed kappa
for itry = 1:MO.maxfittry
[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, MO.tauBoundary);
end
[val, pos] = min(normT);
parf2 = parf{pos};
outPar(1,:) = struct2array(parf2);
norm(1,1) = normT(pos);
end
function [thetaT, tauT] = fitThetaTtauTArray(MO, dataC, kappa, iC)
% FITTHETATTAUTARRAY find optimal thetaT and tauT for several data simultaneously with fixed kappa
% INPUT:
% dataC - a cell array containing [tau corrCh1 corrCh1_e corrCh2 corrCh2_e Xcorr Xcorr_e] arrays
% kappa - fixed kappa value
% iC - index of channel to be used 1 - Ch1, 2 - Ch2, 3 - XCorr
% OUTPUT:
% thetaT - mean optimal thetaT for several measurements
% tauT - mean optimal tauT for several measurements
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 iC*2 iC*2+1]);
[thetTfit(idata,:), normCorr_thetaT(:,idata)] = MO.fitThetaTtauT(data, kappa);
end
tauT = mean(thetTfit(:,3));
thetaT = mean(thetTfit(:,2));
end
function [outPar, norm] = fitKappaVal(MO, data, kappaVal)
% FITKAPPAVAL compute N and tauD1 for different values of kappa. All other parameters remain constant
% data - an array [tau, corr]
% kappaVal - a vector containing values for kappa
for i = 1:length(kappaVal)
pari = MO.par;
pari.kappa = kappaVal(i);
clear('parf')
for itry = 1:MO.maxfittry
[N, tauD1, tauD2] = MO.initialGuess(data, 2, 10000);
[parf{itry}, normT(itry)] = MO.fitmodel([N, tauD1], {'N', 'tauD1'}, pari, data, MO.tauBoundary);
end
[val, pos] = min(normT);
parf1 = parf{pos};
[parf, norm(i,1)] = MO.fitmodel([parf1.N, parf1.tauD1], {'N', 'tauD1'}, parf1, data, MO.tauBoundary);
outPar(i,:) = struct2array(parf);
end
end
function [outPar, norm] = fitKappaValThetaT(MO, data, kappaVal)
% FITKAPPAVALTHETAT compute N and tauD1 for different values of kappa. All other parameters remain constant
% INPUT:
% data - an array [tau, corr]
% kappaVal - a vector containing values for kappa
% OUTPUT:
% outPar - N, tauD1, thetaT
% norm - norm of fit
for i = 1:length(kappaVal)
pari = MO.par;
pari.kappa = kappaVal(i);
%display(sprintf('Run fit for kappa %.1f ',kappaVal(i)));
% get initial guess for N and tauD1
for itry = 1:MO.maxfittry
[N, tauD1, tauD2] = MO.initialGuess(data, MO.tauBoundary);
[parf{itry}, normT(itry)] = MO.fitmodel([N, tauD1], {'N', 'tauD1'}, pari, data, MO.tauBoundary);
end
[val, pos] = min(normT);
parf1 = parf{pos};
%fit N, tauD1 and thetaT using initial guess of N and tauD1, tauR remains fixed for all measurements
for itry = 1:MO.maxfittry
[parf{itry}, normT(itry)] = MO.fitmodel([parf1.N, parf1.tauD1, pari.thetaT+(1-2*rand)*pari.thetaT/10 ], ...
{'N', 'tauD1', 'thetaT'}, parf1, data, MO.tauBoundary);
end
[val, pos] = min(normT);
parf2 = parf{pos};
norm(i,1) = normT(pos);
outPar(i,:) = struct2array(parf2);
end
end
end
end
This diff is collapsed.