Upload files to ''

MATLAB Scripts/CellReg_classification.m

@@ -0,0 +1,100 @@
+clear;
+
+[files, path] = uigetfile('*.mat', 'Choose CellReg cell_registered file');
+    load(fullfile(path,files));
+ 
+dir = uigetdir('Choose parent folder for first session');
+
+if(exist(fullfile(dir, 'classifications.mat'),'file'))
+        load(fullfile(dir, 'classifications.mat'))
+end
+if(exist(fullfile(dir, 'Cd.mat'),'file'))
+        load(fullfile(dir, 'Cd.mat'))
+end
+
+Cd = cell2mat(Cd);
+Cd_N = [];
+for i = 1:size(Cd,1)
+    Cd_N(i,:) = (Cd(i,:) - min(Cd(i,:))) / (max(Cd(i,:)) - min(Cd(i,:)));
+end
+FOV1 = Cd_N;
+one_classifications = classifications;
+
+dir = uigetdir('Choose parent folder for second session');
+
+if(exist(fullfile(dir, 'classifications.mat'),'file'))
+        load(fullfile(dir, 'classifications.mat'))
+end
+if(exist(fullfile(dir, 'Cd.mat'),'file'))
+        load(fullfile(dir, 'Cd.mat'))
+end
+
+Cd = cell2mat(Cd);
+for i = 1:size(Cd,1)
+    Cd_N(i,:) = (Cd(i,:) - min(Cd(i,:))) / (max(Cd(i,:)) - min(Cd(i,:)));
+end
+FOV2 = Cd_N;
+two_classifications = classifications;
+
+    
+for i = 1:length(cell_registered_struct.cell_to_index_map)
+    for j = 1:2
+        if cell_registered_struct.cell_to_index_map(i,j) == 0
+           cell_registered_struct.cell_to_index_map(i,:) = 0;
+        end
+    end
+end
+
+pairs = cell_registered_struct.cell_to_index_map;
+pairs(pairs < 1) = [];
+pairs = reshape(pairs, [length(pairs)/2, 2]);
+
+pre_class = pairs;
+name_class = char(15,2);
+
+for i = 1:size(pairs,1)
+    if any(one_classifications.active == pairs(i,1))
+        pre_class(i,1) = 1;
+    elseif any(one_classifications.quiescent == pairs(i,1))
+        pre_class(i,1) = 2;
+    else 
+        pre_class(i,1) = 3;
+    end
+    if any(two_classifications.active == pairs(i,2))
+        pre_class(i,2) = 1;
+    elseif any(two_classifications.quiescent == pairs(i,2))
+        pre_class(i,2) = 2;
+    else 
+        pre_class(i,2) = 3;
+    end
+end
+
+same_class = cell(size(pre_class,1),2);
+
+for i = 1:size(pre_class,1)
+    if pre_class(i,1) == 1
+        same_class{i,1} = 'active';
+    elseif pre_class(i,1) == 2
+        same_class{i,1} = 'quiescent';
+    else
+        same_class{i,1} = 'indiscriminate';
+    end
+    if pre_class(i,2) == 1
+        same_class{i,2} = 'active';
+    elseif pre_class(i,2) == 2
+        same_class{i,2} = 'quiescent';
+    else
+        same_class{i,2} = 'indiscriminate';
+    end
+end
+
+
+for i = 1:size(pairs,1)
+    FOV1_dff(i,:) = FOV1(pairs(i,1),:);
+    FOV2_dff(i,:) = FOV2(pairs(i,2),:);
+end
+
+save(fullfile(path,'FOV1_dff'),'FOV1_dff','-v7.3');
+save(fullfile(path,'FOV2_dff'),'FOV2_dff','-v7.3');
+save(fullfile(path,'same_class'),'same_class','-v7.3');
+       

MATLAB Scripts/free_move_correlation.m

@@ -0,0 +1,168 @@
+%% Freely Moving Kinetics 
+clear;
+%% Import file and extract full signal
+
+[files, path] = uigetfile({'*.MotoTrak;*.ArdyMotor'},...
+    'Select Mototrak File', ...
+    'multiselect','on');
+try 
+    cd(path);
+catch err
+    return
+end
+if ischar(files)                                                      %If only one file was selected...
+    files = {files};                                                    %Convert the string to a cell array.
+end
+for i = 1:length(files)
+    D = dir(files{1,i});
+    Bytes(i) = D.bytes;
+end
+[row, col] = find(Bytes>500);                             
+files = files(col);
+knob_data.trial_length = 500;
+knob_data.num_sessions = length(files);
+knob_data.combined_sessions_length = 0;
+knob_data.session_length = nan(1, length(files));
+knob_data.max_session_trials = 0;
+for s = 1:knob_data.num_sessions
+    try                                                                     %Try to read in the data file...
+        [~,filename,ext] = fileparts(files{s});
+        switch ext
+            case '.ArdyMotor'
+                data = ArdyMotorFileRead(files{s});
+            case '.MotoTrak'
+                data = MotoTrakFileRead(files{s});
+                data = MotoTrak_to_ArdyMotor(data);
+                data.rat = data.subject;
+        end
+        %         temp = ArdyMotorFileRead(files{f});                                 %Read in the data from each file.
+    catch err                                                               %If an error occurs...
+        warning(['ERROR READING: ' files{s}]);                              %Show which file had a read problem...
+        warning(err.message);                                               %Show the actual error message.
+        continue
+    end
+%     data = ArdyMotorFileRead(files{s});
+    temp = length(data.trial);
+    knob_data.combined_length = temp+knob_data.combined_sessions_length;
+    knob_data.session_length(s) = temp;
+    
+    if knob_data.session_length(s) > knob_data.max_session_trials
+        knob_data.max_session_trials = knob_data.session_length(s);
+    end
+end
+knob_data.combined_trials = nan(knob_data.combined_length, 500);
+knob_data.trial = nan(knob_data.max_session_trials, 500, knob_data.num_sessions);
+
+% for t = 1:length(data.trial)                                                %Step through each trial.
+%         data.trial(t).timestamps = ...
+%             double(data.trial(t).sample_times)/86400000 + ...
+%             data.trial(t).starttime;                                        %Convert the sample times to serial date numbers.
+% end
+% 
+% time = vertcat(data.trial.timestamps);
+for i = 1:length(files)
+    signal(:,i) = vertcat(data.trial.signal);
+end
+
+%% Find and classify peaks in the force signal
+
+lower = 15;
+upper = 18;
+three = 3;
+nine = 9;
+
+[pks, locs, widths] = findpeaks(signal, 'MinPeakHeight', three, 'MinPeakDistance', 8);  %Finds peaks and their corresponding locations for pulls over 3 grams at least 0.4s apart
+
+peakTimes = locs;
+pks = pks(1: length(peakTimes));
+NumPeaks = length(pks); 
+trial_start = 1:1:NumPeaks;
+trial_end = 1:1:NumPeaks;
+
+width_full = widths;
+
+avg_pull = mean(width_full);
+
+for i = 1:NumPeaks
+    trial_start(1,i) = peakTimes(i,1) - avg_pull;
+    trial_end(1,i) = peakTimes(i,1) + avg_pull;
+end
+
+R_trial_start = floor(trial_start);
+R_trial_end = ceil(trial_end);
+
+pullTimes = reshape([R_trial_start(:) R_trial_end(:)]', [1, 2* NumPeaks]); 
+
+behavior = cell(1, NumPeaks);
+
+for j = 1:2:length(pullTimes)
+    behavior{1,j} = signal(pullTimes(j) : pullTimes(j+1), 1);
+end
+
+behavior(cellfun('isempty',behavior)) = [];
+
+[pks, locs, widths] = findpeaks(signal, 'MinPeakHeight', lower, 'MinPeakDistance', 8);
+
+if pks > 0
+    for i = 1:length(pks)
+        if pks(i,:) > upper
+           locs(i,:) = 1;
+           widths(i,:) = 1;
+        end
+    end
+
+    locs(locs == 1) = [];
+    widths(widths == 1) = [];
+
+    successPeaks = locs;
+end
+
+[~, index_successPeaks, index_peakTimes] = intersect(successPeaks,peakTimes,'rows');
+successloc = zeros(1, length(peakTimes)); 
+
+successloc(:,index_peakTimes) = 1;
+successloc(1,1) = 0;
+  
+for j = 1:size(behavior,2) - 1
+    temp = corrcoef(behavior{1,j}, behavior{1,j+1});
+    pull_corr{1,j} = temp(1,2);
+end
+
+avg_pull_corr = mean(abs(cell2mat(pull_corr)));
+
+for i = 1:size(behavior,2) - 1
+    if successloc(:,i) == 1
+        temp = corrcoef(behavior{1,i}, behavior{1,i-1});
+        suc_prevcorr{1,i} = temp(1,2);
+    else
+        temp = NaN;
+        suc_prevcorr{1,i} = NaN;
+    end
+end
+
+for i = 1:size(behavior,2) - 1
+    if successloc(:,i) == 1
+        temp = corrcoef(behavior{1,i}, behavior{1,i+1});
+        suc_postcorr{1,i} = temp(1,2);
+    else
+        temp = NaN;
+        suc_postcorr{1,i} = NaN;
+    end
+end
+
+avg_presuccess_corr = nanmean(abs(cell2mat(suc_prevcorr)));
+avg_postsuccess_corr = nanmean(abs(cell2mat(suc_postcorr)));
+
+% %%
+% numNeurons = size(Cd_N,1);
+% activity = cell(numNeurons, length(pullTimes) / 2);
+% 
+% for i = 1:numNeurons
+%     for j = 1:2:length(pullTimes)
+%         activity{i,j} = Cd_N(i, pullTimes(j) : pullTimes(j+1));
+%     end
+% end
+% 
+% activity(cellfun('isempty',activity)) = [];
+% activity = reshape(activity, [numNeurons, length(meanFrames)/2]);
+% activity2 = cell2mat(activity);

MATLAB Scripts/processDFFInitVars.m

@@ -0,0 +1,113 @@
+function [ newNeurons, fluorescenceData, classifications, binaryPullTimes,pulls,options] = processDFFInitVars(dir, pullFrames, fr, autoClassifyNeurons, pTA)                                                                                                                                                                                                                                                                                                                                                                                                              
+% This function initializes variables needed for processDFFPipeline.m
+%   Detailed explanation goes here
+
+    %% Find the json file and Load in variables for dff extraction from same directory.
+    if isempty(dir)
+        disp('Pick Centroid json file');
+        [foldername, dir] = uigetfile('.json', 'Pick Centroid json file');     
+        jsonFilePath = fullfile(dir,foldername); % Set the file name using this variable
+    else
+        jsonFilePath = fullfile(dir, 'centroids.json');
+    end
+
+    if(exist(fullfile(dir, 'Fdf.mat'),'file'))
+        load(fullfile(dir, 'Fdf.mat'))
+    end
+    if(exist(fullfile(dir, 'Cd.mat'),'file'))
+        load(fullfile(dir, 'Cd.mat'))
+    end
+    if(exist(fullfile(dir, 'Sp.mat'),'file'))
+        load(fullfile(dir, 'Sp.mat'))
+    end
+
+    %% Create a new struct with neuron data (nid, dff, Cd, and Sp) concatenated by time (optional)
+    Fdf_concat = [];
+    Sp_concat = [];
+    Cd_concat = [];
+    for i = 1:length(F_df)
+        if(exist('F_df','var'))
+            Fdf_concat = horzcat(Fdf_concat, cell2mat(F_df(i)));
+        end
+        if(exist('Sp','var'))
+            Sp_concat = horzcat(Sp_concat, cell2mat(Sp(i)));
+        end
+        if(exist('Cd','var'))
+            Cd_concat = horzcat(Cd_concat, cell2mat(Cd(i)));
+        end
+    end
+    
+    fluorescenceData = struct('Fdf',Fdf_concat,'Cd',Cd_concat,'Sp',Sp_concat);
+    
+    newNeurons = struct('nid',[],'dff',[],'Cd',[],'Sp',[]);
+    neurons = jsonread(jsonFilePath);
+    numNeurons = length(neurons.jmesh);
+    for i = 1:numNeurons
+        newNeurons(i).nid = i;
+        %newNeurons(i).dff = Fdf_concat(i,:)';
+        newNeurons(i).Sp = Sp_concat(i,:)';
+        newNeurons(i).Cd = Cd_concat(i,:)';
+    end
+
+    %% Initialize Variables
+    numFrames = length(Cd_concat);
+    xpoints = 1:numFrames;
+
+    if isempty(pTA)
+        pTA = 1; % Default Value - frames before and after pull to average
+    end
+    if isempty(fr)
+        fr = 3; %% If no framerate set, just use frame numbers.
+    end
+    options = struct('numFrames',numFrames,'numNeurons',numNeurons,'pTA',pTA,'xpoints',xpoints,'framerate',fr);
+    
+    %% Pull Time Data
+
+    binaryPullTimes = zeros(1,numFrames);
+    for i = 1:2:length(pullFrames)
+        binaryPullTimes(pullFrames(i) - pTA :pullFrames(i+1) + pTA) = 1;
+    end
+    pulls= struct('pullNum',[],'pullFrames',[],'average',[]);
+    pullNum = 1;
+    for i = 1:2:length(pullFrames)
+        thisPull = Cd_concat(:,pullFrames(i) - pTA : pullFrames(i+1) + pTA);
+        meanPull = mean(thisPull,1);
+        pulls(pullNum).pullNum = pullNum;
+        pulls(pullNum).pullFrames = [pullFrames(i) pullFrames(i+1)];
+        pulls(pullNum).average = meanPull;
+        pullNum = pullNum + 1;
+    end
+    %% Initialize Data Frame for Classifying Cells as Active or Quiescent Active
+    % data(3).im(2).roi_trace_thresh(10,:) % Third Animal on second days 10th roi
+    % Data Struct - first input
+    data=struct('im',[]);
+    data.im = struct('roi_trace_thresh',Sp_concat,'roi_trace_df',Sp_concat);
+
+    % Analysis Struct - second input
+    % analysis(3).lever(2).lever_move_frames(:,1) % Third Animal on the Second
+    % day - binarized movement frames
+    analysis = struct('lever',[]);
+    analysis.lever = struct('lever_move_frames',[]);
+    analysis(1).lever(1).lever_move_frames = binaryPullTimes';
+    
+    [classified_rois, classified_p] = AP_classify_movement_cells_continuous(data,analysis); % Seems to work pretty well
+
+    %% Neuron Classification Variables
+    if isempty(autoClassifyNeurons)
+        autoClassifyNeurons = true;
+    end
+    if autoClassifyNeurons
+        active = find(classified_rois.movement);
+        quiesc = find(classified_rois.quiescent);
+        indisc = find(classified_rois.unclassified_active);
+    else % Have a csv file with the data for
+        disp('Pick neuron classification file');
+        nCfile = uigetfile(fullfile(dir,'*.csv'),'Pick neuron classification file');
+        neuronClass = csvread(fullfile(dir,nCfile));
+        active = find(neuronClass==1);
+        quiesc = find(neuronClass==2);
+        indisc = find(neuronClass==3);
+    end
+    classifications = struct('classified_rois',classified_rois,'classified_p',classified_p,'active',active,'quiescent',quiesc,'indisc',indisc);
+end
+

MATLAB Scripts/processDFF_Adaptive.m

@@ -0,0 +1,508 @@
+clear
+
+%% IMPORTANT NOTES
+
+% Voltage conversion equation: Force = m * ([Voltage] - b)
+% These constants will change if you use a different pull module or
+% recalibrate it (use Check_Pull_Calibration_Values script to recalculate)
+
+% paper position for saving .pdfs
+% x = -0.8, y = 1.0, width = 11.0, height = 7.1875
+
+%% Set Variables
+
+fr = 29.6;                                     %Framerate
+frames = 10000;                                %Number of Frames in the Session
+m = 25.045;                                    %Contants from Pull Calibration Script:
+b = 1.628;                                     %Force = m * ([Voltage] - b) 
+
+% m = 21.973, b = 0.049 (Newer values, but recalibrate for future experiments)
+%% Extract and Classify Pulls from Mototrak and Thorsync Files
+
+%Extract adaptive threshold from the selected MotoTrak file
+[file,path] = uigetfile('*MotoTrak', 'Select Mototrak file');
+cd(path);
+D = dir(file);
+[~,~,ext] = fileparts(file);
+data = MotoTrakFileRead(file);
+threshold = [data.trial.parameters];
+lower_thresh = threshold(:,2:2:end);
+volt_thresh = (lower_thresh / 25.045) + 1.628;
+
+%Load ThorSync file and find times the trials start 
+foldername = uigetdir('Choose folder to save results');
+[filename, pathname] = uigetfile('*.h5', 'Pick a Sync Episode file');
+if isequal(filename,0) || isequal(pathname,0)
+   disp('User pressed cancel')
+   return
+else
+   disp(['User selected ', fullfile(pathname, filename)])
+end
+
+fnam = [pathname,filename];
+
+dataOut = ThorsyncLoadSimple(fnam);
+
+time = dataOut.time;
+trial = dataOut.FrameIn; 
+
+time(trial(1,:) == 0) = [];
+
+Frame_In = uniquetol(time', 0.0001) - 0.05; %subtract 0.05s to account for the delay between MotoTrak and ThorSync
+
+%Find frames where the animal is pulling and classify them based on adaptive threshold
+end_frames = frames - 1;
+
+lower = (15/m) + b;
+upper = (18/m) + b;
+three = (3/m) + b;
+nine = (9/m) + b;
+
+%segment voltage data into cell by trial
+Pull_trial = cell(1, length(Frame_In));
+trial_time = cell(1, length(Frame_In));
+
+segment = floor(Frame_In * 100000);
+
+for i = 1:length(Frame_In) - 1
+    Pull_trial{1,i} = dataOut.Pull(segment(i):segment(i+1));
+    trial_time{1,i} = dataOut.time(segment(i):segment(i+1));
+    for j = length(Frame_In)
+        Pull_trial{1,j} = dataOut.Pull(segment(j):end);
+        trial_time{1,j} = dataOut.time(segment(j):end);
+    end
+end
+
+%find peaks in each block greater than 3 grams
+Peaks = cell(1, length(Frame_In));
+locations = cell(1, length(Frame_In));
+width = cell(1, length(Frame_In));
+
+
+for i = 1:length(Frame_In)
+    [pks, locs, widths] = findpeaks(Pull_trial{1,i}, trial_time{1,i}, 'MinPeakHeight', three, 'MinPeakDistance', 0.1 , 'MaxPeakWidth', 1);
+    Peaks{1,i} = pks;
+    locations{1,i} = locs;
+    width{1,i} = widths;
+end
+
+peakTimes = cell2mat(locations);
+pks = cell2mat(Peaks);
+peakFrames = cell2mat(locations) * fr;
+peakFrames(peakFrames > end_frames) = [];
+NumPeaks = length(peakFrames);
+trial_start = 1:1:NumPeaks;
+trial_end = 1:1:NumPeaks;
+widthFrames = cell2mat(width) * fr;
+avg_pull_frame = mean(widthFrames);
+avg_pull_frame = round(avg_pull_frame, 0);
+intFrames = round(peakFrames);
+intTimes = round((peakFrames/fr) * 100000);
+avg_pull_time = round((avg_pull_frame/fr) * 100000);
+
+for i = 1:NumPeaks
+    trial_start(1,i) = peakFrames(1,i) - widthFrames(1,i);
+    trial_end(1,i) = peakFrames(1,i) + widthFrames(1,i);
+end
+
+R_trial_start = floor(trial_start);
+R_trial_end = ceil(trial_end);
+
+
+%Find start and end times of pulls based on the width of the peak 
+pullFrames = reshape([R_trial_start(:) R_trial_end(:)]', [1, 2* NumPeaks]); 
+
+if pullFrames(1,1) < 3   %excludes first trial if there are an insufficient number of pull Frames for it
+   pullFrames(:,1) = [];  
+   pullFrames(:,1) = [];
+end
+
+pullForce = m * (pks - b);
+
+save(fullfile(foldername,'pullForce'),'pullForce','-v7.3');
+
+%excludes pulls that are cut off when ThorSync recording stops 
+pullFrames(pullFrames > end_frames) = [];
+if mod(length(pullFrames), 2) > 0
+   pullFrames(end) = [];
+end
+  
+[row, col] = find(pks > three & pks < nine);
+  threetonine = peakTimes(:,col);
+  threetonineFrames = threetonine * fr;
+  
+[row, col] = find(pks > nine & pks < lower);
+  ninetofifteen = peakTimes(:,col);
+  ninetofifteenFrames = ninetofifteen * fr;
+  
+[row, col] = find(pks > lower & pks < upper);
+  fifteentoeighteen = peakTimes(:,col);
+  fifteentoeighteenFrames = fifteentoeighteen * fr;
+  
+[row, col] = find(pks > upper);
+  eighteenplus = peakTimes(:,col);
+  eighteenplusFrames = eighteenplus * fr;
+%% Classify Successful and Unsuccessful Peaks based on Adaptive Threshold
+
+sPeaks = cell(1,length(Frame_In));
+uPeaks = cell(1,length(Frame_In));
+slocs = cell(1,length(Frame_In));
+ulocs = cell(1,length(Frame_In));
+swidths = cell(1,length(Frame_In));
+uwidths = cell(1,length(Frame_In));
+
+if length(Frame_In) < length(volt_thresh)
+    q = length(Frame_In);
+else
+    q = length(volt_thresh);
+end
+
+%if script throws index exceeds matrix error here replace length(Frame_In)
+%with the length of volt_thresh, make sure to change it back
+for i = 1:q
+    [col] = find(Peaks{1,i} > volt_thresh(i));
+    sPeaks{1,i} = Peaks{1,i}(col);
+    slocs{1,i} = locations{1,i}(col);
+    swidths{1,i} = width{1,i}(col);
+    [col] = find(Peaks{1,i} < volt_thresh(i));
+    uPeaks{1,i} = Peaks{1,i}(col);
+    ulocs{1,i} = locations{1,i}(col);
+    uwidths{1,i} = width{1,i}(col);
+end
+
+%find pull frames for each peak based on width and location
+successPeaks = cell2mat(slocs) * fr;
+successWidth = cell2mat(swidths) * fr;
+
+unsuccessPeaks = cell2mat(ulocs) * fr;
+unsuccessWidth = cell2mat(uwidths) * fr;
+
+trial_start = 1:1:length(successPeaks);
+trial_end = 1:1:length(successPeaks);
+
+for i = 1:length(successPeaks)
+    trial_start(1,i) = successPeaks(1,i) - successWidth(1,i);
+    trial_end(1,i) = successPeaks(1,i) + successWidth(1,i);
+end
+
+    R_trial_start = floor(trial_start);
+    R_trial_end = ceil(trial_end);
+
+    successFrames = reshape([R_trial_start(:) R_trial_end(:)]', [1, 2* length(successPeaks)]); 
+if successFrames > 0
+    if successFrames(1,1) < 3   %excludes first trial if there are an insufficient number of pull Frames for it
+       successFrames(:,1) = [];  
+       successFrames(:,1) = [];
+    end
+
+    %excludes pulls that are cut off when ThorSync recording stops 
+    successFrames(successFrames > frames) = [];
+    if mod(length(successFrames), 2) > 0
+       successFrames(end) = [];
+    end
+end
+    
+trial_start = 1:1:length(unsuccessPeaks);
+trial_end = 1:1:length(unsuccessPeaks);
+
+    for i = 1:length(unsuccessPeaks)
+        trial_start(1,i) = unsuccessPeaks(1,i) - unsuccessWidth(1,i);
+        trial_end(1,i) = unsuccessPeaks(1,i) + unsuccessWidth(1,i);
+    end
+
+    R_trial_start = floor(trial_start);
+    R_trial_end = ceil(trial_end);
+
+    unsuccessFrames = reshape([R_trial_start(:) R_trial_end(:)]', [1, 2* length(unsuccessPeaks)]); 
+
+    if unsuccessFrames(1,1) < 3   %excludes first trial if there are an insufficient number of pull Frames for it
+       unsuccessFrames(:,1) = [];  
+       unsuccessFrames(:,1) = [];
+    end
+
+    %excludes pulls that are cut off when ThorSync recording stops 
+    unsuccessFrames(unsuccessFrames > frames) = [];
+    if mod(length(unsuccessFrames), 2) > 0
+       unsuccessFrames(end) = [];
+    end
+    
+for i = 1:NumPeaks
+    mean_start(1,i) = intFrames(1,i) - avg_pull_frame;
+    mean_end(1,i) = intFrames(1,i) + avg_pull_frame;
+    time_start(1,i) = intTimes(1,i) - avg_pull_time;
+    time_end(1,i) = intTimes(1,i) + avg_pull_time;
+end
+
+meanFrames = reshape([mean_start(:) mean_end(:)]', [1, 2*NumPeaks]);
+meanTimes = reshape([time_start(:) time_end(:)]', [1, 2*NumPeaks]);
+
+if meanFrames(1,1) < 3
+    meanFrames(:,1) = [];
+    meanFrames(:,1) = [];
+end
+
+meanFrames(meanFrames > end_frames) = [];
+if mod(length(meanFrames), 2) > 0
+    meanFrames(end) = [];
+end
+
+if meanTimes(1,1) < (3/fr * 100000)
+    meanTimes(:,1) = [];
+    meanTimes(:,1) = [];
+end
+
+meanTimes(meanTimes > (end_frames/fr * 100000)) = [];
+if mod(length(meanTimes), 2) > 0
+    meanTimes(end) = [];
+end
+
+window = 5;                      %# of seconds in each trial
+Frame_In = round(Frame_In * 100000);
+Frame_Out = Frame_In + (window * 100000);
+
+Frame_times = reshape([Frame_In(:) Frame_Out(:)]', [1, 2*length(Frame_In)]);
+
+Frame_times(Frame_times > ((end_frames/fr)*100000)) = [];
+if mod(length(Frame_times), 2) > 0
+    Frame_times(end) = [];
+end
+%% Initialize variables for Correlations and Classification
+ 
+dir = '';
+autoClassifyNeurons = true;
+pTA = 2; % frames before and after pull that should be included in the average
+
+[newNeurons,fluorescenceData,classifications,binaryPullTimes,pulls,options] = processDFFInitVars(dir,pullFrames,fr,autoClassifyNeurons,pTA);
+
+% Unpack Data from function call
+numFrames = options.numFrames;
+numNeurons = options.numNeurons;
+xpoints = options.xpoints;
+framerate = options.framerate;
+%
+Fdf = fluorescenceData.Fdf;
+Cd = fluorescenceData.Cd;
+Sp = fluorescenceData.Sp;
+
+
+if pks > 0
+    binarySuccessTimes = zeros(1,numFrames);
+        for i = 1:2:length(successFrames)
+            binarySuccessTimes(successFrames(i) - pTA:successFrames(i+1) + pTA) = 1;
+        end
+end
+
+binarySuccessTimes(binarySuccessTimes > frames) = [];
+
+binaryUnsuccessTimes = zeros(1,numFrames);
+    for i = 1:2:length(unsuccessFrames)
+        binaryUnsuccessTimes(unsuccessFrames(i) - pTA:unsuccessFrames(i+1) + pTA) = 1;
+    end
+
+binaryUnsuccessTimes(binaryUnsuccessTimes > frames) = [];
+    
+Cd_N = [];
+for i = 1:numNeurons
+    Cd_N(i,:) = (Cd(i,:) - min(Cd(i,:))) / (max(Cd(i,:)) - min(Cd(i,:)));
+end
+
+fulldff = Cd_N';
+fullCd = Cd;
+
+activeCd = Cd_N(classifications.active,:);
+quiescentCd = Cd_N(classifications.quiescent,:);
+indiscCd = Cd_N(classifications.indisc,:);
+
+cells = cat(1, classifications.active, classifications.quiescent, classifications.indisc);
+cells = sortrows(cells, 'ascend');
+
+Cd_N = Cd_N(cells,:);
+Cd = Cd(cells,:);
+
+numNeurons = size(Cd_N,1);
+%activity_Cd = Cd_N + 1;
+
+behavior = cell(1, length(meanFrames) / 2);
+activity = cell(numNeurons, length(meanFrames) / 2);
+
+for i = 1:numNeurons
+    for j = 1:2:length(meanFrames)%j = 1:2:length(meanFrames)
+        behavior{1,j} = dataOut.Pull(1, meanTimes(j) : meanTimes(j+1));
+        activity{i,j} = Cd_N(i, pullFrames(j) : pullFrames(j+1));
+    end
+end
+
+activity(cellfun('isempty',activity)) = [];
+activity = reshape(activity, [numNeurons, length(meanFrames)/2]);
+activity_cat = cell2mat(activity);
+behavior(cellfun('isempty',behavior)) = [];
+
+%[coeff,score,latent,tsquared,explained,mu] = pca(activity_cat);
+%Preliminary pca analysis for activity
+
+for i = 1:numNeurons
+    for j = 1:size(behavior,2) - 1
+        temp = corrcoef(behavior{1,j}, behavior{1,j+1});
+        pull_corr{1,j} = temp(1,2);
+    end
+end
+
+for i = 1:numNeurons - 1
+    temp2 = corrcoef(activity_cat(i,:), activity_cat(i+1,:));
+    activity_corr(i,:) = temp2(1,2);
+end
+
+avg_pull_corr = mean(abs(cell2mat(pull_corr)));
+avg_act_corr = mean(activity_corr);
+Correlations = struct('pull_corr',avg_pull_corr, 'activity_corr', avg_act_corr);
+
+%% Preliminary Analysis of Peaks during Pulls vs Queiscent Periods
+
+flourFrames = cell(1,numNeurons);
+behaviorPeaks = zeros(1,numNeurons);
+nonBehaviorPeaks = zeros(1,numNeurons);
+SbehaviorPeaks = zeros(1,numNeurons);
+UbehaviorPeaks = zeros(1,numNeurons);
+allPeaks = zeros(1,numNeurons);
+
+for i = 1:numNeurons
+    ROI = Cd_N(i,:)';
+    [pks, locs, widths] = findpeaks(ROI, xpoints/framerate, 'MinPeakHeight', 0.1, 'MinPeakDistance', 0.35);
+    temp = locs * fr;
+    temp = locs * fr;
+    temp2 = locs * fr;
+    flourFrames{1,i} = round(temp,0);
+    for k = 1:length(locs)
+        if binaryPullTimes(1,flourFrames{1,i}(1,k)) > 0 
+            temp(1,k) = 1;
+        else
+            temp(1,k) = 0;
+        end
+    end
+    for k = 1:length(locs)
+        if binarySuccessTimes(1,flourFrames{1,i}(1,k)) > 0 
+            temp(1,k) = 1;
+        end
+    end
+    for k = 1:length(locs)
+        if binaryUnsuccessTimes(1,flourFrames{1,i}(1,k)) > 0 
+            temp2(1,k) = 1;
+        end
+    end
+    behaviorPeaks(1,i) = length(temp(temp == 1));
+    nonBehaviorPeaks(1,i) = length(temp(temp == 0));
+    SbehaviorPeaks(1,i) = length(temp(temp == 1));
+    UbehaviorPeaks(1,i) = length(temp2(temp2 == 1));
+    allPeaks(1,i) = length(temp);
+end
+
+peakPercent = behaviorPeaks ./ allPeaks;
+peakPercent(allPeaks < 5) = [];
+
+[plotPercent,index] = sortrows(peakPercent', 'descend'); 
+x = 1:length(peakPercent);
+% figure
+% scatter(x, plotPercent)
+
+behavior_ca = struct('behaviorPeaks', behaviorPeaks, 'nonBehaviorPeaks', nonBehaviorPeaks, 'SbehaviorPeaks', SbehaviorPeaks, 'UbehaviorPeaks', UbehaviorPeaks, 'allPeaks', allPeaks);
+
+% Choose which data we want to use for everything.
+dff = Cd_N'; 
+%Cd_N = normalized, Cd = not normalized
+
+figure;
+K = heatmap(Cd_N);
+K.GridVisible = 'off';
+K.Colormap = hot;
+
+%% Save files for segmentation
+
+save(fullfile(foldername,'Cd_N'),'Cd_N','-v7.3');
+save(fullfile(foldername,'classifications'),'classifications','-v7.3');
+save(fullfile(foldername,'activeCd'),'activeCd','-v7.3');
+save(fullfile(foldername,'quiescentCd'),'quiescentCd','-v7.3');
+save(fullfile(foldername,'indiscCd'),'indiscCd','-v7.3');
+save(fullfile(foldername,'peakFrames'),'peakFrames','-v7.3');
+save(fullfile(foldername,'successPeaks'),'successPeaks','-v7.3');
+save(fullfile(foldername,'unsuccessFrames'),'unsuccessPeaks','-v7.3');
+save(fullfile(foldername,'threetonineFrames'),'threetonineFrames','-v7.3');
+save(fullfile(foldername,'ninetofifteenFrames'),'ninetofifteenFrames','-v7.3');
+save(fullfile(foldername,'fifteentoeighteenFrames'),'fifteentoeighteenFrames','-v7.3');
+save(fullfile(foldername,'eighteenplusFrames'),'eighteenplusFrames','-v7.3');
+save(fullfile(foldername,'behavior_ca'),'behavior_ca','-v7.3');
+save(fullfile(foldername,'Correlations'),'Correlations','-v7.3');
+
+%% Start by looking at all neurons plotted on same plot and stacked 
+% Colors
+co = ...
+    [0        0.4470    0.7410;
+    0.8500    0.3250    0.0980;
+    0.9290    0.6940    0.1250;
+    0.4940    0.1840    0.5560;
+    0.4660    0.6740    0.1880;
+    0.3010    0.7450    0.9330;
+    0.6350    0.0780    0.1840];
+
+%% Plot shows all neuron plots stacked.
+
+figure;
+plot(xpoints/framerate, bsxfun(@plus,dff,0:numNeurons-1))
+hold  on
+
+p = patch(xpoints/framerate, binaryUnsuccessTimes * numNeurons, 'k');
+set(p, 'FaceAlpha', 0.2)
+set(p, 'EdgeAlpha', 0)
+
+if successpeaks > 0
+   m = patch(xpoints/framerate, binarySuccessTimes * numNeurons, 'r');
+   set(m, 'FaceAlpha', 0.4)
+   set(m, 'EdgeAlpha', 0)
+end
+
+nvs = gca;
+nvs.XTick = 0:50:max(xpoints/framerate);
+nvs.YTick = 0:10:numNeurons;
+if framerate ~= 1
+    xlabel('Time (seconds)');
+else
+    xlabel('Frame');
+end
+ylabel('Neuron ID');
+title('All Neurons Stacked');
+set(gca,'fontsize',24)
+set(gca,'LooseInset',get(gca,'TightInset'));
+
+saveas(gcf,strcat(foldername,'/dff','.fig'));
+
+
+%% Population, Active, Quiescent, Indiscriminant Averages
+active = classifications.active;
+quiesc = classifications.quiescent;
+indisc = classifications.indisc;
+
+figure;
+plot(xpoints/framerate, mean(fulldff,2)+3, 'col', co(1,:)) % Population Average
+hold on;
+plot(xpoints/framerate, mean(fulldff(:,active),2)+2, 'col', co(2,:)) % Active Average
+
+plot(xpoints/framerate, mean(fulldff(:,quiesc),2)+1, 'col', co(3,:)) % Quiescent Average
+
+plot(xpoints/framerate, mean(fulldff(:,indisc),2), 'col', co(4,:)) % Indiscriminant Active Average
+
+lgd = legend(['Population Average (n=' num2str(length(newNeurons)) ')'], ...
+    ['Active Average (n=' num2str(length(active)) ')'], ...
+    ['Quiescent Average (n=' num2str(length(quiesc)) ')'], ...
+    ['Indiscriminant Average (n=' num2str(length(indisc)) ')'], ...
+    'Location', 'eastoutside');
+% plot(xpoints/framerate, mean(dffActive,2)+2,'color', [0,0,0]+0.8)
+% plot(xpoints/framerate, mean(dffActive,2)+2, 'col', co(2,:)) % Active Average
+
+lgd.FontSize = 24;
+p = patch(xpoints/framerate, binaryPullTimes * 4, 'b');
+set(p, 'FaceAlpha', 0.4)
+set(p, 'EdgeAlpha', 0) % Pull Bars
+
+set(gca,'YTick',[])
+xlabel('Time (seconds)');
+set(gca,'fontsize',24)
+set(gca,'LooseInset',get(gca,'TightInset'));

MATLAB Scripts/processDFF_Isopull.m

@@ -0,0 +1,550 @@
+clear;
+
+%% IMPORTANT NOTES
+
+% Voltage conversion equation: Force = m * ([Voltage] - b)
+% These constants will change if you use a different pull module or
+% recalibrate it (use Check_Pull_Calibration_Values script to recalculate)
+
+% paper position for saving .pdfs
+% x = -0.8, y = 1.0, width = 11.0, height = 7.1875
+
+%% Set Variables
+
+fr = 29.6;                                     %Framerate
+frames = 10000;                                %Number of Frames in the Session
+m = 25.045;                                    %Contants from Pull Calibration Script:
+b = 1.628;                                     %Force = m * ([Voltage] - b) 
+
+% m = 21.973, b = 0.049 (Newer values, but recalibrate for future experiments)
+
+%% Select one h5 file:
+
+foldername = uigetdir('C:\','Choose folder to save results');
+[filename, pathname] = uigetfile('*.h5', 'Pick a Sync Episode file');
+if isequal(filename,0) || isequal(pathname,0)
+   disp('User pressed cancel')
+   return
+else
+   disp(['User selected ', fullfile(pathname, filename)])
+end
+
+fnam = [pathname,filename];
+
+dataOut = ThorsyncLoadSimple(fnam);
+time = dataOut.time;
+trial = dataOut.FrameIn; 
+
+time(trial(1,:) == 0) = [];
+
+Frame_In = uniquetol(time', 0.0001) - 0.05;
+pTA = 2;
+end_frames = frames - 1;
+
+lower = (15/m) + b;
+upper = (18/m) + b;
+three = (3/m) + b;
+nine = (9/m) + b;
+
+%% Extract and Classify Pulls from Thorsync File
+    
+[pks, locs, widths] = findpeaks(dataOut.Pull, dataOut.time, 'MinPeakHeight', three, 'MinPeakDistance', 0.4, 'MaxPeakWidth', 1);  %Finds peaks and their corresponding locations for pulls over 3 grams at least 0.4s apart
+peakTimes = locs;
+peakFrames = peakTimes * fr;
+pks = pks(1: length(peakTimes));
+NumPeaks = length(peakFrames); 
+trial_start = 1:1:NumPeaks;
+trial_end = 1:1:NumPeaks;
+widthFrames = (widths/2) * fr;
+widthFrames_full = widths * fr;
+avg_pull_frame = mean(widthFrames_full);
+avg_pull_frame = round(avg_pull_frame, 0);
+
+if avg_pull_frame == 0
+    avg_pull_frame = 1;
+end
+
+intFrames = round(peakFrames);
+intTimes = round((peakFrames/29.6) * 100000);
+avg_pull_time = round((avg_pull_frame/29.6) * 100000);
+
+if avg_pull_time == 0
+    avg_pull_time = 1;
+end
+
+
+for i = 1:NumPeaks
+    trial_start(1,i) = peakFrames(1,i) - widthFrames(1,i);
+    trial_end(1,i) = peakFrames(1,i) + widthFrames(1,i);
+end
+
+R_trial_start = floor(trial_start);
+R_trial_end = ceil(trial_end);
+
+
+%Find start and end times of pulls based on the width of the peak 
+pullFrames = reshape([R_trial_start(:) R_trial_end(:)]', [1, 2* NumPeaks]); 
+
+if pullFrames(1,1) < 3   %excludes first trial if there are an insufficient number of pull Frames for it
+   pullFrames(:,1) = [];  
+   pullFrames(:,1) = [];
+   pks(:,1) = [];
+end
+
+%excludes pulls that are cut off when ThorSync recording stops 
+pullFrames(pullFrames > end_frames - pTA) = [];
+if mod(length(pullFrames), 2) > 0
+   pullFrames(end) = [];
+   pks(end) = [];
+end
+
+pullForce = m * (pks - b);
+
+save(fullfile(foldername,'pullForce'),'pullForce','-v7.3');
+
+[row, col] = find(pks > lower & pks < upper);
+  temp = peakTimes;
+  temp(:,unique(col)) = [];
+  widthFrames(:,unique(col)) = [];
+  unsuccessWidth = widthFrames;
+  unsuccessPeaks = temp * fr;
+  
+[row, col] = find(pks > three & pks < nine);
+  threetonine = peakTimes(:,col);
+  threetonineFrames = threetonine * fr;
+  
+[row, col] = find(pks > nine & pks < lower);
+  ninetofifteen = peakTimes(:,col);
+  ninetofifteenFrames = ninetofifteen * fr;
+  
+[row, col] = find(pks > upper);
+  eighteenplus = peakTimes(:,col);
+  eighteenplusFrames = eighteenplus * fr;
+  
+%% Classify Successful and Unsuccessful Peaks based on Thorsync Voltage
+  
+trial_start = 1:1:length(unsuccessPeaks);
+trial_end = 1:1:length(unsuccessPeaks);
+
+for i = 1:length(unsuccessPeaks)
+    trial_start(1,i) = unsuccessPeaks(1,i) - unsuccessWidth(1,i);
+    trial_end(1,i) = unsuccessPeaks(1,i) + unsuccessWidth(1,i);
+end
+
+R_trial_start = floor(trial_start);
+R_trial_end = ceil(trial_end);
+
+unsuccessFrames = reshape([R_trial_start(:) R_trial_end(:)]', [1, 2* length(unsuccessPeaks)]); 
+
+if unsuccessFrames(1,1) < 3   %excludes first trial if there are an insufficient number of pull Frames for it
+   unsuccessFrames(:,1) = [];  
+   unsuccessFrames(:,1) = [];
+end
+
+%excludes pulls that are cut off when ThorSync recording stops 
+unsuccessFrames(unsuccessFrames > frames - pTA) = [];
+if mod(length(unsuccessFrames), 2) > 0
+   unsuccessFrames(end) = [];
+end
+
+[pks, locs, widths] = findpeaks(dataOut.Pull, dataOut.time, 'MinPeakHeight', lower , 'MinPeakDistance', 0.4, 'MaxPeakWidth', 1);
+
+
+if pks > 0
+    for i = 1:length(pks)
+        if pks(:,i) > upper
+           locs(:,i) = 1;
+           widths(:,i) = 1;
+        end
+    end
+
+    locs(locs == 1) = [];
+    widths(widths == 1) = [];
+if locs > 0
+    successPeaks = locs * fr;
+    successWidth = (widths/2) * fr;
+
+    trial_start = 1:1:length(locs);
+    trial_end = 1:1:length(locs);
+
+    for i = 1:length(locs)
+        trial_start(1,i) = successPeaks(1,i) - successWidth(1,i);
+        trial_end(1,i) = successPeaks(1,i) + successWidth(1,i);
+    end
+
+    R_trial_start = floor(trial_start);
+    R_trial_end = ceil(trial_end);
+
+    successFrames = reshape([R_trial_start(:) R_trial_end(:)]', [1, 2* length(locs)]); 
+
+    if successFrames(1,1) < 3   %excludes first trial if there are an insufficient number of pull Frames for it
+       successFrames(:,1) = [];  
+       successFrames(:,1) = [];
+    end
+
+    %excludes pulls that are cut off when ThorSync recording stops 
+    successFrames(successFrames > frames - pTA) = [];
+    if mod(length(successFrames), 2) > 0
+       successFrames(end) = [];
+    end
+    else
+    successPeaks = [];
+    successFrames = [];
+end
+else
+    successPeaks = [];
+    successFrames = [];
+end
+%outputs an array of pulls between 15 and 18 grams
+
+for i = 1:NumPeaks
+    mean_start(1,i) = intFrames(1,i) - avg_pull_frame;
+    mean_end(1,i) = intFrames(1,i) + avg_pull_frame;
+    time_start(1,i) = intTimes(1,i) - avg_pull_time;
+    time_end(1,i) = intTimes(1,i) + avg_pull_time;
+end
+
+meanFrames = reshape([mean_start(:) mean_end(:)]', [1, 2*NumPeaks]);
+meanTimes = reshape([time_start(:) time_end(:)]', [1, 2*NumPeaks]);
+
+if meanFrames(1,1) < 3
+    meanFrames(:,1) = [];
+    meanFrames(:,1) = [];
+end
+
+meanFrames(meanFrames > end_frames) = [];
+if mod(length(meanFrames), 2) > 0
+    meanFrames(end) = [];
+end
+
+if meanTimes(1,1) < (3/29.6 * 100000)
+    meanTimes(:,1) = [];
+    meanTimes(:,1) = [];
+end
+
+meanTimes(meanTimes > (end_frames/29.6 * 100000)) = [];
+if mod(length(meanTimes), 2) > 0
+    meanTimes(end) = [];
+end
+
+
+window = 5;                      %# of seconds in each trial
+Frame_In = round(Frame_In * 100000);
+Frame_Out = Frame_In + (window * 100000);
+
+Frame_times = reshape([Frame_In(:) Frame_Out(:)]', [1, 2*length(Frame_In)]);
+
+Frame_times(Frame_times > ((end_frames/29.6)*100000)) = [];
+if mod(length(Frame_times), 2) > 0
+    Frame_times(end) = [];
+end
+%% Initialize variables for Correlations and Classification
+ 
+dir = '';
+autoClassifyNeurons = true;
+pTA = 2; % frames before and after pull that should be included in the average
+
+[newNeurons,fluorescenceData,classifications,binaryPullTimes,pulls,options] = processDFFVars(dir,pullFrames,fr,autoClassifyNeurons,pTA);
+
+% Unpack Data from function call
+numFrames = options.numFrames;
+numNeurons = options.numNeurons;
+xpoints = options.xpoints;
+framerate = options.framerate;
+%
+Fdf = fluorescenceData.Fdf;
+Cd = fluorescenceData.Cd;
+Sp = fluorescenceData.Sp;
+
+if locs > 0
+    binarySuccessTimes = zeros(1,numFrames);
+        for i = 1:2:length(successFrames)
+            binarySuccessTimes(successFrames(i) - pTA:successFrames(i+1) + pTA) = 1;
+        end
+else
+    binarySuccessTimes = [];
+end
+
+binaryUnsuccessTimes = zeros(1,numFrames);
+    for i = 1:2:length(unsuccessFrames)
+        binaryUnsuccessTimes(unsuccessFrames(i) - pTA:unsuccessFrames(i+1) + pTA) = 1;
+    end
+
+Cd_N = [];
+for i = 1:numNeurons
+    Cd_N(i,:) = (Cd(i,:) - min(Cd(i,:))) / (max(Cd(i,:)) - min(Cd(i,:)));
+end
+
+fulldff = Cd_N';
+fullCd = Cd;
+
+activeCd = Cd_N(classifications.active,:);
+quiescentCd = Cd_N(classifications.quiescent,:);
+indiscCd = Cd_N(classifications.indisc,:);
+
+cells = cat(1, classifications.active, classifications.quiescent, classifications.indisc);
+cells = sortrows(cells, 'ascend');
+
+Cd_N = Cd_N(cells,:);
+Cd = Cd(cells,:);
+
+numNeurons = size(Cd_N,1);
+
+numPulls = length(peakFrames);
+peakFrames_2 = round(peakFrames');
+n = 10;
+peakFrames(peakFrames < n + 1) = [];
+
+start = 1:1:numPulls;
+finish = 1:1:numPulls;
+
+for i = 1:numPulls    % modify it to look at n frames before and after pull
+start(1,i) = peakFrames_2(i,1) - n;
+finish(1,i) = peakFrames_2(i,1) + n;
+end
+
+pullFrames_2 = reshape([start(:) finish(:)]', [1, 2* numPulls]);
+pullFrames_2(pullFrames_2 > length(Cd_N)) = [];
+if mod(length(pullFrames_2), 2) > 0
+   pullFrames_2(end) = [];
+end
+
+behavior = cell(1, length(meanFrames) / 2);
+activity = cell(numNeurons, length(meanFrames) / 2);
+
+for i = 1:numNeurons
+    for j = 1:2:length(meanFrames)%j = 1:2:length(Frame_times)
+        behavior{1,j} = dataOut.Pull(1, meanTimes(j) : meanTimes(j+1));
+        activity{i,j} = Cd_N(i, pullFrames_2(j) : pullFrames_2(j+1));
+    end
+end
+
+activity(cellfun('isempty',activity)) = [];
+activity = reshape(activity, [numNeurons, length(meanFrames)/2]);
+behavior(cellfun('isempty',behavior)) = [];
+activity_cat = cell2mat(activity);
+
+% [coeff,score,latent,tsquared,explained,mu] = pca(activity_cat);
+%Preliminary pca analysis for activity
+
+for i = 1:numNeurons
+    for j = 1:size(behavior,2) - 1
+        temp = corrcoef(behavior{1,j}, behavior{1,j+1});
+        pull_corr{1,j} = temp(1,2);
+    end
+end
+
+for i = 1:numNeurons - 1
+    temp2 = corrcoef(activity_cat(i,:), activity_cat(i+1,:));
+    activity_corr(i,:) = temp2(1,2);
+end
+
+avg_pull_corr = mean(abs(cell2mat(pull_corr)));
+avg_act_corr = mean(activity_corr);
+Correlations = struct('pull_corr',avg_pull_corr, 'activity_corr', avg_act_corr);
+
+%% Code for activity correlation between successful pulls and previous or following pull
+
+% [~, index_successPeaks, index_peakTimes] = intersect(successPeaks',peakFrames','rows');
+% successloc = zeros(1, length(peakTimes)); 
+% 
+% successloc(:,index_peakTimes) = 1;
+% successloc(1,1) = 0;
+% 
+% for i = 1:size(behavior,2) - 1
+%     if successloc(:,i) == 1
+%         temp = corrcoef(behavior{1,i}, behavior{1,i-1});
+%         suc_prevcorr{1,i} = temp(1,2);
+%     else
+%         temp = NaN;
+%         suc_prevcorr{1,i} = NaN;
+%     end
+% end
+% 
+% for i = 1:size(behavior,2) - 1
+%     if successloc(:,i) == 1
+%         temp = corrcoef(behavior{1,i}, behavior{1,i+1});
+%         suc_postcorr{1,i} = temp(1,2);
+%     else
+%         temp = NaN;
+%         suc_postcorr{1,i} = NaN;
+%     end
+% end
+% 
+% avg_presucc_corr = nanmean(abs(cell2mat(suc_prevcorr)));
+% avg_postsucc_corr = nanmean(abs(cell2mat(suc_postcorr)));
+% 
+% 
+% x1 = [1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]';
+% temp = repmat(x1, size(activity,2));
+% x1 = temp(:,1);
+% x2 = [0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0]';
+% temp = repmat(x2, size(activity,2));
+% x2 = temp(:,1);
+% x3 = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1]';
+% temp = repmat(x3, size(activity,2));
+% x3 = temp(:,1);
+% 
+% activity_cat = cell2mat(activity);
+% 
+% coefficients = cell(numNeurons,1);
+% for i = 1:numNeurons
+% x4 = activity_cat(i,:)';
+% tbl = table(x1,x2,x3,x4, 'VariableNames', {'pre', 'pull', 'post', 'activity'});
+% mdl = fitlm(tbl);
+% coefficients{i,1} = mdl.Coefficients;
+% %r_square{i,1} = mdl.Rsquared;
+% end
+
+%% Preliminary Analysis of Peaks during Pulls vs Queiscent Periods
+
+flourFrames = cell(1,numNeurons);
+behaviorPeaks = zeros(1,numNeurons);
+nonBehaviorPeaks = zeros(1,numNeurons);
+SbehaviorPeaks = zeros(1,numNeurons);
+UbehaviorPeaks = zeros(1,numNeurons);
+allPeaks = zeros(1,numNeurons);
+ 
+for i = 1:numNeurons
+    ROI = Cd_N(i,:)';
+    [pks, locs, widths] = findpeaks(ROI, xpoints/framerate, 'MinPeakHeight', 0.1);
+    temp = locs * fr;
+    temp2 = locs * fr;
+    temp2 = locs * fr;
+    flourFrames{1,i} = round(temp,0);
+    for k = 1:length(locs)
+        if binaryPullTimes(1,flourFrames{1,i}(1,k)) > 0 
+            temp(1,k) = 1;
+        else
+            temp(1,k) = 0;
+        end
+    end
+    for k = 1:length(locs)
+        if isempty(binarySuccessTimes)
+            temp2(1,k) = 0;
+        elseif binarySuccessTimes(1,flourFrames{1,i}(1,k)) > 0
+            temp2(1,k) = 1;
+        end
+    end
+    for k = 1:length(locs)
+        if binaryUnsuccessTimes(1,flourFrames{1,i}(1,k)) > 0 
+            temp2(1,k) = 1;
+        end
+    end
+    behaviorPeaks(1,i) = length(temp(temp == 1));
+    nonBehaviorPeaks(1,i) = length(temp(temp == 0));
+    SbehaviorPeaks(1,i) = length(temp2(temp2 == 1));
+    UbehaviorPeaks(1,i) = length(temp2(temp2 == 1));
+    allPeaks(1,i) = length(temp);
+end
+
+peakPercent = behaviorPeaks ./ allPeaks;
+peakPercent(allPeaks < 5) = [];
+
+[plotPercent,index] = sortrows(peakPercent', 'descend'); 
+x = 1:length(peakPercent);
+% figure
+% scatter(x, plotPercent)
+
+behavior_ca = struct('behaviorPeaks', behaviorPeaks, 'nonBehaviorPeaks', nonBehaviorPeaks, 'SbehaviorPeaks', SbehaviorPeaks, 'UbehaviorPeaks', UbehaviorPeaks, 'allPeaks', allPeaks);
+
+% Choose which data we want to use for everything.
+dff = Cd_N';
+%Cd_N = normalized, Cd = not normalized
+
+figure;
+K = heatmap(Cd);
+K.GridVisible = 'off';
+K.Colormap = hot;
+
+%% Save files for segmentation
+
+save(fullfile(foldername,'Cd_N'),'Cd_N','-v7.3');
+save(fullfile(foldername,'classifications'),'classifications','-v7.3');
+save(fullfile(foldername,'activeCd'),'activeCd','-v7.3');
+save(fullfile(foldername,'quiescentCd'),'quiescentCd','-v7.3');
+save(fullfile(foldername,'indiscCd'),'indiscCd','-v7.3');
+save(fullfile(foldername,'peakFrames'),'peakFrames','-v7.3');
+save(fullfile(foldername,'successPeaks'),'successPeaks','-v7.3');
+save(fullfile(foldername,'unsuccessFrames'),'unsuccessPeaks','-v7.3');
+save(fullfile(foldername,'threetonineFrames'),'threetonineFrames','-v7.3');
+save(fullfile(foldername,'ninetofifteenFrames'),'ninetofifteenFrames','-v7.3');
+save(fullfile(foldername,'eighteenplusFrames'),'eighteenplusFrames','-v7.3');
+save(fullfile(foldername,'behavior_ca'),'behavior_ca','-v7.3');
+save(fullfile(foldername,'Correlations'),'Correlations','-v7.3');
+
+%% Start by looking at all neurons plotted on same plot and stacked 
+% Colors
+co = ...
+    [0        0.4470    0.7410;
+    0.8500    0.3250    0.0980;
+    0.9290    0.6940    0.1250;
+    0.4940    0.1840    0.5560;
+    0.4660    0.6740    0.1880;
+    0.3010    0.7450    0.9330;
+    0.6350    0.0780    0.1840];
+
+%% Plot shows all neuron plots stacked.
+
+figure;
+plot(xpoints/framerate, bsxfun(@plus,dff,0:numNeurons-1))
+hold  on
+
+p = patch(xpoints/framerate, binaryUnsuccessTimes * numNeurons, 'k');
+set(p, 'FaceAlpha', 0.2)
+set(p, 'EdgeAlpha', 0)
+
+if successPeaks > 0 
+   m = patch(xpoints/framerate, binarySuccessTimes * numNeurons, 'r');
+   set(m, 'FaceAlpha', 0.4)
+   set(m, 'EdgeAlpha', 0)
+end
+
+nvs = gca;
+nvs.XTick = 0:50:max(xpoints/framerate);
+nvs.YTick = 0:10:numNeurons;
+if framerate ~= 1
+    xlabel('Time (seconds)');
+else
+    xlabel('Frame');
+end
+ylabel('Neuron ID');
+title('All Neurons Stacked');
+set(gca,'fontsize',24)
+set(gca,'LooseInset',get(gca,'TightInset'));
+
+saveas(gcf,strcat(foldername,'/dff','.fig')); 
+
+
+%% Population, Active, Quiescent, Indiscriminant Averages
+active = classifications.active;
+quiesc = classifications.quiescent;
+indisc = classifications.indisc;
+
+figure;
+plot(xpoints/framerate, mean(fulldff,2)+3, 'col', co(1,:)) % Population Average
+hold on;
+plot(xpoints/framerate, mean(fulldff(:,active),2)+2, 'col', co(2,:)) % Active Average
+
+plot(xpoints/framerate, mean(fulldff(:,quiesc),2)+1, 'col', co(3,:)) % Quiescent Average
+
+plot(xpoints/framerate, mean(fulldff(:,indisc),2), 'col', co(4,:)) % Indiscriminant Active Average
+
+lgd = legend(['Population Average (n=' num2str(length(newNeurons)) ')'], ...
+    ['Active Average (n=' num2str(length(active)) ')'], ...
+    ['Quiescent Average (n=' num2str(length(quiesc)) ')'], ...
+    ['Indiscriminant Average (n=' num2str(length(indisc)) ')'], ...
+    'Location', 'eastoutside');
+% plot(xpoints/framerate, mean(dffActive,2)+2,'color', [0,0,0]+0.8)
+% plot(xpoints/framerate, mean(dffActive,2)+2, 'col', co(2,:)) % Active Average
+
+lgd.FontSize = 24;
+p = patch(xpoints/framerate, binaryPullTimes * 4, 'b');
+set(p, 'FaceAlpha', 0.4)
+set(p, 'EdgeAlpha', 0) % Pull Bars
+
+set(gca,'YTick',[])
+xlabel('Time (seconds)');
+set(gca,'fontsize',24)
+set(gca,'LooseInset',get(gca,'TightInset'));
+