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

fix to properly process CC data

parent 626768d0
function master_workflow_fit_dye_subdirectories(wd)
% MASTER_WORKFLOW_FIT_DYE_SUBDIRECTORIES call workflow_fit_dye for all
% subdirectories in a main directory (wd)
% The program search only for one level down wd. Subdirectories should contain correlation traces generated by FA
% MASTER_WORKFLOW_FIT_DYE_SUBDIRECTORIES() - prompt for a main directory
% MASTER_WORKFLOW_FIT_DYE_SUBDIRECTORIES(wd) - use wd as main directory
%
% Antonio Politi, EMBL, April 2017
if nargin < 1
wd = uigetdir('.', 'Specify main directory where to search for subfolders');
cd(wd);
end
% get all subdirectories
[indirCell] = getAllFolders(wd, '.*');
indirCell = {wd, indirCell{:}};
% output subdiretoy to store results
outputdir = 'optimisedFit';
% range to fit kappa
kappaVal = [2:0.1:9]';
force = 0;
alexa488 = [464, 464, 464];
alexa568 = [521, 521, 521];
atto488 = [509, 509, 509];
FAsession = '1c';
for idir = 1:length(indirCell)
indir = indirCell{idir};
% check for xml files in this subdirectory otherwise move to next
% directory
if isempty(getAllFiles(indir, 'xml',['\.' FAsession '\.'], 0))
continue
end
display(['master_workflow_fit_dye processing directory ' num2str(idir) ' ' indir])
try
workflow_fit_dye(indir, kappaVal, outputdir, FAsession, force, alexa488);
catch ME
display(['failed master_workflow_fit_dye processing directory ' num2str(idir) ' ' indir ' ']);
display(getReport(ME))
end
end
\ No newline at end of file
...@@ -7,7 +7,6 @@ session = '2c' ...@@ -7,7 +7,6 @@ session = '2c'
force = 1; force = 1;
if nargin <1 if nargin <1
wd = uigetdir('.', 'Specify main directory where to search protein folders'); wd = uigetdir('.', 'Specify main directory where to search protein folders');
%wd = 'Z:\AntonioP_elltier1\CelllinesImaging\MitoSysPipelines'; %wd = 'Z:\AntonioP_elltier1\CelllinesImaging\MitoSysPipelines';
cd(wd); cd(wd);
end end
...@@ -16,7 +15,7 @@ if ~exist(fullfile(wd, 'ProcessingHelpFiles')); ...@@ -16,7 +15,7 @@ if ~exist(fullfile(wd, 'ProcessingHelpFiles'));
end end
if ~exist(fullfile(wd, 'ProcessingHelpFiles', 'protein2cfiles.mat')) || force if ~exist(fullfile(wd, 'ProcessingHelpFiles', 'protein2cfiles.mat')) || force
[indirCell] = getAllFolders(wd, '\Mito(s|S)ys\d'); [indirCell] = getAllFolders(wd, '(\Mito(s|S)ys\d)|(LSM880)');
resfiles = []; resfiles = [];
for i = 1:length(indirCell) for i = 1:length(indirCell)
resfiles = [resfiles; getAllFiles(indirCell{i}, 'res', [session '.res'], 0)]; resfiles = [resfiles; getAllFiles(indirCell{i}, 'res', [session '.res'], 0)];
......
function [focVolA, focVolfile] = workflow_fit_dye(indir, kappaVal, outputdir, FAsession, force) function [focVolA, focVolfile] = workflow_fit_dye(indir, kappaVal, outputdir, FAsession, force, diffCoeff)
% COMPUTEDIRECTORY perform fitting of dye in a specific directory. program % COMPUTEDIRECTORY perform fitting of dye in a specific directory. program
% assumes that traces of one experiment are all in one directory! % assumes that traces of one experiment are all in one directory!
% INPUT: % INPUT:
...@@ -9,18 +9,18 @@ function [focVolA, focVolfile] = workflow_fit_dye(indir, kappaVal, outputdir, FA ...@@ -9,18 +9,18 @@ function [focVolA, focVolfile] = workflow_fit_dye(indir, kappaVal, outputdir, FA
% OUTPUT: % OUTPUT:
% focVolA - [w0_Ch1 Vol_Ch1 kappa_Ch1; w0_Ch2 Vol_Ch2 kappa_Ch2; % focVolA - [w0_Ch1 Vol_Ch1 kappa_Ch1; w0_Ch2 Vol_Ch2 kappa_Ch2;
% w0_XC Vol_XC kappa_XC] if 2 channels have been acquired the mean of kappa values is used to compute focal volumes. % w0_XC Vol_XC kappa_XC] if 2 channels have been acquired the mean of kappa values is used to compute focal volumes.
% %
% focVolA is also saved in indir/focalVolume.txt % focVolA is also saved in indir/focalVolume.txt
% images and mat file are saved in outputdir % images and mat file are saved in outputdir
% %
% Workflow is: % Workflow is:
% 1. Find best tauT for all dye measurements of one experiment at a given kappa % 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. % 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) % 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) % 4. Repeat 1-3 with new kappa (this is repeated ~ 3 times)
% %
% Antonio Politi, EMBL, January 2017 % Antonio Politi, EMBL, January 2017
replot = 1;
%% default values %% default values
if nargin < 1 if nargin < 1
% prompt for a directory % prompt for a directory
...@@ -50,23 +50,31 @@ end ...@@ -50,23 +50,31 @@ end
clear('MO'); clear('MO');
MO = dyeFcsmodel(); MO = dyeFcsmodel();
% make outputdir if required if nargin == 6
outdir = fullfile(indir, outputdir); MO.diffCoeff = diffCoeff;
if ~exist(outdir)
mkdir(outdir);
end end
% find files from ZEN and FA % find files from ZEN and FA
zenfiles = getAllFiles(indir, 'fcs'); zenfiles = getAllFiles(indir, 'fcs');
xmlFAfiles = getAllFiles(indir, 'xml',['\.' FAsession '\.']); xmlFAfiles = getAllFiles(indir, 'xml',['\.' FAsession '\.']);
%corFiles = getAllFiles(indir,'cor', '1c'); this is obsolete as we do extract data from xml %corFiles = getAllFiles(indir,'cor', '1c'); this is obsolete as we do extract data from xml
%% %%
if isempty(xmlFAfiles) & isempty(zenfiles)
display('No FA xml or .fcs found stop here');
return
end
if isempty(xmlFAfiles) if isempty(xmlFAfiles)
display('No FA xml found use .fcs file instead'); display('No FA xml found use .fcs file instead');
cond = 3; cond = 3;
else else
cond = 1; cond = 1;
end end
% make outputdir if required
outdir = fullfile(indir, outputdir);
if ~exist(outdir)
mkdir(outdir);
end
switch cond switch cond
case 1 %fit using data from FA weighted with standard deviation case 1 %fit using data from FA weighted with standard deviation
...@@ -109,33 +117,44 @@ if exist(fullfile(outdir, [prefix '.mat'])) %load result file if it exists alrea ...@@ -109,33 +117,44 @@ if exist(fullfile(outdir, [prefix '.mat'])) %load result file if it exists alrea
process =1; process =1;
end end
end end
if (process || force)
%% load data into cell array dataC
if ZEN if ZEN
dataC = MO.readZeissfcs(zenfiles{1}); dataC = MO.readZeissfcs(zenfiles{1});
else else
dataC = MO.readFAxml(xmlFAfiles); dataC = MO.readFAxml(xmlFAfiles);
end end
if isempty(dataC)
return
end
if (process || force)
%% load data into cell array dataC
if ZEN
dataC = MO.readZeissfcs(zenfiles{1});
else
dataC = MO.readFAxml(xmlFAfiles);
end
if isempty(dataC)
return
end
% starting value of kappa % starting value of kappa
kappa = 5; kappa = 5;
% number of repeats for each fit (starts from random values) % number of repeats for each fit (starts from random values)
%MO.maxfittry = 1; % MO.maxfittry = 1;
%maxRounds = 1; % maxRounds = 1;
MO.maxfittry = 5; MO.maxfittry = 5;
maxRounds = 3; maxRounds = 3;
% struct to store data for each channel % struct to store data for each channel
Chfits = struct('outPar', [], 'kappaVal', kappaVal, 'kappa', [],'posBest', [], 'focRad', [], 'mfocRad', [], 'focVol', [], ... Chfits = struct('outPar', [], 'kappaVal', kappaVal, 'kappa', [],'posBest', [], 'focRad', [], 'mfocRad', [], 'focVol', [], ...
'mfocVol', [], 'conc', [], 'mConc', [], 'normN', [], 'mNorm', []); 'mfocVol', [], 'conc', [], 'mConc', [], 'normN', [], 'mNorm', []);
Chfits(2) = Chfits; Chfits(2) = Chfits(1);
Chfits(3) = Chfits(1);
% for each channel % for each channel
for iC = 1:2 for iC = 1:3
if (dataC{1}(1,iC*2) == 0) % do not process if channel has not been acquired if (dataC{1}(1,iC*2) == 0) % do not process if channel has not been acquired
continue continue
end end
for irounds = 1:maxRounds % repeat workflow 3 times with updated optimal kappavalue at each round for irounds = 1:maxRounds % repeat workflow 3 times with updated optimal kappavalue at each round
%% find optimal tauT and thetaT %% find optimal tauT and thetaT
display(sprintf('Alexa dye fit optimal tauT and thetaT for Channel %d round %d/%d', iC, irounds, maxRounds)); display(sprintf('Alexa dye fit optimal tauT and thetaT for Channel %d round %d/%d', iC, irounds, maxRounds));
%default starting values %default starting values
MO.par.tauT = 7; MO.par.tauT = 7;
...@@ -206,25 +225,26 @@ end ...@@ -206,25 +225,26 @@ end
if (dataC{1}(1,4) == 0) if (dataC{1}(1,4) == 0)
Chfits(2) = Chfits(1); Chfits(2) = Chfits(1);
end end
if (dataC{1}(1,6) == 0)
Chfits(3) = Chfits(1);
end
save(fullfile(outdir, [prefix '.mat']), 'Chfits'); save(fullfile(outdir, [prefix '.mat']), 'Chfits');
end end
% avoid cases were boundary of kappaVal is touched and use a value of kappa % avoid cases where boundary of kappaVal is touched and use a value of kappa
% = 5.8 instead this is the average values obtained for many experiments % = 5.8 instead this is the average values obtained for many experiments
posB = find(kappaVal==5.8); posB = find(kappaVal==5.8);
if all([Chfits.posBest] == [1 1]) || all([Chfits.posBest] == [lk lk]) if all([Chfits.posBest] == ones(1, length(Chfits))) || all([Chfits.posBest] == lk*ones(1, length(Chfits)))
posB = find(kappaVal==5.8); posB = find(kappaVal==5.8);
elseif (Chfits(1).posBest == 1) || (Chfits(1).posBest == lk) elseif (Chfits(1).posBest == 1) || (Chfits(1).posBest == lk)
posB = Chfits(2).posBest; posB = Chfits(2).posBest;
elseif (Chfits(2).posBest == 1) || (Chfits(2).posBest == lk) elseif (Chfits(2).posBest == 1) || (Chfits(2).posBest == lk)
posB = Chfits(1).posBest; posB = Chfits(1).posBest;
else else
posB = round(mean([Chfits(1).posBest,Chfits(2).posBest])); posB = round(mean([Chfits(1).posBest,Chfits(2).posBest]));
end end
focVolA = [Chfits(1).mfocRad(posB) Chfits(1).mfocVol(posB) Chfits(1).kappaVal(posB);... focVolA = [Chfits(1).mfocRad(posB) Chfits(1).mfocVol(posB) Chfits(1).kappaVal(posB);...
Chfits(2).mfocRad(posB) Chfits(2).mfocVol(posB) Chfits(2).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)]; sqrt(Chfits(1).mfocRad(posB)*Chfits(2).mfocRad(posB)) sqrt(Chfits(1).mfocVol(posB)*Chfits(2).mfocVol(posB)) Chfits(2).kappaVal(posB)];
focVolfile = fullfile(indir, 'focalVolume.txt'); focVolfile = fullfile(indir, 'focalVolume.txt');
fid = fopen(focVolfile, 'w'); fid = fopen(focVolfile, 'w');
fprintf(fid,'w0_um\tVol_fl\tkappa\r\n'); fprintf(fid,'w0_um\tVol_fl\tkappa\r\n');
...@@ -235,115 +255,115 @@ fclose(fid); ...@@ -235,115 +255,115 @@ fclose(fid);
%% save all rounds of fit %% save all rounds of fit
if process || force if process || force || replot
for iC = 1:2 for iC = 1:3
if (dataC{1}(1,iC*2) == 0) % do not process if channel has not been acquired if (dataC{1}(1,iC*2) == 0) % do not process if channel has not been acquired
continue continue
end
fid = fopen(fullfile(outdir, [prefix '_Ch' num2str(iC) '.txt']), 'w');
struct('N', 10, 'thetaT', 0.2, 'tauT', 7, 'f1', 1, 'tauD1', 20, 'alpha1', 1, 'tauD2', 5000, 'alpha2', 1, 'kappa', 5.5);
fprintf(fid,'N\tthetaT\ttauT\tf1\ttauD1\talpha1\ttauD2\talpha2\tkappa\tnorm\tR2\tNpt\tid_exp\r\n');
outPar = Chfits(iC).outPar;
for i = 1:size(outPar,1)
for j = 1:size(outPar,2)-1
fprintf(fid,'%.4e\t', outPar(i,j));
end end
fprintf(fid,'%.4e\r\n', outPar(i,end)); fid = fopen(fullfile(outdir, [prefix '_Ch' num2str(iC) '.txt']), 'w');
end struct('N', 10, 'thetaT', 0.2, 'tauT', 7, 'f1', 1, 'tauD1', 20, 'alpha1', 1, 'tauD2', 5000, 'alpha2', 1, 'kappa', 5.5);
fclose(fid); fprintf(fid,'N\tthetaT\ttauT\tf1\ttauD1\talpha1\ttauD2\talpha2\tkappa\tnorm\tR2\tNpt\tid_exp\r\n');
end outPar = Chfits(iC).outPar;
for i = 1:size(outPar,1)
%% update xml files and write a res file for j = 1:size(outPar,2)-1
if ~ZEN fprintf(fid,'%.4e\t', outPar(i,j));
for iD = 1:length(dataC)
node = xmlFA.readnode(xmlFAfiles{iD});
[path, fout, ext] = fileparts(xmlFAfiles{iD});
fname_out = fullfile(path, [fout '_opt' ext]);
idxs = MO.getDataIndexes(dataC{iD}, MO.tauBoundary);
Nfits = length(idxs(1):length(dataC{iD})-idxs(2));
fit = zeros(Nfits, 3);
par = repmat(struct2array(MO.par),3,1);
nrPar = length(par);
for iC = 1:2
if (dataC{iD}(1,iC*2) == 0) % do not process if channel has not been acquired
continue
end end
posExp = find(round(10*Chfits(iC).outPar(:,9))==round(10*Chfits(iC).kappaVal(posB))); fprintf(fid,'%.4e\r\n', outPar(i,end));
try
fit(:,iC) = MO.Gcomp(MO.parvecToparstruct(Chfits(iC).outPar(posExp(iD), 1:nrPar)), dataC{iD}(idxs(1):end-idxs(2),1));
catch
display('error')
end
par(iC, :) = Chfits(iC).outPar(posExp(iD), 1:nrPar);
end end
xmlFA.writepar(node, par, struct2array(MO.lb), struct2array(MO.hb)); fclose(fid);
xmlFA.writefit(node, fit, idxs, MO.model);
xmlFA.writenode(fname_out, node);
end
end
%% 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 (fl) 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 end
hf = setupFigure(iC, [100 100 400 800], 'off' ); %% update xml files and write a res file
clf if ~ZEN
% normalized Chi2 for iD = 1:length(dataC)
subplot(4,1,1) node = xmlFA.readnode(xmlFAfiles{iD});
hold('on'); [path, fout, ext] = fileparts(xmlFAfiles{iD});
posB = Chfits(iC).posBest; fname_out = fullfile(path, [fout '_opt' ext]);
errorbar(kappaVal, mean(Chfits(iC).normN,2), std(Chfits(iC).normN, 1,2), '-') idxs = MO.getDataIndexes(dataC{iD}, MO.tauBoundary);
plot(Chfits(iC).kappa, Chfits(iC).mNorm(posB), 'ro', 'MarkerFaceColor', 'r' ) Nfits = length(idxs(1):length(dataC{iD})-idxs(2));
ylim([min(Chfits(iC).normN(:)) max(Chfits(iC).normN(:))]) fit = zeros(Nfits, 3);
ylabel('normalized Chi2'); par = repmat(struct2array(MO.par),3,1);
xlabel('kappa') nrPar = length(par);
title(sprintf('Channel %d kappa %.1f; Vol %.3f fl; w0 %.3f um', iC, Chfits(iC).kappa, Chfits(iC).mfocVol(posB), Chfits(iC).mfocRad(posB))) for iC = 1:3
if (dataC{iD}(1,iC*2) == 0) % do not process if channel has not been acquired
continue
end
posExp = find(round(10*Chfits(iC).outPar(:,9))==round(10*Chfits(iC).kappaVal(posB)));
try
fit(:,iC) = MO.Gcomp(MO.parvecToparstruct(Chfits(iC).outPar(posExp(iD), 1:nrPar)), dataC{iD}(idxs(1):end-idxs(2),1));
catch
display('error')
end
par(iC, :) = Chfits(iC).outPar(posExp(iD), 1:nrPar);
end
xmlFA.writepar(node, par, struct2array(MO.lb), struct2array(MO.hb));
xmlFA.writefit(node, fit, idxs, MO.model);
xmlFA.writenode(fname_out, node);
end
end
% focal volume
subplot(4,1,2)
hold('on');
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('on');
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 %% plot results
subplot(4,1,4) for iC = 1:3
hold('on'); if (dataC{1}(1,iC*2) == 0) % do not process if channel has not been acquired
errorbar(kappaVal,mean(Chfits(iC).conc,2), std(Chfits(iC).conc,1,2), '-') continue
plot(Chfits(iC).kappa, Chfits(iC).mConc(posB), 'ro', 'MarkerFaceColor', 'r' ) end
ylabel('Concentration (nM)') % means are computed from mean tauD1 and mean N. Otherwise error
xlabel('kappa') % increases
saveas(hf,fullfile(outdir,[prefix '_Ch' num2str(iC) '.png']), 'png') for i = 1:lk
end % find indexes for given kappaValue
posExp = find(round(10*Chfits(iC).outPar(:,9))==round(10*kappaVal(i)));
% compute mean focal volume (fl) 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], 'off' );
clf
% normalized Chi2
subplot(4,1,1)
hold('on');
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 fl; w0 %.3f um', iC, Chfits(iC).kappa, Chfits(iC).mfocVol(posB), Chfits(iC).mfocRad(posB)))
% focal volume
subplot(4,1,2)
hold('on');
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('on');
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('on');
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
end end
end 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