Commit 78abb88d authored by Antonio Politi's avatar Antonio Politi

example for single pore data analysis

parent 5a2e2377
......@@ -43,3 +43,4 @@ matlabcode/FCS/example_data/Alexa/
expdata
.Rproj.user
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
Copyright (c) 2017, YoonOh Tak
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
function oimg = loadtiff(path)
% Copyright (c) 2012, YoonOh Tak
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% * Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
% * Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in
% the documentation and/or other materials provided with the distribution
% * Neither the name of the Gwangju Institute of Science and Technology (GIST), Republic of Korea nor the names
% of its contributors may be used to endorse or promote products derived
% from this software without specific prior written permission.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
tStart = tic;
warn_old = warning('off', 'all'); % To ignore unknown TIFF tag.
%% Check directory and file existence
path_parent = pwd;
[pathstr, ~, ~] = fileparts(path);
if ~isempty(pathstr) && ~exist(pathstr, 'dir')
error 'Directory is not exist.';
end
if ~exist(path, 'file')
error 'File is not exist.';
end
%% Open file
file_opening_error_count = 0;
while ~exist('tiff', 'var')
try
tiff = Tiff(path, 'r');
catch
file_opening_error_count = file_opening_error_count + 1;
pause(0.1);
if file_opening_error_count > 5 % automatically retry to open for 5 times.
reply = input('Failed to open the file. Do you wish to retry? Y/n: ', 's');
if isempty(reply) || any(upper(reply) == 'Y')
file_opening_error_count = 0;
else
error(['Failed to open the file ''' path '''.']);
end
end
end
end
%% Load image information
tfl = 0; % Total frame length
tcl = 1; % Total cell length
while true
tfl = tfl + 1; % Increase frame count
iinfo(tfl).w = tiff.getTag('ImageWidth');
iinfo(tfl).h = tiff.getTag('ImageLength');
iinfo(tfl).spp = tiff.getTag('SamplesPerPixel');
iinfo(tfl).color = iinfo(tfl).spp > 2; % Grayscale: 1(real number) or 2(complex number), Color: 3(rgb), 4(rgba), 6(rgb, complex number), or 8(rgba, complex number)
iinfo(tfl).complex = any(iinfo(tfl).spp == [2 6 8]);
if tfl > 1
% If tag information is changed, make a new cell
if iinfo(tfl-1).w ~= iinfo(tfl).w || ...
iinfo(tfl-1).h ~= iinfo(tfl).h || ...
iinfo(tfl-1).spp ~= iinfo(tfl).spp || ...
iinfo(tfl-1).color ~= iinfo(tfl).color || ...
iinfo(tfl-1).complex ~= iinfo(tfl).complex
tcl = tcl + 1; % Increase cell count
iinfo(tfl).fid = 1; % First frame of this cell
else
iinfo(tfl).fid = iinfo(tfl-1).fid + 1;
end
else
iinfo(tfl).fid = 1; % Very first frame of this file
end
iinfo(tfl).cid = tcl; % Cell number of this frame
if tiff.lastDirectory(), break; end;
tiff.nextDirectory();
end
%% Load image data
if tcl == 1 % simple image (no cell)
for tfl = 1:tfl
tiff.setDirectory(tfl);
temp = tiff.read();
if iinfo(tfl).complex
temp = temp(:,:,1:2:end-1,:) + temp(:,:,2:2:end,:)*1i;
end
if ~iinfo(tfl).color
oimg(:,:,iinfo(tfl).fid) = temp; % Grayscale image
else
oimg(:,:,:,iinfo(tfl).fid) = temp; % Color image
end
end
else % multiple image (multiple cell)
oimg = cell(tcl, 1);
for tfl = 1:tfl
tiff.setDirectory(tfl);
temp = tiff.read();
if iinfo(tfl).complex
temp = temp(:,:,1:2:end-1,:) + temp(:,:,2:2:end,:)*1i;
end
if ~iinfo(tfl).color
oimg{iinfo(tfl).cid}(:,:,iinfo(tfl).fid) = temp; % Grayscale image
else
oimg{iinfo(tfl).cid}(:,:,:,iinfo(tfl).fid) = temp; % Color image
end
end
end
%% Close file
tiff.close();
cd(path_parent);
warning(warn_old);
display(sprintf('The file was loaded successfully. Elapsed time : %.3f s.', toc(tStart)));
end
This diff is collapsed.
%% This file shows how to use 'saveastiff' function.
clearvars;
[~,~,Z] = peaks(100);
Z = single(Z);
Z_index = uint8((Z - min(Z(:))) * (255 / (max(Z(:)) - min(Z(:)))));
Z_color = uint8(ind2rgb(Z_index, hsv(256)*256));
Z_alpha = Z_color; Z_alpha(:,:,4) = repmat(0:1/100:1-1/100, 100, 1)*255;
Z_color_multiframe = reshape([Z_color(:)*0.2 Z_color(:)*0.6 Z_color(:)], 100, 100, 3, 3);
disp('8-bit, grayscale image');
clear options;
saveastiff(uint8(Z_index), 'Z_uint8.tif');
disp('Lossless LZW compression');
clear options;
options.comp = 'adobe';
saveastiff(uint8(Z_index), 'Z_uint8_compress.tif', options);
disp('Overwrite to an existing file');
clear options;
options.overwrite = true;
options.compress = 'adobe';
saveastiff(uint8(Z_index), 'Z_uint8_compress.tif', options);
disp('Disable message printing.');
clear options;
options.message = false;
saveastiff(uint8(Z_index), 'Z_uint8_Message.tif', options);
disp('16 bit, grayscale image');
clear options;
saveastiff(uint16(single(Z_index)/255*(2^16-1)), 'Z_uint16.tif');
disp('32 bit, grayscale image');
clear options;
saveastiff(uint32(single(Z_index)/255*(2^32-1)), 'Z_uint32.tif');
disp('32 bit single, grayscale image');
clear options;
saveastiff(Z, 'Z_single.tif');
disp('Color image');
clear options;
options.color = true;
saveastiff(Z_color, 'Z_rgb.tif', options);
disp('Color image with alpha channel');
clear options;
options.color = true;
saveastiff(Z_alpha, 'Z_alpha.tif', options);
disp('Save each R, G and B chanels of the color image, separately.');
clear options;
saveastiff(Z_color, 'Z_rgb_channel.tif');
disp('Save the multi-frame RGB color image');
clear options;
options.color = true;
saveastiff(Z_color_multiframe, 'Z_rgb_multiframe.tif', options);
disp('32 bit single, 50x50x50 volume data');
clear options;
saveastiff(single(rand(50, 50, 50)), 'volume_50x50x50.tif');
disp('Append option is ignored if path dose not exist.');
clear options;
options.append = true;
saveastiff(Z_index, 'Z_uint8_append.tif', options);
disp('You can append any type of image to an existing tiff file.');
clear options;
options.append = true;
saveastiff(single(rand(10, 10, 3)), 'Z_uint8_append.tif', options);
options.color = true;
saveastiff(Z_color_multiframe, 'Z_uint8_append.tif', options);
disp('Save image to a sub directory');
clear options;
saveastiff(uint8(Z_index), 'sub1/sub2/Z_uint8.tif');
options.append = true;
saveastiff(uint8(Z_index), 'sub1/sub2/Z_uint8.tif', options);
disp('Save complex number images');
saveastiff(fft2(Z_index), 'Z_fft2_gray.tif');
saveastiff(fft2(Z_color), 'Z_fft2_color.tif', options);
clear options;
options.color = true;
options.append = true;
saveastiff(fft2(Z_index), 'Z_fft2_append.tif');
saveastiff(fft2(Z_color), 'Z_fft2_append.tif', options);
saveastiff(fft2(Z_color), 'Z_fft2_append.tif', options);
saveastiff(fft2(Z_alpha), 'Z_fft2_append.tif', options);
disp('Load multiframe tiff');
multiframe = loadtiff('volume_50x50x50.tif');
complexframe = loadtiff('Z_uint8_append.tif');
subdir = loadtiff('sub1/sub2/Z_uint8.tif');
fftgray = loadtiff('Z_fft2_gray.tif');
fftcolor = loadtiff('Z_fft2_color.tif');
fftappend = loadtiff('Z_fft2_append.tif');
% Big Tiff File (64 bit)
[~,maxArraySize]=computer;
is64bitMatlab = maxArraySize> 2^31;
if is64bitMatlab
% Use 'big' option to save larget than 4GB files .
disp('24000 by 24000, 32 bit single, 2GB image');
clear options;
options.big = true; % Use BigTIFF format
saveastiff(single(zeros(24000, 24000)), 'BigTiff(2GB+2GB).btf', options);
disp('Append 2GB image to the existing file.');
options.append = true;
saveastiff(single(zeros(24000, 24000)), 'BigTiff(2GB+2GB).btf', options); % 4GB Big TIFF file
bigtiff = loadtiff('BigTiff(2GB+2GB).btf');
end
......@@ -20,7 +20,7 @@
% [prctile(ratio, 50) prctile(ratio, 25) prctile(ratio, 75)]
%
% % ratio 0.385 - 0.617
% % in elife use 0.5
% % in elife used 0.5
% %% ratio mature to total in outer core 160229-EMdata-density-precursor.xlsx at 19.2 min
% m = 11.3;
% v = 2.3^2;
......
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MCMC fitting/simulated annealing fitting %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function MCMC_fun_classes(MO, par_idx, time, time_data, data, filename, Tmax, maxstep)
MCMCSteps = 50000;
T_lambda = 0.995;
t = clock;
J = mod(round(cumprod([10000,60,60,24,31,12])*t(6:-1:1)'),2^32);
rng(J);
%maximal size of changes in one parameter
minstep = 1e-5;
DISPLAY = 1;
Npar = length(par_idx);
%initiate MCMC run
partofit_values = MO.get_parvec;
parbest = partofit_values(par_idx);
partofit_values = partofit_values(par_idx);
lbound = MO.get_parvec_lb;
hbound = MO.get_parvec_hb;
lbound = lbound(par_idx);
hbound = hbound(par_idx);
[res, mlogL] = MO.dist_par(partofit_values, par_idx, time, time_data, data);
%-log(Likelihood)
mlogL = 1/2*mlogL;
bestmL = mlogL;
poww = 1;
%compute width of steps
%width = width_comp(dist, partofit, partofit_names, par, data, bg, dataname,vara modelfun, mlogL, varargin{:});
widthp = parbest;%.*width;
widthl = 0.9^poww*widthp;
widthl(widthl>maxstep) = maxstep;
T = Tmax;
if T_lambda > 1
error(' lamba, rate for deacreae in temperature < 1')
end
acceptance = 0;
acceptance50 = 0;
fittedpar = [];
fid = fopen(filename, 'w');
fprintf(fid, '%% MCMC run for model %s\n', MO.name);
fprintf(fid, '%% idx fitted %d', par_idx(1));
for i = 2:length(par_idx)
fprintf(fid, ' %d', par_idx(i));
end
fprintf(fid, '\n');
fclose(fid);
for i = 1:MCMCSteps
T = Tmax*T_lambda^i;
if(floor(i/100)-i/100+1==1) %compute step size every 1000 points
%width = width_comp(dist, partofit, partofit_names, par, data, bg, dataname, modelfun, mlogL, varargin);
end
%widthp = partofit_values;%.*width;
%acceptance should remain between 10 and 30%
if(floor(i/50)-i/50+1==1)
acceptance50=acceptance/50;
llim_accept=0.2;%+0.6/(par.fit.Tmax + 1)*par.fit.T;
hlim_accept=0.3;%+0.5/(par.fit.Tmax + 1)*par.fit.T;
if DISPLAY
sprintf('step %i Chi2 %.2f acceptance %.2f width_m %.2f poww %i llim_accept %.2f hlim_accept %.2f Temperature %f\n',...
i, 2*mlogL, acceptance/50, widthl(1),poww,...
llim_accept,hlim_accept, T)
end
if(acceptance/50 < llim_accept) & T < 0.1
poww = poww+1;
if poww > 5
poww = 5;
end
end
if(acceptance/50 > hlim_accept) & T < 0.1
poww = poww - 1;
if poww < -5
poww = -5;
end
end
acceptance=0;
end
widthl = widthp*0.9^poww;
widthl(widthl>maxstep) = maxstep;
widthl(widthl<minstep) = minstep;
%check that boundaries are respected
success = 0;
while ( ~success )
partofit_values_new = partofit_values + widthl.*randn(size(widthl,1),size(widthl,2));
if ~isempty(lbound) && ~isempty(hbound)
partofit_values_new(partofit_values_new < lbound) = 1*lbound(partofit_values_new < lbound) + 2*rand*lbound(partofit_values_new < lbound);
partofit_values_new(partofit_values_new > hbound) = 1*hbound(partofit_values_new > hbound) - 0.2*rand*hbound(partofit_values_new > hbound);
else
error('MCMC_fun provide boundaries for parameters');
end
try
[res, mlogL_new] = MO.dist_par(partofit_values_new, par_idx, time, time_data, data);
success = 1;
catch
display(['MCMC_fun.m: integration of ' MO.type ' failed at MCMC_step ' num2str(i) ' with present set of parameters try a different one']);
success = 0;
partofit_values(partofit_values < 100*lbound) = 100*lbound(partofit_values < 100*lbound);
end
end
mlogL_new = 1/2*mlogL_new;
if DISPLAY
%display([i; partofit_new; mlogL_new]);
end
%Update best set
if bestmL > mlogL_new
bestmL = mlogL_new;
parbest = partofit_values_new;
end
if(mlogL >= mlogL_new||(log(rand)<1/( 1 + T )*(mlogL - mlogL_new))) %accept if parameters are better or with certain probability rand
partofit_values = partofit_values_new;
mlogL = mlogL_new;
acceptance = acceptance+1;
end
if((floor(i/10)-i/10+1==1)&& T < 1e-3 )%save every 10 points
%sprintf('step i %d mlogL %f',i,mlogL)
fid = fopen(filename, 'A');
par = MO.get_parvec();
par(par_idx) = partofit_values;
for ipar = 1:length(par)
fprintf(fid, '%.2f\t', par(ipar));
end
fprintf(fid, '%.2f\t%.2f\t%.2f\t%.2f\n', 2*mlogL, acceptance50, T, poww);
fclose(fid);
end
end
% print best set at the end in case it has been lost between updates
fid = fopen(filename, 'A');
par = MO.get_parvec();
par(par_idx) = parbest;
for ipar = 1:length(par)
fprintf(fid, '%.2f\t', par(ipar));
end
fprintf(fid, '%.2f\t%.2f\t%.2f\t%.2f\n', 2*bestmL, acceptance50, T, poww);
fclose(fid);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute relative width of steps in the parameterspace %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function width = width_comp(dist, partofit, partofit_names, par, data, bg, dataname, modelfun, mlogL, varargin)
Npar = length(partofit);
jac = zeros(1,Npar);
for j = 1:Npar
partofit_width = partofit;
partofit_width(j) = partofit_width(j) + 0.05*partofit_width(j);
[res, mlogL_width] = dist(partofit_width, partofit_names, par, data, bg, dataname, modelfun, varargin{:});
mlogL_width = 1/2*mlogL_width;
jac(j) = abs((mlogL-mlogL_width)/0.05);
end
width=1-jac/(sum(jac));
%Normalize vector
width=width'/norm(width);
end
classdef abs_single_npc < handle
properties (Access = protected)
par_ini;
par_kin;
end
methods
function MO = abs_single_npc(par_ini, par_kin)
MO.par_ini = par_ini;
MO.par_kin = par_kin;
end
function setpar_ini(MO, par)
assert(MO.assert_par_ini(par), 'Number of parameters or range for par_ini is not correct');
MO.par_ini = par;
end
function setpar_kin(MO, par)
assert(MO.assert_par_kin(par), 'Number of parameters or range for par_kin is not correct');
MO.par_kin = par;
end
function X = kin_ini_conv_mom(MO, time, imom)
% This is X(t) = int_0^\infty g(t') \int_0^\infty f(x, t - t') x^imom d x dt'
X = conv(MO.ini_pdf(time), MO.kin_mom_pdf(time, imom));
X = X(1:length(time));
end
function plot_ini_cdf(MO, time, ifig)
figure(ifig);
clf;
plot(time, MO.ini_cdf(time));
xlabel('Time (s)');
ylabel('Initiation cdf');
end
function plot_ini_pdf(MO, time, ifig)
figure(ifig)
clf;
plot(time, MO.ini_pdf(time));
xlabel('Time (s)');
ylabel('Initiation pdf');
end
function plot_kin_mom_pdf(MO, time, imom, ifig)
figure(ifig);
clf;
plot(time, MO.kin_mom_pdf(time, imom));
xlabel('Time (s)');
ylabel(sprintf('NUP kinietics moment %d', imom));
end
function plot_combined_dist(MO, time, ifig)
dt = time(2)-time(1);
figure(ifig);
clf;
subplot(3,1,1)
plot(time, MO.kin_mom_pdf(time, 1));
xlabel('Time (min)');
ylabel('NUP kinetics 1st moment');
subplot(3,1,2)
plot(time, MO.ini_pdf(time));
xlabel('Time (min)');
ylabel('Initiation pdf');
subplot(3,1,3)
kin_ini_conv = MO.kin_ini_conv_mom(time, 1);
plot(time, kin_ini_conv(1:length(time))*dt);
xlabel('Time (min)');
ylabel('Average NPCs occupacy');
end
function par = get_par_ini(MO)
par = MO.par_ini;
end
function par = get_par_kin(MO)
par = MO.par_kin;
end
end
methods(Abstract)
bool = assert_par_kin(MO, par);
bool = assert_par_ini(MO, par);
vec = ini_cdf(MO, time) % int_0^t g(t) dt cumulative distribution function for the initiation. Returns a vector
vec = ini_pdf(MO, time) % g(t) probability density function for the initiation. Returns a vector
mat = kin_pdf(MO, time) % f(x,t) probability density function for the kinetics where x is occupancy. Should return a matrix
vec = kin_mom_pdf(MO, time, imom) % Moments of the distribution int_0^\infty f(x,t) x^imom. F(t) = kin_mom_pdf(time, 1)
end
end
\ No newline at end of file
This diff is collapsed.
classdef single_npc_logn_hill < abs_single_npc
% par_ini(1) : Mean of distribution
% par_ini(2) : std of distribution
properties
% Set of boundaries for the parameters
par_kin_lb = [0.5; 1];
par_kin_hb = [10; 200];
par_ini_lb = [1; 2];
par_ini_hb = [200; 200];
name = 'single_npc_logn_hill';
end
methods
function MO = single_npc_logn_hill(par_ini, par_kin)
display('single_npc_logn_hill: Log normal initiation, hill sigmoidal kinetics')
MO = MO@abs_single_npc(par_ini, par_kin);
MO.assert_par_kin(par_kin);
MO.assert_par_ini(par_ini);
end
function bool = assert_par_ini(MO, par)
bool = (length(par) == 2);
if bool == 0
display('par_ini length == 2, mean and std of initiation lognormal distribution.');
end
end
function bool = assert_par_kin(MO, par)
bool = (length(par) == 2);
if bool == 0
display( 'par_kin length == 2, slope and offset of lognormal distribution');
end
end
function vec = ini_pdf(MO, time)
%% g(t) probability density function for the initiation. Returns a vector
% convert parmaters to entries for log-normal distribution
mu = log((MO.par_ini(1)^2)/sqrt(MO.par_ini(2)+MO.par_ini(1)^2));
sigma = sqrt(log(MO.par_ini(2)/(MO.par_ini(1)^2)+1));
% [M, V]= lognstat(mu, sigma); M == par_ini(1), V == par_ini(2)
pd = makedist('Lognormal',mu, sigma); % Create object for distribution
vec = pdf(pd, time);
end
function vec = ini_cdf(MO, time)
%% int_0^t g(t) dt cumulative distribution function for the initiation. Returns a vector
% convert parmaters to entries for log-normal distribution
mu = log((MO.par_ini(1)^2)/sqrt(MO.par_ini(2)+MO.par_ini(1)^2));
sigma = sqrt(log(MO.par_ini(2)/(MO.par_ini(1)^2)+1));
% Create object for distribution
pd = makedist('Lognormal',mu, sigma);
vec = cdf(pd, time);
end
function mat = kin_pdf(MO, time)
%% NOT IMPLEMENTED f(x,t) probability density function for the kinetics where x is occupancy. Should return a matrix
display('Distribution of NUP kinetics has not been defined')
mat = [];