Merge branch 'master' of https://github.com/ayalab1/neurocode

SharpWaveRipples/eventSpikingTreshold.m

@@ -40,7 +40,7 @@
 winSize = p.Results.winSize;
 eventSize = p.Results.eventSize;
 figOpt = p.Results.figOpt;
-EVENTFILE = p.Results.figOpt;
+EVENTFILE = p.Results.EVENTFILE;
 
 % 
 if isempty(spikes)

SharpWaveRipples/removeArtifactsFromEvents.m

@@ -63,13 +63,20 @@
 
 % 
 validEvents = find(stdEvents<=cutpoint);
-events.timestamps = events.timestamps(validEvents,:);
-try 
-    events.peaks = events.peaks(validEvents,:);
-    events.peakNormedPower = events.peakNormedPower(validEvents,:);
+fields = fieldnames(events);
+matchLength = size(events.timestamps,1);
+for i = 1:length(fields)
+    currentField = fields{i};
+    if size(events.(currentField),1)==matchLength
+        events.(currentField)= events.(currentField)(validEvents,:);
+    elseif size(events.(currentField),2)==matchLength
+        events.(currentField) = evetns.(currentField)(:,validEvents);
+    end
 end
+
 events.artifactsRemovalParameters.cutpoint = cutpoint;
 events.artifactsRemovalParameters.winSize = winSize;
+events.removed = ~ismember((1:length(stdEvents)),validEvents);
 
 fprintf('Keeping %4.0f of %4.0f events \n',length(validEvents),length(stdEvents));
 

behavior/general_behavior_file.m

@@ -124,7 +124,7 @@ function general_behavior_file(varargin)
 behavior.stateNames = stateNames;
 behavior.notes = notes;
 behavior.epochs = session.epochs;
-behavior.processinginfo.date = date;
+behavior.processinginfo.date = datetime("today");
 behavior.processinginfo.function = 'general_behavioral_file.mat';
 behavior.processinginfo.source = source;
 
@@ -141,14 +141,11 @@ function general_behavior_file(varargin)
         disp('you need to specify your environments in session')
         session = gui_session(session);
     end
-    start = [];
-    stop = [];
-    for ep = 1:length(session.epochs)
-        if ~contains(session.epochs{ep}.environment,'sleep')
-            start = [start,session.epochs{ep}.startTime];
-            stop = [stop,session.epochs{ep}.stopTime];
-        end
-    end
+    % load epochs and locate non-sleep epochs
+    epoch_df = load_epoch('basepath',basepath);
+    epoch_df = epoch_df(epoch_df.environment ~= "sleep",:);
+    start = epoch_df.startTime;
+    stop = epoch_df.stopTime;
 
     good_idx = manual_trackerjumps(behavior.timestamps,...
         behavior.position.x,...
@@ -191,7 +188,10 @@ function general_behavior_file(varargin)
                 end
                 convert_pix_to_cm_ratio = (pos_range / maze_sizes(maze_sizes_i));
                 maze_sizes_i = maze_sizes_i + 1;
-
+
+                % add convert_pix_to_cm_ratio to epochs in behavior file
+                behavior.epochs{ep}.pix_to_cm_ratio = convert_pix_to_cm_ratio;
+
                 % iterate over each tracker point
                 for pos_fields_i = 1:length(pos_fields)
                     behavior.position.(pos_fields{pos_fields_i})(idx) =...
@@ -215,7 +215,14 @@ function general_behavior_file(varargin)
             pos_range = maze_distance_gui(fullfile(files.folder,files.name));
         end
         convert_pix_to_cm_ratio = (pos_range / maze_sizes);
-
+
+        % add convert_pix_to_cm_ratio to epochs in behavior file
+        for ep = 1:length(session.epochs)
+            if ~contains(session.epochs{ep}.environment,'sleep')
+                behavior.epochs{ep}.pix_to_cm_ratio = convert_pix_to_cm_ratio;
+            end
+        end
+
         % iterate over each tracker point
         for pos_fields_i = 1:length(pos_fields)
             behavior.position.(pos_fields{pos_fields_i}) =...

behavior/manual_trackerjumps.m

@@ -1,9 +1,9 @@
 function good_idx = manual_trackerjumps(ts,x,y,StartofRec,EndofRec,basepath,varargin)
-% 
-% Manually cut out xy coordinates that are outside the bounds of your maze. 
-% 
-% These can be caused by unplugs or if the rat jumps out. 
-% If you do not remove these points, your ratemap will be messed up. 
+%
+% Manually cut out xy coordinates that are outside the bounds of your maze.
+%
+% These can be caused by unplugs or if the rat jumps out.
+% If you do not remove these points, your ratemap will be messed up.
 %
 % Input:
 %       ts
@@ -22,13 +22,15 @@
 addParameter(p,'axis_equal',true,@islogical);
 addParameter(p,'alpha',.2,@isnumeric);
 addParameter(p,'add_scatter',true,@islogical);
+addParameter(p,'save_boundary_file',false,@islogical);
 
 parse(p,varargin{:});
 darkmode = p.Results.darkmode;
 restic_dir = p.Results.restic_dir;
 axis_equal = p.Results.axis_equal;
 alpha = p.Results.alpha;
 add_scatter = p.Results.add_scatter;
+save_boundary_file = p.Results.save_boundary_file;
 
 savets=[];
 for i=1:length(StartofRec)
@@ -37,35 +39,40 @@
     ytemp=y(ts>=StartofRec(i) & ts<=EndofRec(i));
     tstemp=ts(ts>=StartofRec(i) & ts<=EndofRec(i));
     % use the gui to cut out points
-    [~,~,in]=restrictMovement(xtemp,ytemp,restic_dir,darkmode,axis_equal,alpha,add_scatter);
+    [~,~,in]=restrictMovement(xtemp,ytemp,restic_dir,darkmode,axis_equal,...
+        alpha,add_scatter);
     % save the ts where the tracker error exists
     savets=[savets,tstemp(in)];
 end
+
 % locate the index for each tracker error
 good_idx=ismember(ts,savets);
+
 % save that index to your session folder so you won't have to do this again
 % each time you run your data
-basename = basenameFromBasepath(basepath);
-save(fullfile(basepath,[basename,'.restrictxy.mat']),'good_idx')
+if save_boundary_file
+    basename = basenameFromBasepath(basepath);
+    save(fullfile(basepath,[basename,'.restrictxy.mat']),'good_idx')
+end
 end
 
 function [x,y,in]=restrictMovement(x,y,direction,darkmode,axis_equal,alpha,add_scatter)
 % restrictMovement allows you to draw a line around xy coordinates in order
 % to eliminate certain points you don't want...ie when the rat jumps out of
-% maze or tracker errors. 
+% maze or tracker errors.
 %
 % Also, I have included a direction input argument so you have either
 % restrict outside or inside points. This is valuable if you have a maze
 % like a circular track where the rat could jump out away or towards the
-% center of the maze. 
+% center of the maze.
 %
-% Input         x,y: coordinates 
+% Input         x,y: coordinates
 %         direction: 1 (default) to remove outside points; 0 to remove inside points
-%         
+%
 %
 % Output        x,y: retained coordinates
 %                in: logical of which coordinates were retained (so you can index ts etc.)
-%       
+%
 %
 % Ryan Harvey
 

behavior/process_and_sync_dlc.m

@@ -34,20 +34,20 @@
         if ~isempty(dir(fullfile(basepath,MergePoints.foldernames{ii},'*DLC*csv')))
             % locate file
             file = dir(fullfile(basepath,MergePoints.foldernames{ii},'*DLC*csv'));
-            
+
             % if there are more than 1 csv, then it is likely to be a dlc
             % filtered file. If so, take the filtered csv
             if length(file) > 1
                 file = dir(fullfile(basepath,MergePoints.foldernames{ii},'*DLC*filtered.csv'));
             end
-            
+
             video_file = dir(fullfile(file(1).folder,'*.avi'));
             obj = VideoReader(fullfile(video_file.folder,video_file.name));
             fs = obj.FrameRate;
-            
+
             % load csv with proper header
             df = load_dlc_csv(fullfile(file(1).folder,file(1).name));
-            
+
             % locate columns with [x,y,likelihood]
             field_names = df.Properties.VariableNames;
 
@@ -63,10 +63,10 @@
             ts = df{:,1}/fs;
             %             x = df{:,x_col(primary_coords)};
             %             y = df{:,y_col(primary_coords)};
-            
+
             x = df{:,x_col};
             y = df{:,y_col};
-            
+
             tempTracking{count} = sync_ttl(file.folder,x,y,ts,fs,pulses_delta_range);
             trackFolder(count) = ii;
             count = count + 1;
@@ -118,14 +118,14 @@
     % determine shape of data to insert
     n_col = size(tempTracking{ii}.position.x,2);
     n_row = size(tempTracking{ii}.position.x,1);
-    
+
     % insert data using above shape and offset
     x(row_offset+1:n_row+row_offset,1:n_col) = tempTracking{ii}.position.x;
     y(row_offset+1:n_row+row_offset,1:n_col) = tempTracking{ii}.position.y;
-    
+
     % change offset to reflect last inserted data
     row_offset = n_row + row_offset;
-    
+
     % metadata
     folder{ii} = tempTracking{ii}.folder;
     samplingRate = [samplingRate; tempTracking{ii}.samplingRate];
@@ -151,11 +151,26 @@
 
 if ~exist(fullfile(folder,'digitalIn.events.mat'),'file')
     if ~exist(fullfile(folder,'digitalin.dat'),'file')
-        error([fullfile(folder,'digitalin.dat'),'   does not exist'])
+
+        warning([fullfile(folder,'digitalin.dat'),'   does not exist'])
+
+        disp('decide if you want to use tracking without ttls (not recommended), if so type dbcont')
+        keyboard
+
+        [~,folder_name] = fileparts(folder);
+        tracking.position.x = x;
+        tracking.position.y = y;
+        tracking.timestamps = ts;
+        tracking.originalTimestamps = ts;
+        tracking.folder = folder_name;
+        tracking.samplingRate = fs;
+        tracking.description = '';
+        tracking.notes = 'no ttls';
+        return
     end
     digitalIn = getDigitalIn('all','folder',folder);
 end
-load(fullfile(folder,'digitalIn.events.mat'))
+load(fullfile(folder,'digitalIn.events.mat'),'digitalIn')
 
 Len = cellfun(@length, digitalIn.timestampsOn, 'UniformOutput', false);
 [~,idx] = max(cell2mat(Len));
@@ -166,9 +181,9 @@
 bazlerTtl(extra_pulses) = [];
 
 if isempty(bazlerTtl)
-   warning('no real ttls, something went wrong')
-   disp('fix the issue and type dbcont to continue')
-   keyboard 
+    warning('no real ttls, something went wrong')
+    disp('fix the issue and type dbcont to continue')
+    keyboard
 end
 
 basler_intan_diff = length(bazlerTtl) - size(x,1);
@@ -194,7 +209,7 @@
 if (length(bazlerTtl) == size(x,1)) || abs(basler_intan_diff)<=2 %assumes 1 frame could be cut at 0 and 1 frame at end
     notes = 'N of frames match!!';
     warning(notes);
-    
+
 elseif basler_intan_diff>0 && abs(basler_intan_diff)<fs
     notes = [num2str(abs(length(bazlerTtl) - size(x,1))) ' of frames dont match, probably at the end of the recording'];
     warning(notes);
@@ -210,15 +225,15 @@
     warning(notes);
     x = x(1:length(bazlerTtl),:);
     y = y(1:length(bazlerTtl),:);
-    
+
 elseif abs(basler_intan_diff)>60*fs
     notes = 'More than 1 minute missalignment in total in this session...will interpolate ts';
     warning(notes);
     simulated_ts = linspace(min(bazlerTtl),max(bazlerTtl),length(x));
     bazlerTtl_new = interp1(bazlerTtl,bazlerTtl,simulated_ts);
     bazlerTtl = bazlerTtl_new';
     fs = mode(diff(bazlerTtl));
-    
+
 elseif abs(basler_intan_diff)>2*fs
     notes = 'More than 2 seconds missalignment in total in this session...will adjust to the closer one...';
     warning(notes);
@@ -338,7 +353,7 @@
         % take timestamp in seconds
         digitalIn.timestampsOn{ii} = digital_on{ii}/fs;
         digitalIn.timestampsOff{ii} = digital_off{ii}/fs;
-        
+
         % intervals
         d = zeros(2,max([size(digitalIn.timestampsOn{ii},1) size(digitalIn.timestampsOff{ii},1)]));
         d(1,1:size(digitalIn.timestampsOn{ii},1)) = digitalIn.timestampsOn{ii};
@@ -349,7 +364,7 @@
         if d(2,end) == 0; d(2,end) = nan; end
         digitalIn.ints{ii} = d;
         digitalIn.dur{ii} = digitalIn.ints{ii}(2,:) - digitalIn.ints{ii}(1,:); % durantion
-        
+
         clear intsPeriods
         intsPeriods(1,1) = d(1,1); % find stimulation intervals
         intPeaks =find(diff(d(1,:))>lag);
@@ -366,13 +381,13 @@
     xt = linspace(0,size(data,2)/fs,size(data,2));
     data = flip(data);
     data = data(1:size(digitalIn.intsPeriods,2),:);
-    
+
     h=figure;
     imagesc(xt,1:size(data,2),data);
     xlabel('s'); ylabel('Channels'); colormap gray
     mkdir(fullfile(folder,'Pulses'));
     saveas(h,fullfile(folder,'Pulses','digitalIn.png'));
-    
+
     save(fullfile(folder,'digitalIn.events.mat'),'digitalIn');
 else
     digitalIn = [];

lfp/FindThetaCycles.m

@@ -7,7 +7,7 @@
 % as the troughs take theta asymmetry into account.
 %
 % It's recommended to provide lfp restricted to a behavior session (excluding
-% sleep sessions). The function will take into account theta asymmetry to find 
+% sleep sessions). The function will take into account theta asymmetry to find
 % the exact peak and trough timestamps as described by Belluscio et al (2012).
 %
 % USAGE
@@ -25,14 +25,14 @@
 %     'artefactThreshold' threshold to use to pass on to CleanLFP to exclude
 %                       artifacts when computing theta amplitude (default = 5)
 %     'baseline'        interval(s) of the behavior session, excluding sleep sessions,
-%                       (there is no need to restrict to running epochs) provided in [start stop] 
+%                       (there is no need to restrict to running epochs) provided in [start stop]
 %                       matrix format. These intervals will be used to estimate
 %                       the expected theta amplitude and compute amplitude thresholds.
 %                       (default = [0 Inf]);
 %     'marginAmplitude' the minimum theta amplitude (in sd-s). If the theta amplitude
 %                       during a cycle is lower than this amplitude, the cycle will
 %                       be discarded (default = -1). Note that this should be low
-%                       because during a behavioral epoch, theta oscillations are 
+%                       because during a behavioral epoch, theta oscillations are
 %                       expected (as the animal is running) in the majority of the
 %                       session.
 %    =========================================================================
@@ -107,11 +107,15 @@
 badsize = diff(peaktopeak,[],2) > maxDuration | diff(peaktopeak,[],2) < minDuration;
 ok(badsize) = false;
 % apply minimum amplitude threshold
-[~,bad,~] = CleanLFP(lfp,'thresholds',[artefactThreshold Inf],'manual',false); % the derivative threshold is set to Inf because theta is a slow signal and fast artefacts captured by the derivative are not relevant
+% the derivative threshold is set to Inf because theta is a slow signal
+% and fast artefacts captured by the derivative are not relevant
+[~,bad,~] = CleanLFP(lfp,'thresholds',[artefactThreshold Inf],'manual',false);
 amplitude(bad,2) = nan;
-nottheta = amplitude(~(amplitude(:,2)>amplitudeThreshold),1); 
-ok = CountInIntervals(nottheta, peaktopeak)==0; % intervals containing moments of low amplitude theta are not theta cycles
+nottheta = amplitude(~(amplitude(:,2)>amplitudeThreshold),1);
+% intervals containing moments of low amplitude theta are not theta cycles
+ok = CountInIntervals(nottheta, peaktopeak)==0 & ok;
 
 troughs = troughs(ok);
 peaktopeak = peaktopeak(ok,:);
 amplitude = amplitude(ok,2);
+end