Skip to content

GitLab

Commit 7deb3379 authored by Antonio Politi's avatar Antonio Politi
Browse files

complete dyeFit pipeline

parent 88a0eb85
Hide whitespace changes
Inline Side-by-side
matlabcode/FCS/masterFitDye.m
View file @ 7deb3379
function masterFitDye(indirCell)
% MASTERFITDYE perform a fitting procedure for several directories containig dye data
% workflow also finds best kappa.
% workflow is first find best tauT for all measurments of one experiment at a given kappa
% the starting value does not matter much, a value of 5 is a good compromise.
% Then fix tauT and fit N, tauD1, thetaT for different values of kappa
% repeat this procedure several times
% workflow also finds best kappa (ratio of z to xy focal volume).
% indirCell - A cell array containing path to data for dye
% Data is outputted in working directories
% Workflow is:
% 1. Find best tauT for all dye measurements of one experiment at a given kappa
% 2. Fix tauT and fit N, tauD1, and thetaT for different values of kappa.
% 3. Find kappa that yields the lowest mean Chi2 (each Chi2 is normalized)
% 4. Repeat 1-3 with new kappa (this is repeated ~ 3 times)
%
% Antonio Politi, EMBL, January 2017
if nargin < 1
% indirCell = {'W:\Shotaro\Live-cell-imaging\HeLa-FCS\160701-HeLa-Nup107\160701_gfpNUP107z26z31\Calibration\Alexa', ...
......@@ -20,40 +23,59 @@ if nargin < 1
indirCell = {'C:\Users\toni\Dropbox\NPCMaturation\matlabcode\FCS\example_data\'};
end
end
% output subdiretoy to store results
outputdir = 'optimisedFit'
% range to fit kappa
kappaVal = [2:0.1:7]';
for idir = 1:length(indirCell)
indir = indirCell{idir};
computeDirectory(indir, kappaVal)
computeDirectory(indir, kappaVal, outputdir)
end
function computeDirectory(indir, kappaVal)
function computeDirectory(indir, kappaVal, outputdir)
% COMPUTEDIRECTORY perform fitting on a specific directory
% inidir - directory where to extract files
% kappaVal - a vector of kappa values
% outputdir - subdirectory name to store some results of the fit
% Workflow is:
% 1. Find best tauT for all dye measurements of one experiment at a given kappa
% 2. Fix tauT and fit N, tauD1, and thetaT for different values of kappa.
% 3. Find kappa that yields the lowest mean Chi2 (each Chi2 is normalized)
% 4. Repeat 1-3 with new kappa (this is repeated ~ 3 times)
%
% several conditions to run the fit
%% Initialization
% load model class
clear('MO');
MO = dyeFcsmodel();
outdir = fullfile(indir,'optimisedFit');
% make outputdir if required
outdir = fullfile(indir, outputdir);
mkdir(outdir);
zenfiles = getAllFiles(indir, 'fcs');
% find files from ZEN and FA
zenfiles = getAllFiles(indir, 'fcs');
xmlFAfiles = getAllFiles(indir, 'xml', '1c');
corFiles = getAllFiles(indir,'cor', '1c');
%corFiles = getAllFiles(indir,'cor', '1c'); this is obsolete as we do
%extract data from xml
%%
for cond = 1
switch cond
case 1
case 1 %fit using data from FA weighted with standard deviation
MO.weightChi2 = 1;
ZEN = 0;
case 2
case 2 %fit using data from FA unweighted
MO.weightChi2 = 0;
ZEN = 0;
case 3
case 3 %fit using data from .fcs ZEN files
MO.weightChi2 = 0; % this is unweighted by definition as std is not available for .fcs
ZEN = 1;
end
% define file prefix
if ZEN
prefix = 'parZEN';
else
......@@ -63,105 +85,193 @@ for cond = 1
prefix = [prefix '_unw'];
end
% load data into cell array dataC
if ZEN
%%
tic
dataC = MO.readZeissfcs(zenfiles{1});
toc
else
%%
dataC = MO.readFAxml(xmlFAfiles);
%dataC2 = MO.readFAcor(corFiles);
end
norm = [];
%% find optimal tauT at a constant kappa value here = 5
kappa = 5;
MO.maxfittry = 5;
for irounds = 1:3
display(['Find optimal tauT and thetaT round ' irounds] )
[thetaT, tauT] = MO.fitThetaTtauTArray(dataC, kappa);
% associate these values to default parameter set
MO.par.tauT = tauT;
MO.par.thetaT = thetaT;
%%
outParCorr = [];
normCorr = [];
display(['Fit different kappa values ' irounds])
%this fit keep tauT constant for all measurments and fit to each measurement a thetaT
for idata = 1:length(dataC)
data = dataC{idata}(:,[1 2 3]);
[outParl, normCorr] = MO.fitKappaValThetaT(data, kappaVal);
outParCorr = [outParCorr; outParl, normCorr];
%% Workflow
lk = length(kappaVal);
if exist(fullfile(outdir, [prefix '.mat'])) %load result file if it exists already
load(fullfile(outdir, [prefix '.mat']));
else
% starting value of kappa
kappa = 5;
% number of repeats for each fit (starts from random values)
MO.maxfittry = 1;
maxRounds = 1;
% struct to store data for each channel
Chfits = struct('outPar', [], 'kappaVal', kappaVal, 'kappa', [],'posBest', [], 'focRad', [], 'mfocRad', [], 'focVol', [], ...
'mfocVol', [], 'conc', [], 'mConc', [], 'normN', [], 'mNorm', []);
Chfits(2) = Chfits;
% for each channel
for iC = 1:2
if (dataC{1}(1,iC*2) == 0) % do not process if channel has not been acquired
continue
end
for irounds = 1:maxRounds % repeat workflow 3 times
%% find optimal tauT and thetaT
display(sprintf('Find optimal tauT and thetaT for Channel %d round %d/3', iC, irounds));
%default starting values
MO.par.tauT = 7;
MO.par.thetaT = 0.2;
[thetaT, tauT] = MO.fitThetaTtauTArray(dataC, kappa, iC);
% associate these values to default parameter set
MO.par.tauT = tauT;
MO.par.thetaT = thetaT;
%% fit N, tauD1, thetaT for constant tauT and different kappa values
display(sprintf('Fit different kappa values for Channel %d round %d/3', iC, irounds));
% [N_kappaVal*Nexp Npar] array of parameters (all) for each kappa
% value and experiment. Last value is norm
outPar = [];
for idata = 1:length(dataC)
data = dataC{idata}(:,[1 iC*2 iC*2+1]);
[outParl, normCorr] = MO.fitKappaValThetaT(data, kappaVal);
outPar = [outPar; outParl, normCorr];
end
%% compute focal radius and volume
focRad = []; % focal Radius
focVol = []; % focal Volume
normN = []; % normalized Norm for maximal value for a given experiment
for i = 1:length(dataC)
normN(:,i) = outPar(1+(i-1)*lk:lk*i,end)/max(outPar(1+(i-1)*lk:lk*i,end));
focVol(:,i) = MO.computeVolume(outPar(1+(i-1)*lk:lk*i,1:end-1),iC);
focRad(:,i) = MO.computeRadius(outPar(1+(i-1)*lk:lk*i,1:end-1),iC);
conc(:,i) = MO.computeConcentration(outPar(1+(i-1)*lk:lk*i,1:end-1), iC);
end
%% find best fit
mNorm = mean(normN,2); % mean of normalized norms along experiments
[val, pos] = min(mean(normN,2)) % find minimum at pos
kappa = kappaVal(pos); % kappa value
end
%% Compute means from mean tauD1 and mean N. Using directly the data could increase error on focal volume (^3)
for i = 1:lk
% find indexes for given kappaValue
posExp = find(round(10*outPar(:,9))==round(10*kappaVal(i)));
% compute mean focal volume (pl) for all experiments at given
% kappa
mVol(i) = MO.computeVolume(mean(outPar(posExp, 1:end-1),1),iC);
% compute mean focal radius (um) for all experiments at given
% kappa
mw0(i) = MO.computeRadius(mean(outPar(posExp, 1:end-1),1),iC);
% compute mean concentration(nM) for all experiments at given
% kappa
mConc(i) = MO.computeConcentration(mean(outPar(posExp, 1:end-1),1),iC);
end
% store in struct
Chfits(iC) = struct('outPar', outPar, 'kappaVal', kappaVal, 'kappa', kappa, 'posBest', pos, 'focRad', focRad, 'mfocRad', mw0, 'focVol', focVol, ...
'mfocVol', mVol, 'conc', conc, 'mConc', mConc, 'normN', normN, 'mNorm', mNorm);
end
radCorr = [];
volCorr = [];
normCorrN = [];
lk = length(kappaVal);
for i = 1:length(dataC)
normCorrN(:,i) = outParCorr(1+(i-1)*lk:lk*i,end)/max(outParCorr(1+(i-1)*lk:lk*i,end));
volCorr(:,i) = MO.computeVolume(outParCorr(1+(i-1)*lk:lk*i,1:end-1),1);
radCorr(:,i) = MO.computeRadius(outParCorr(1+(i-1)*lk:lk*i,1:end-1),1);
%% save results
% duplicate entry in case only one Channel has been recorded
if (dataC{1}(1,2) == 0) % do not process if channel has not been acquired
Chfits(1) = Chfits(2);
end
mNorm = mean(normCorrN,2);
[val, pos] = min(mean(normCorrN,2))
kappa = kappaVal(pos);
if (dataC{1}(1,4) == 0)
Chfits(2) = Chfits(1);
end
save(fullfile(outdir, [prefix '.mat']), 'Chfits');
end
%%
% save focal volume means use mean of 2 channel if available otherwise
% what is available
posB = Chfits(1).posBest;
focV = [Chfits(1).mfocRad(posB) Chfits(1).mfocVol(posB) Chfits(1).kappaVal(posB);...
Chfits(2).mfocRad(posB) Chfits(2).mfocVol(posB) Chfits(2).kappaVal(posB); ...
sqrt(Chfits(1).mfocRad(posB)*Chfits(2).mfocRad(posB)) sqrt(Chfits(1).mfocVol(posB)*Chfits(2).mfocVol(posB)) Chfits(2).kappaVal(posB)]
posB = Chfits(2).posBest;
focV = [Chfits(1).mfocRad(posB) Chfits(1).mfocVol(posB) Chfits(1).kappaVal(posB);...
Chfits(2).mfocRad(posB) Chfits(2).mfocVol(posB) Chfits(2).kappaVal(posB); ...
sqrt(Chfits(1).mfocRad(posB)*Chfits(2).mfocRad(posB)) sqrt(Chfits(1).mfocVol(posB)*Chfits(2).mfocVol(posB)) Chfits(2).kappaVal(posB)]
posB = round(mean([Chfits(1).posBest,Chfits(2).posBest]));
focV = [Chfits(1).mfocRad(posB) Chfits(1).mfocVol(posB) Chfits(1).kappaVal(posB);...
Chfits(2).mfocRad(posB) Chfits(2).mfocVol(posB) Chfits(2).kappaVal(posB); ...
sqrt(Chfits(1).mfocRad(posB)*Chfits(2).mfocRad(posB)) sqrt(Chfits(1).mfocVol(posB)*Chfits(2).mfocVol(posB)) Chfits(2).kappaVal(posB)]
fid = fopen(fullfile(indir, 'Focus.txt'), 'w');
fprintf(fid,'w0_um\tVol_fl\tkappa\r\n');
for i = 1:3
fprintf(fid,'%.4f\t%.4f\t%.1f\r\n', focV(i,1),focV(i,2),focV(i,3) );
end
fclose(fid);
%save(fullfile(outdir,[prefix '.txt']), 'outParCorr', '-ascii')
% %%
%%
% mean is computed from mean tauD1
for i = 1:lk
posExp = find(round(10*outParCorr(:,9))==round(10*kappaVal(i)));
mVol(i) = MO.computeVolume(mean(outParCorr(posExp, 1:end-1),1),1);
mw0(i) = MO.computeRadius(mean(outParCorr(posExp, 1:end-1),1),1);
end
%%
%
hf = figure(1)
clf
hold
subplot(3,1,1)
hold
errorbar(kappaVal, mean(normCorrN,2), std(normCorrN, 1,2), '-')
[val, pos] = min(mean(normCorrN,2))
kappa = kappaVal(pos);
plot(kappa, mNorm(pos), 'ro', 'MarkerFaceColor', 'r' )
ylim([0.7 1])
ylabel('normalized Chi2');
xlabel('kappa')
title(['kappa ' num2str(kappaVal(pos)) '; Vol ' num2str(mVol(pos)) ' pl; w0 ' num2str(mw0(pos)) ' um'])
subplot(3,1,2)
hold
errorbar(kappaVal, mean(volCorr,2),std(volCorr, 1,2), '-')
plot(kappa, mVol(pos), 'ro', 'MarkerFaceColor', 'r' )
ylabel('Confocal volume (pl)')
xlabel('kappa')
%% plot results
for iC = 1:2
if (dataC{1}(1,iC*2) == 0) % do not process if channel has not been acquired
continue
end
% means are computed from mean tauD1 and mean N. Otherwise error
% increases
for i = 1:lk
% find indexes for given kappaValue
posExp = find(round(10*Chfits(iC).outPar(:,9))==round(10*kappaVal(i)));
% compute mean focal volume (pl) for all experiments at given
% kappa
mVol(i) = MO.computeVolume(mean(Chfits(iC).outPar(posExp, 1:end-1),1),iC);
% compute mean focal radius (um) for all experiments at given
% kappa
mw0(i) = MO.computeRadius(mean(Chfits(iC).outPar(posExp, 1:end-1),1),iC);
% compute mean concentration(nM) for all experiments at given
% kappa
mConc(i) = MO.computeConcentration(mean(Chfits(iC).outPar(posExp, 1:end-1),1),iC);
end
hf = setupFigure(iC, [100 100 400 800] )
clf
% normalized Chi2
subplot(4,1,1)
hold
posB = Chfits(iC).posBest;
errorbar(kappaVal, mean(Chfits(iC).normN,2), std(Chfits(iC).normN, 1,2), '-')
plot(Chfits(iC).kappa, Chfits(iC).mNorm(posB), 'ro', 'MarkerFaceColor', 'r' )
ylim([min(Chfits(iC).normN(:)) max(Chfits(iC).normN(:))])
ylabel('normalized Chi2');
xlabel('kappa')
title(sprintf('Channel %d kappa %.1f; Vol %.3f pl; w0 %.3f um', iC, Chfits(iC).kappa, Chfits(iC).mfocVol(posB), Chfits(iC).mfocRad(posB)))
% focal volume
subplot(4,1,2)
hold
errorbar(kappaVal, mean(Chfits(iC).focVol,2),std(Chfits(iC).focVol, 1,2), '-')
plot(Chfits(iC).kappa, Chfits(iC).mfocVol(posB), 'ro', 'MarkerFaceColor', 'r' )
ylabel('Confocal volume (pl)')
xlabel('kappa')
ylim([min(Chfits(iC).focVol(:)) max(Chfits(iC).focVol(:))])
% focal radius
subplot(4,1,3)
hold
errorbar(kappaVal, mean(Chfits(iC).focRad,2), std(Chfits(iC).focRad,1,2), '-')
plot(Chfits(iC).kappa,Chfits(iC).mfocRad(posB), 'ro', 'MarkerFaceColor', 'r' )
ylabel('Confocal radius w0 (um)')
xlabel('kappa')
ylim([min(Chfits(iC).focRad(:)) max(Chfits(iC).focRad(:))])
% concentration
subplot(4,1,4)
hold
errorbar(kappaVal,mean(Chfits(iC).conc,2), std(Chfits(iC).conc,1,2), '-')
plot(Chfits(iC).kappa, Chfits(iC).mConc(posB), 'ro', 'MarkerFaceColor', 'r' )
ylabel('Concentration (nM)')
xlabel('kappa')
saveas(hf,fullfile(outdir,[prefix '_Ch' num2str(iC) '.png']), 'png')
end
ylim([0.1 0.3])
subplot(3,1,3)
hold
errorbar(kappaVal, mean(radCorr,2), std(radCorr,1,2),std(radCorr,1,2), '-')
plot(kappa, mw0(pos), 'ro', 'MarkerFaceColor', 'r' )
ylabel('Confocal radius w0 (um)')
xlabel('kappa')
ylim([0.15 0.3])
saveas(hf,fullfile(outdir,[prefix '.png']), 'png')
end
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment