main.m

%% 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 = 'D:\data\session_1';
outFolder = 'D:\data\out';
binningFactor = [1, 1, 2];
nSubIm = 10;
CellSortXCorr(inFolder, outFolder, nSubIm, binningFactor);
%% load files
%go to corr folder
outFolder = 'D:\data\out';
%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);
%% load files from alignment script
outFolder = 'D:\data\out';
[data_g, data_r, frames, fn_green, fn_red, fileNamesGreen, fileNamesRed] = CellSortLoadFilesFromAlignmentScript(outFolder);
%% 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,:);
%% bin data again
secondBinningFactor = [1,1,1];
data_g = binning(data_g, secondBinningFactor, 'mean');
data_r = binning(data_r, secondBinningFactor, 'mean');
%% 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);