diff --git a/functions/CellsortApplyFilter.m b/functions/CellsortApplyFilter.m index 825db2c..4eed142 100644 --- a/functions/CellsortApplyFilter.m +++ b/functions/CellsortApplyFilter.m @@ -1,4 +1,4 @@ -function cell_sig = CellsortApplyFilter(fn, ica_segments, flims, movm, subtractmean, data_g) +function cell_sig = CellsortApplyFilter(ica_segments, flims, data_g) % cell_sig = CellsortApplyFilter(fn, ica_segments, flims, movm, subtractmean) % %CellsortApplyFilter @@ -21,60 +21,58 @@ % if (nargin<3)||isempty(flims) - nt = size(data_g, 3)%tiff_frames(fn); + nt = size(data_g, 3);%tiff_frames(fn); flims = [1,nt]; else nt = diff(flims)+1; end -if nargin<5 - subtractmean = 1; -end [pixw,pixh] = size(data_g(:,:,1));%size(imread(fn,1)); -if (nargin<4)||isempty(movm) - movm = ones(pixw,pixh); -else - movm = double(movm); -end -movm(movm==0) = 1; % Just in case there are black areas in the average image -k=0; +% k=0; cell_sig = zeros(size(ica_segments,1), nt); -ica_segments = reshape(ica_segments, [], pixw*pixh); - -fprintf('Loading %5g frames for %d ROIs.\n', nt, size(ica_segments,1)) -tic -while k0),2); + temp2=temp2/quotient; + cell_sig(seg,:) = temp2; end - function j = tiff_frames(fn) - % - % n = tiff_frames(filename) - % - % Returns the number of slices in a TIFF stack. - % - % Modified April 9, 2013 for compatibility with MATLAB 2012b +% for i = 1:size(data_g,2) +% for seg = 1:size(ica_segments,1) +% sig = data_g(:,i).*cast(squeeze(ica_segments(seg,:)),dataType)'; +% temp = mean(sig(sig~=0)); +% cell_sig(seg,i) = temp; +% end +% end - j = length(imfinfo(fn)); - end -end +end \ No newline at end of file diff --git a/functions/alternativePlot.m b/functions/alternativePlot.m index e38bcf8..55f03c4 100644 --- a/functions/alternativePlot.m +++ b/functions/alternativePlot.m @@ -26,7 +26,12 @@ 'value',0, 'min',0, 'max', ceil(size(cell_sig_merged{1},2)/frameSteps),... 'SliderStep', [1 1]/(ceil(size(cell_sig_merged{1},2)/frameSteps)),'Callback', @redrawPlot,... 'units','normalized','Position',[0.6 0.01 0.3 0.03]); - +imagetxt = uicontrol('units','normalized','Style','text',... + 'Position',[0.1 0.04 0.3 0.03],... + 'String','Showing All ROIs'); +plottxt = uicontrol('units','normalized','Style','text',... + 'Position',[0.6 0.04 0.3 0.03],... + 'String','Showing All Frames'); function redrawPlot(~,~) val1 = round(imageSlider.Value); @@ -56,7 +61,8 @@ function redrawPlot(~,~) kStart = 1; kEnd = size(cell_sig_merged{1},2); end - + imagetxt.String = ['Showing ROI ' num2str(iStart) ' to ' num2str(iEnd) ' of ' num2str(subPlotCols)]; + plottxt.String = ['Showing Frame ' num2str(kStart) ' to ' num2str(kEnd) ' of ' num2str(size(cell_sig_merged{1},2))]; drawEverything(iStart, iEnd, kStart, kEnd); end @@ -176,5 +182,7 @@ function drawEverything(imageStart, imageEnd, plotStart, plotEnd) vline(frames(i)+sweepDelay+s*sweepInterval+sweepDuration, 'r'); end end + xlabel('Frames'); + ylabel('Intensity + Arbitrary Offset'); end end \ No newline at end of file diff --git a/main.asv b/main.asv new file mode 100644 index 0000000..f308ac1 --- /dev/null +++ b/main.asv @@ -0,0 +1,144 @@ +%% load files and do cross correlation alignment +%Choose input folder, where all trial folders are and output folder where +%shift corrected tiff stacks for each trial shall be generated. You can +%also choose a binning factor [x, y, t]. nSubIm ist the number of +%iterations the shift correction will do. The more the better, but of +%course it is more time consuming. The standard value here is 10. +inFolder = 'F:\Arbeit\testData\PCAICA'; +outFolder = 'F:\Arbeit\testData\PCAICA\corr'; +binningFactor = [1, 1, 2]; +nSubIm = 10; +CellSortXCorr(inFolder, outFolder, nSubIm, binningFactor); +%% load files +%go to corr folder +outFolder = 'F:\Arbeit\testData\PCAICA'; +%shifft correction on all frames +%[data_g, data_r, frames, fn_green, fn_red, fileNamesGreen, fileNamesRed] = CellSortLoadTrials(outFolder); +%shift correction only on trials +[data_g, data_r, frames, fn_green, fn_red, fileNamesGreen, fileNamesRed] = CellSortLoadTrialsShiftCorrOnTrials(outFolder); +%If you want to load single files, use the functions below. +% fn_green = 'data_g.tif'; +% fn_red = 'data_r.tif'; +% data_g = readTiff(fn_green); +% data_r = readTiff(fn_red); +%% Crop Data +%If you want to crop the data to get rid of the shift correction artifacts, +%type in the intervals of pixels in x and y directions you want to keep and +%evaluate this section. Use matVis(data_r) to look at your data. +xInterval = 40:380; +yInterval = 40:380; + +data_g = data_g(xInterval,yInterval,:); +data_r = data_r(xInterval,yInterval,:); +%% 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, data_g); + +%OR specify the PCAs to use manually: +%PCuse = [7,8,18,20,21,25,27,30,32,35,36,46]; +%% 2b. Plot PC spectrum +%optionally, you can display the PC spectrum.. + +CellsortPlotPCspectrum(fn_green, CovEvals, PCuse, data_g) + +%% 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:size(ica_sig,1)]); + +%% 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 = 1.3; +arealims = [200]; +plotting = 1; + +[ica_segments, segmentlabel, segcentroid] = CellsortSegmentation(ica_filters, smwidth, thresh, arealims, plotting); + +%% 4b. CellsortApplyFilter +%filter signal and extract regions signals +cell_sig = CellsortApplyFilter(ica_segments, flims, , 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); +%% caluclate correlation matrix for cell_sig +clusterThreshold = 0.9; +[ica_segments, cluster, subPlotCols, subPlotWidth, sortedIdx] = calcCorrEff(cell_sig, ica_segments, clusterThreshold); +%% plot +%acustic sweep parameters: +%number of frames after trial begin when first sweep occurrs +sweepDelay = 100; +%number of frames between sweeps +sweepInterval = 30; +%number of sweeps +numberOfSweeps = 5; +%number of frames a sweep lasts +sweepDuration = 10; + +%number of time courses to show when slide bar is used +clusterSize = 10; +frameSteps = 20; + +cell_sig_merged = alternativePlot(ica_segments, cluster, cell_sig, subPlotCols, sortedIdx, data_r, frames, fileNamesGreen, clusterSize, sweepDelay, sweepInterval, numberOfSweeps, sweepDuration, frameSteps); + +%% 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 +cell_sig_merged = alternativePlot(ica_segments, cluster_new, cell_sig, subPlotCols, sortedIdx, data_r); diff --git a/main.m b/main.m index a29b2b5..b0f7613 100644 --- a/main.m +++ b/main.m @@ -86,10 +86,7 @@ %% 4b. CellsortApplyFilter %filter signal and extract regions signals -subtractmean = 1; - -cell_sig = CellsortApplyFilter(fn_green, ica_segments, flims, movm, subtractmean, data_g); - +cell_sig = CellsortApplyFilter(ica_segments, flims, data_g); %% 5. CellsortFindspikes %optionally, find spikes in signal deconvtau = 0;