diff --git a/CellSortExportTifStack.m b/CellSortExportTifStack.m new file mode 100644 index 0000000..3e850c2 --- /dev/null +++ b/CellSortExportTifStack.m @@ -0,0 +1,5 @@ +function CellSortExportTifStack(data, fn) +for i=1:size(data,3) + imwrite(data(:,:,i),fn, 'WriteMode', 'append'); +end +end \ No newline at end of file diff --git a/CellSortLoadSingleFiles.m b/CellSortLoadSingleFiles.m new file mode 100644 index 0000000..8313231 --- /dev/null +++ b/CellSortLoadSingleFiles.m @@ -0,0 +1,24 @@ +function [data_g, data_r] = CellSortLoadSingleFiles(green, red) +fileListGreen = dir(green); +fileListRed = dir(red); + +numFilesGreen = size(fileListGreen,1); +numFilesRed = size(fileListRed,1); +if numFilesGreen ~= numFilesRed +disp('Green and Red Channel have different number of Images, take the smaller value..'); +end +numFiles = min([numFilesRed, numFilesGreen]); + +%read first frame to get resolution fily type etc.. +temp = readTiff(fileListGreen(1).name); +data_g = zeros(size(temp,1), size(temp,2), numFiles, class(temp)); +temp = readTiff(fileListRed(1).name); +data_r = zeros(size(temp,1), size(temp,2), numFiles, class(temp)); + +for i=1:numFiles + temp = readTiff(fileListGreen(i).name); + data_g(:,:,i) = temp; + temp = readTiff(fileListRed(i).name); + data_r(:,:,i) = temp; +end +end \ No newline at end of file diff --git a/CellSortShiftCorrection.m b/CellSortShiftCorrection.m new file mode 100644 index 0000000..c06ea8f --- /dev/null +++ b/CellSortShiftCorrection.m @@ -0,0 +1,13 @@ +function [data_g, data_r] = CellSortShiftCorrection(data_g, data_r) +%% shiftcorrection +%import green channel +disp('Performing Shift Correction..'); +%shift correction +firstRedFrame = data_r(:,:,1); +parfor r = 2:size(data_g, 3) + shiftVector = findshift(firstRedFrame, data_r(:,:,r), 'iter'); + data_r(:,:,r) = shift(data_r(:,:,r),shiftVector); + data_g(:,:,r) = shift(data_g(:,:,r),shiftVector); +end +disp('finished'); +end \ No newline at end of file diff --git a/CellSortTemporalBinning.m b/CellSortTemporalBinning.m new file mode 100644 index 0000000..0e6844c --- /dev/null +++ b/CellSortTemporalBinning.m @@ -0,0 +1,4 @@ +function [data_g, data_r] = CellSortTemporalBinning(data_g, data_r, binningFactor) +data_g = binning(data_g, bininngFactor, 'mean'); +data_r = binning(data_r, bininngFactor, 'mean'); +end diff --git a/binning.m b/binning.m new file mode 100644 index 0000000..adcb241 --- /dev/null +++ b/binning.m @@ -0,0 +1,149 @@ +function A = binning(A,bin,fun_str) +% +% A=binning(A,bin,fun) +% +% This function rearranges the array A according to the specifications +% of the vector bin. In the first dimension of the generated temporary +% array the bin elements from all dimension are arranged. By applying +% functions like 'sum', 'mean' or 'std' to the first dimension, a real +% binning ('sum'), a averaging ('mean') and error calculations can be +% done. +% +% A=binning(A,bin) The input array 'A' is binned by the vector 'bin'. +% If 'bin' is a scalar, 'bin' elements in each dimenion of 'A' are +% summarized. +% +% Note: The input array 'A' is cut to the maxial integer factor of +% 'bin'. +% +% A=binning(A,bin,fun) The function 'fun' is applied to the first +% dimension of the rearranged array. Beside the real binning with the +% function 'sum' a number of additional functions can be applied. Up +% to now the following functions are implemented: +% 'sum', 'mean', 'prod' +% +% Special comment: +% The functionality can be easily extended by implementing additional +% functions, or due to the time consuming permutations by the +% combination of appropriate functions like 'mean' and 'std'. +% +% Author: A. Zeug (last modified 19.07.2007) +% Version: 2.0.0 + +% comments to: azeug@gwdg.de + +tic; t1=toc; +%% Input validation +if nargin < 3 + fun_str='sum'; + fun_type = 1; +elseif sum(strcmp(fun_str,{'sum';'mean';'prod'}))>0 + fun_type = 1; +elseif sum(strcmp(fun_str,{'std';'var'}))>0 + fun_type = 2; +else + error('Function not implemented yet') +end +fun=str2func(fun_str); % defines function type + +if nargin < 1 + help binning + return +elseif nargin < 2 + error('Please specify binning vector.') +end + +%% bin vector dimension +trans = false; +dims = ndims(A); +if isvector(A) + dims = 2; + if size(A,2)==1; trans = true; end +end +if length(bin)==1; + bin = bin * ones(1, dims); +elseif length(bin)~=dims; + error('Wrong length of binning vector.') +end + +%% resizing of A to fit bin +sz=size(A); +if dims==1; sz=length(A); end +bin=min([sz;bin]); +for i=1:dims + index{i}=1:(bin(i)*floor(sz(i)/bin(i))); +end +A = A(index{:}); + +%% reshaping A +sz = size(A); if dims==1; sz = length(A); end +v = []; sz_new = sz./bin; +for i=1:dims + v=[v, bin(i), sz_new(i)]; +end +A=reshape(A,v); + +%% Calculate binning +if fun_type==1 + % standard type 'sum', 'mean','prod' + for i=1:dims + if isequal(fun, @sum) + A=fun(A,(2*i-1),'native'); + else + A=fun(A,(2*i-1)); + end + end +else + % type 'std', 'var' + % need to reshape again to calc function at ones + v = []; + for i=1:dims + v=[v, i, dims+i]; + end + A = permute(A, v); + A = reshape(A, [prod(bin) sz_new]); + try + A = fun( double(A) ); + catch + t=toc; fprintf('(%07.3fs) binning: direct calculation failed, start subsequential method (%0.5fs)\n',t,t-t1); t1=t; + try + A1 = zeros(sz_new,'double'); + fun1 = @double; + catch + t=toc; fprintf('(%07.3fs) binning: calculating %s with double precession failed, try single (%0.5fs)\n',t,func2str(fun),t-t1); t1=t; + A1 = zeros(sz_new, 'single'); + fun1 = @single; + end + h = waitbar(0,'Please wait...','Name','binning: Plese Wait!'); + % subsequential calculation + sub_dim = find(sz_new == max(sz_new)); + for i = 1:dims; cc{i}=':'; end + step = floor(1E6 / (prod(sz_new(sz_new < max(sz_new)))*prod(bin))); + for ii = 0:step:sz_new(sub_dim) % subsequentially... + waitbar(ii/sz_new(sub_dim),h,['Calculated binning: ' num2str(100*ii/sz_new(sub_dim),'%0.2f') '%']); + ii_end = min(sz_new(sub_dim), ii+step) - ii; + cc{sub_dim}=(1:ii_end)+ii; + A1(cc{:}) = fun(fun1(A(:,cc{:}))); + end + A = A1; + delete(h); drawnow; + t=toc; fprintf('(%07.3fs) binning: calculation with %s precession susseed (%0.5fs)\n',t,func2str(fun1),t-t1); t1=t; + end +end + +% +%A = reshape(A, sz_new); + +if dims==1 + A=squeeze(A); + if trans==true + A=A'; + end +else + A = reshape(A, sz_new); + +end + + + + \ No newline at end of file diff --git a/main.m b/main.m index f785465..1d08362 100644 --- a/main.m +++ b/main.m @@ -1,6 +1,18 @@ -%% filename -fn_green = 'green_13.tiff'; -fn_red = 'red_13.tiff'; +%% load and prepare all single tiffs +%importing +[data_g,data_r] = CellSortLoadSingleFiles('ChanA*.tif', 'ChanB*.tif'); +%% shiftCorrection +%generates shift corrected files for the green and red channels, needs +%DIPlib (available for free) +[data_g, data_r] = CellSortShiftCorrection(data_g,data_r); +%temporal binning +binningFactor = [1 1 2]; %x, y, t +[data_g, data_r] = CellSortTemporalBinning(data_g, data_r, binningFactor); +%% export +fileName_g = 'green.tif'; +fileName_r = 'red.tif'; +CellSortExportTifStack(data_g, fileName_g); +CellSortExportTifStack(data_r, fileName_r); %% merge files fn_green = 'merged_green.tiff'; fn_red = 'merged_red.tiff'; @@ -8,11 +20,6 @@ CellSortMergeTiff(fileList, fn_green); fileList = {'red_13.tiff', 'red_14.tiff','red_15.tiff','red_1.tiff','red_2.tiff','red_3.tiff'}; CellSortMergeTiff(fileList, fn_red); - -%% shiftCorrection -%generates shift corrected files for the green and red channels, needs -%DIPlib (available for free) -[fn_green, fn_red, data_g, data_r] = shiftCorrection(fn_green,fn_red, true); %% load files %if you haven't used the shiftCorrection you have to load the files extra fn_green = 'green_13_corr.tif'; diff --git a/main_old.m b/main_old.m new file mode 100644 index 0000000..f785465 --- /dev/null +++ b/main_old.m @@ -0,0 +1,134 @@ +%% filename +fn_green = 'green_13.tiff'; +fn_red = 'red_13.tiff'; +%% merge files +fn_green = 'merged_green.tiff'; +fn_red = 'merged_red.tiff'; +fileList = {'green_13.tiff', 'green_14.tiff','green_15.tiff','green_1.tiff','green_2.tiff','green_3.tiff'}; +CellSortMergeTiff(fileList, fn_green); +fileList = {'red_13.tiff', 'red_14.tiff','red_15.tiff','red_1.tiff','red_2.tiff','red_3.tiff'}; +CellSortMergeTiff(fileList, fn_red); + +%% shiftCorrection +%generates shift corrected files for the green and red channels, needs +%DIPlib (available for free) +[fn_green, fn_red, data_g, data_r] = shiftCorrection(fn_green,fn_red, true); +%% load files +%if you haven't used the shiftCorrection you have to load the files extra +fn_green = 'green_13_corr.tif'; +fn_red = 'red_13_corr.tif'; +data_g = readTiff(fn_green); +data_r = readTiff(fn_red); +%% 1. PCA +%calculate PCAs +nPCs = 100; +flims = []; +dSamp = [1 1]; %temporal, spatial downsampling +[mixedsig, mixedfilters, CovEvals, covtrace, movm, movtm] = ... + CellsortPCA(fn_green, flims, nPCs, dSamp, [], [], data_g); + +%% 2a. Choose PCs +%choose PCAs, type f/b in command line to go back or forward. type two number +%to select the range "1 ENTER 25 ENTER" would select the PCAs from 1 to 25 +[PCuse] = CellsortChoosePCs(fn_green, mixedfilters); + +%% 2b. Plot PC spectrum +%optionally, you can display the PC spectrum.. +CellsortPlotPCspectrum(fn_green, CovEvals, PCuse) + +%% 3a. ICA +%calculate ICAs, set mu to a value between 0 and 1, where 1 means pure +%temporal identification and 0 pure spatial. A value of 0.1-0.2 seems to be +%best.. make sure that the algorithm converges (see command line). If not +%run it again +nIC = length(PCuse); +mu = 0.2; + +[ica_sig, ica_filters, ica_A, numiter] = CellsortICA(mixedsig, mixedfilters, CovEvals, PCuse, mu, nIC); + +%% 3b. Plot ICs +%optionally, plot ICs (you have to use at least 20.. or change the number +%on the buttom +tlims = []; +dt = 0.1; + +figure(2) +CellsortICAplot('series', ica_filters, ica_sig, movm, tlims, dt, [], [], [1:20]); + +%% 4a. Segment contiguous regions within ICs +%search for indepent segments, you have to adjust these values well: +%smwidth: width of a smoothing filter (1 should be okay) +%thresh: threshold for normalized standard deviation in ICs (verey +% sensitive, standard: 2) +%arealims: limit for the detected areas in pixeks: [min max] or [min] +%plotting: 1 = yes or 0 = no.. shall a plot be shown or not? + +smwidth = 1; +thresh = 2; +arealims = [250]; +plotting = 1; + +[ica_segments, segmentlabel, segcentroid] = CellsortSegmentation(ica_filters, smwidth, thresh, arealims, plotting); + +%% 4b. CellsortApplyFilter +%filter signal and extract regions signals +subtractmean = 1; + +cell_sig = CellsortApplyFilter(fn_green, ica_segments, flims, movm, subtractmean, data_g); + +%% 5. CellsortFindspikes +%optionally, find spikes in signal +deconvtau = 0; +spike_thresh = 2; +normalization = 1; + +[spmat, spt, spc] = CellsortFindspikes(ica_sig, spike_thresh, dt, deconvtau, normalization); + +%% Show results +%optionally, plot signals +figure(2) +CellsortICAplot('series', ica_filters, ica_sig, movm, tlims, dt, 1, 2, [1:20], spt, spc); +%% prepare georgis data for visualization +%optionally, if you want to use the alternative evaluation georgi was using, import the +%mask and run this section +fn='green_13.tiff'; +info = imfinfo(fn); +num_images = numel(info); +data_g = readTiff(fn); +%make mask logical +mask1 = logical(mask.*(mask > 0)); +ica_segments = permute(mask, [3 1 2]); +mask1 = repmat(mask, [1 1 1 size(data_g,3)]); +%generate cell_sig for every feature +cellsig = zeros(size(mask1,3), size(data_g,3)); +for i=1:size(mask1,3) + cell_sig(i,:) = squeeze(mean(mean(squeeze(mask1(:,:,i,:)).*data_g(:,:,:),1),2)); +end +clear mask +%% show signal of cells (cell_sig) +%caluclate correlation matrix for cell_sig +[ica_segments, cluster, subPlotCols, subPlotWidth, sortedIdx] = calcCorrEff(cell_sig, ica_segments); +alternativePlot(ica_segments, cluster, cell_sig, subPlotCols, sortedIdx, data_r); + +%% let user deselect and merge rois +%merge -> change cluster ->recereate plots +merge = {[1,2,3], [5,6]}; +delete = [8,4,11]; +figure; +cluster_new = cluster; +%delete +for i = 1:size(delete,2) + cluster_new{delete(i)} = []; +end +%merge +for i = 1:size(merge, 2) + for j = 2:size(merge{i},2) + cluster_new{merge{i}(1)} = [cluster_new{merge{i}(1)}, cluster_new{merge{i}(j)}]; + cluster_new{merge{i}(j)} = []; + end +end +cluster_new = cluster_new(~cellfun(@isempty, cluster_new)); +subPlotCols = size(cluster_new, 2); + +%redraw +alternativePlot(ica_segments, cluster_new, cell_sig, subPlotCols, sortedIdx, data_r); \ No newline at end of file diff --git a/readTiff.m b/readTiff.m index 33fcc4f..b7a8e05 100644 --- a/readTiff.m +++ b/readTiff.m @@ -1,7 +1,15 @@ function data = readTiff(filename) info = imfinfo(filename); + +if info.BitDepth == 16 || info.BitDepth == 12 + filetype = 'uint16'; +elseif info.BitDepth == 8 + filetype = 'uint8'; +else + filetype = 'uint16'; +end num_images = numel(info); -data = zeros(info(1).Width, info(1).Height, num_images); +data = zeros(info(1).Width, info(1).Height, num_images, filetype); t_g = Tiff(filename, 'r'); for k = 1:num_images t_g.setDirectory(k);