Improved delimiting central arm and the rim

calc_CellMetrics/getTrials_thetamaze.m

@@ -9,7 +9,7 @@
 %
 % maze : maza parameters
 %
-% plots : show summary plot? 
+% plots : show summary plot?
 %
 % circular_track - Added fields:
 %    .trials.alternation.start
@@ -19,18 +19,19 @@
 %    .trials.alternation.stateName
 %
 %    .states.left_right
-%    .states.error 
+%    .states.error
 %    .stateNames.left_right
 %    .stateNames.error
 
 if nargin < 3
     plots = 1;
 end
 if ~isfield(circular_track,'nSamples')
-	circular_track.nSamples = numel(circular_track.timestamps);
+    circular_track.nSamples = numel(circular_track.timestamps);
 end
 
-% Determining polar coordinates
+% Determining polar coordinates (angle in radians, distance from the
+% reference)
 [circular_track.position.polar_theta,circular_track.position.polar_rho] = cart2pol(circular_track.position.y,circular_track.position.x);
 
 % Changing from polar angle to position along circle by multiplying with the radius (inner radius)
@@ -39,53 +40,53 @@
 
 disp('Defining trials for the behavior')
 
-% Determining spatial limits 
+% Determining spatial limits
 
 % Onset of central arm
-central_arm_onset = find(diff(circular_track.position.y > maze.pos_y_limits(1))==1 ...
-    & circular_track.position.x(1:end-1) < maze.pos_x_limits(2)-5 ...
-    & circular_track.position.x(1:end-1) > maze.pos_x_limits(1)+5 ...
+central_arm_onset = find(diff(circular_track.position.y > maze.pos_y_limits(1))==1 ... % Delimited by x-delimiter (+/-0), lower y-delimiter, and negative y-coord
+    & circular_track.position.x(1:end-1) < maze.pos_x_limits(2) ...
+    & circular_track.position.x(1:end-1) > maze.pos_x_limits(1) ...
     & circular_track.position.y(1:end-1) < 0);
 
 % End of central arm
-central_arm_end = find(diff(circular_track.position.polar_rho > maze.pos_y_limits(2))==1 ...
-    & circular_track.position.x(1:end-1) < maze.pos_x_limits(2)-5 ...
-    & circular_track.position.x(1:end-1) > maze.pos_x_limits(1)+5 ...
+central_arm_end = find(diff(circular_track.position.polar_rho > maze.pos_y_limits(2))==1 ... % Delimited by x-delimiterm (+/-0), upper y-delimiter, and positive y-coord
+    & circular_track.position.x(1:end-1) < maze.pos_x_limits(2) ...
+    & circular_track.position.x(1:end-1) > maze.pos_x_limits(1) ...
     & circular_track.position.y(1:end-1) > 0);
 
 % Start of left arm
-left_rim_onset = find(diff(circular_track.position.x < maze.pos_x_limits(1)-5)==1 ...
-    & circular_track.position.y(1:end-1) > maze.pos_y_limits(2)-10);
+left_rim_onset = find(diff(circular_track.position.x < maze.pos_x_limits(1)-5)==1 ... % Delimited by upper left delimiter (-5/0)
+    & circular_track.position.y(1:end-1) > maze.pos_y_limits(2));
 
 % Start of right arm
-right_rim_onset = find(diff(circular_track.position.x > maze.pos_x_limits(2)+5)==1 ...
-    & circular_track.position.y(1:end-1) > maze.pos_y_limits(2)-10);
+right_rim_onset = find(diff(circular_track.position.x > maze.pos_x_limits(2)+5)==1 ... % Delimited by upper right delimiter (+5/0)
+    & circular_track.position.y(1:end-1) > maze.pos_y_limits(2));
 
 % End of left rim
-left_rim_end = find(diff(circular_track.position.polar_theta < -maze.polar_theta_limits(2))==1 ...
+left_rim_end = find(diff(circular_track.position.polar_theta < -maze.polar_theta_limits(2))==1 ... % Delimited by linearised track left delimiter, abs(x-coord) = 10, and distance delimiter -5
     & abs(circular_track.position.x(1:end-1)) > 10 ...
     & circular_track.position.polar_rho(1:end-1) > maze.polar_rho_limits(1)-5);
 
 % End of right rim
-right_rim_end = find(diff(circular_track.position.polar_theta > maze.polar_theta_limits(2))==1 ...
+right_rim_end = find(diff(circular_track.position.polar_theta > maze.polar_theta_limits(2))==1 ... % Delimited by linearised track right delimiter, abs(x-coord) = 10, and distance delimiter -5
     & abs(circular_track.position.x(1:end-1)) > 10 ...
     & circular_track.position.polar_rho(1:end-1) > maze.polar_rho_limits(1)-5);
 
-% All 
+% All
 pos7home = sort([left_rim_end,right_rim_end]);
 
-central_arm_end(find(diff(central_arm_end)<circular_track.sr)+1) = [];
+central_arm_end(find(diff(central_arm_end)<circular_track.sr)+1) = []; % Exclude trials shorter than a second (a precaution?)
 
 if plots == 1
     disp('Plotting position')
-    
+
     figure
-    
+
     plot(circular_track.position.x,circular_track.position.y,'color',[0.5 0.5 0.5]), hold on
     plot([-10,10],maze.pos_y_limits(1)*[1,1],'k') % Onset of central arm
     plot([-10,10],maze.pos_y_limits(2)*[1,1],'k') % End of central arm
-    plot([1,1]*maze.pos_x_limits(1)-5,maze.pos_y_limits(2)+[-10,20],'k') 
-    plot([1,1]*maze.pos_x_limits(2)+5,maze.pos_y_limits(2)+[-10,20],'k')
+    plot([1,1]*maze.pos_x_limits(1)-5,maze.pos_y_limits(2)+[0,30],'k')
+    plot([1,1]*maze.pos_x_limits(2)+5,maze.pos_y_limits(2)+[0,30],'k')
     [y1,x1] = pol2cart(maze.polar_theta_limits(2)/maze.radius_in,maze.polar_rho_limits(1)-5);
     [y2,x2] = pol2cart(maze.polar_theta_limits(2)/maze.radius_in,maze.polar_rho_limits(2));
     plot([x1,x2],[y1,y2],'-k')
@@ -97,7 +98,7 @@
     plot(circular_track.position.x(left_rim_end),circular_track.position.y(left_rim_end),'xc') % Left rim
     plot(circular_track.position.x(right_rim_end),circular_track.position.y(right_rim_end),'xy') % Right rim
     title('Position of the animal'), xlabel('X'), ylabel('Y'), zlabel('Z'),axis tight, % view(2)
-    if exist('plot_ThetaMaze.m')
+    if exist('plot_ThetaMaze.m','file')
         plot_ThetaMaze(maze)
     end
 end
@@ -116,29 +117,29 @@
 trials.start = 0; % Start time of
 trials.end = [];
 
-% Preparing trials matric
+% Preparing trials matrix
 trials.trials = nan(1,circular_track.nSamples);
-trials_states = zeros(1,circular_track.nSamples);
+trials.states = zeros(1,circular_track.nSamples);
 i = 0;
 
 % Behavioral scoring
 for j = 1:length(central_arm_end)
     test1 = find(central_arm_onset < central_arm_end(j));
     test2 = find(pos7home > central_arm_end(j));
-    if ~isempty(test2) & ~isempty(test1)
-        if (pos7home(test2(1))- central_arm_onset(test1(end)))/circular_track.sr < 50
+    if ~isempty(test2) && ~isempty(test1)
+        if (pos7home(test2(1))- central_arm_onset(test1(end)))/circular_track.sr < 50 % Does the whole track in less than 50 secs (corresponds to a single trial)
             if trials.start(end)-central_arm_onset(test1(end)) ~= 0
                 i = i+1;
                 trials.start(i) = central_arm_onset(test1(end));
                 trials.end(i) = pos7home(test2(1));
                 trials.trials(trials.start(i):trials.end(i)) = i;
-                trials_states(trials.start(i):trials.end(i)) = 1;
+                trials.states(trials.start(i):trials.end(i)) = 1;
                 if sum(ismember(left_rim_end,trials.end(i)))
                     % Left trial
                     circular_track.states.left_right(i) = 1;
                     if sum(ismember(right_rim_onset,trials.start(i):trials.end(i)))
                         circular_track.states.error(i) = true;
-                        else
+                    else
                         circular_track.states.error(i) = false;
                     end
                 elseif sum(ismember(right_rim_end,trials.end(i)))
@@ -158,9 +159,16 @@
 end
 
 % Changning from units of samples to units of time
-trials.start = circular_track.timestamps(trials.start);
-trials.end = circular_track.timestamps(trials.end);
-trials.nTrials = numel(trials.start);
+if isempty(trials.end)
+    trials.start = [];
+    trials.end = [];
+    trials.nTrials = 0;
+    plots = 0;
+else
+    trials.start = circular_track.timestamps(trials.start);
+    trials.end = circular_track.timestamps(trials.end);
+    trials.nTrials = numel(trials.start);
+end
 
 trials.stateName = 'Alternative running on track';
 
@@ -169,22 +177,22 @@
 
 if plots == 1
     figure,
-    
+
     subplot(1,2,1)
     plot(circular_track.timestamps-circular_track.timestamps(1),trials.trials,'.k','linewidth',2), xlabel('Time (sec)'), ylabel('Trials')
-    
+
     subplot(3,2,2)
-    stairs(circular_track.timestamps-circular_track.timestamps(1),trials_states,'.-k','linewidth',1), xlabel('Time (sec)'), ylabel('Trial')
-    
+    stairs(circular_track.timestamps-circular_track.timestamps(1),trials.states,'.-k','linewidth',1), xlabel('Time (sec)'), ylabel('Trial')
+
     subplot(3,2,4)
-    stairs(circular_track.states.left_right,'.-b','linewidth',1), xlabel('Trials'), ylabel('Left/Right'), 
+    stairs(circular_track.states.left_right,'.-b','linewidth',1), xlabel('Trials'), ylabel('Left/Right'),
     yticks(1:numel(circular_track.stateNames.left_right)), yticklabels(circular_track.stateNames.left_right)
-    
+
     subplot(3,2,6)
     stairs(circular_track.states.error,'.-k','linewidth',1), xlabel('Trials'), ylabel('Errors')
-    
+
     figure
-    
+
     subplot(1,2,1)
     plot(circular_track.position.x,circular_track.position.y,'.k','markersize',2), hold on
     idx_left = ismember(trials.trials, find(circular_track.states.left_right==1));
@@ -200,4 +208,4 @@
     idx_right = ismember(trials.trials, find(circular_track.states.left_right==2));
     plot3(circular_track.position.x(idx_right),circular_track.position.y(idx_right),circular_track.timestamps(idx_right),'.r','markersize',6)
     title('Trials and time'), xlabel('X'), ylabel('Y'), zlabel('Time (sec)')
-end
+end

calc_CellMetrics/linearize_theta_maze.m

@@ -12,12 +12,17 @@
 behavior.states.arm_rim = nan(1,length(behavior.timestamps)); % Numeric: 1, 2 or nan (outside)
 behavior.stateNames.arm_rim = {'Arm','Rim'};
 
-idx_arm = behavior.position.polar_rho < maze.polar_rho_limits(1) & behavior.position.y > maze.pos_y_limits(1) & behavior.position.y < maze.pos_y_limits(2);
+%idx_arm = behavior.position.polar_rho < maze.polar_rho_limits(1) & behavior.position.y > maze.pos_y_limits(1) & behavior.position.y < maze.pos_y_limits(2);
+idx_arm = behavior.position.polar_rho < maze.polar_rho_limits(1) & behavior.position.y > maze.pos_y_limits(1) & behavior.position.y < maze.pos_y_limits(2) ...
+  & behavior.position.x > -maze.radius_in/2 & behavior.position.x < maze.radius_in/2;
 behavior.states.arm_rim(idx_arm) = 1;
 
-idx_rim = behavior.position.polar_rho > maze.polar_rho_limits(1) & behavior.position.polar_rho < maze.polar_rho_limits(2)+10 & abs(behavior.position.polar_theta) < maze.polar_theta_limits(2);
+%idx_rim = behavior.position.polar_rho > maze.polar_rho_limits(1) & behavior.position.polar_rho < maze.polar_rho_limits(2)+10 & abs(behavior.position.polar_theta) < maze.polar_theta_limits(2);
+idx_rim = behavior.position.polar_rho > maze.radius_in/2 & behavior.position.polar_rho < maze.polar_rho_limits(2)+10 & abs(behavior.position.polar_theta) < maze.polar_theta_limits(2) & ~idx_arm;
 behavior.states.arm_rim(idx_rim) = 2;
 
+%behavior.position.polar_rho > maze.arm_half_width+((maze.radius_in-maze.arm_half_width)/2) ...
+
 % First the central arm is linearized
 pos_linearized = nan(size(behavior.timestamps));
 pos_linearized(behavior.states.arm_rim==1) = behavior.position.y(behavior.states.arm_rim==1)-maze.pos_y_limits(1);

calc_CellMetrics/loadOpenEphysSettingsFile.m

@@ -18,9 +18,9 @@
     session.extracellular.precision = 'int16';
     session.extracellular.fileName = '';
     session.extracellular.srLfp = openEphys_metadata.continuous(2).sample_rate;
-    session.extracellular.electrodeGroups.channels = 1:384;
+    session.extracellular.electrodeGroups.channels{1} = 1:384;
     session.extracellular.nElectrodeGroups = 1;
-    session.extracellular.spikeGroups.channels = 1:384;
+    session.extracellular.spikeGroups.channels{1} = 1:384;
     session.extracellular.nSpikeGroups = 1;
 
     session.timeSeries.dig.fileName = [openEphys_metadata.events(1).folder_name,'timestamps.npy'];

calc_CellMetrics/loadOptitrack.m

@@ -93,9 +93,20 @@
 optitrack_temp.TotalExportedFrames = str2double(dataArray{14}(1));
 optitrack_temp.RotationType = dataArray{16}(1);
 optitrack_temp.LenghtUnit = dataArray{18}(1);
-optitrack_temp.CoorinateSpace = dataArray{20}(1);
+optitrack_temp.CoordinateSpace = dataArray{20}(1);
 optitrack_temp.FrameRate = str2double(dataArray{6}{1});
 
+% Interpolate
+if sum(isnan(optitrack_temp.Xr))
+    optitrack_temp.Xr = interp1(optitrack_temp.Time(~isnan(optitrack_temp.Xr)),optitrack_temp.Xr(~isnan(optitrack_temp.Xr)),optitrack_temp.Time);
+    optitrack_temp.Yr = interp1(optitrack_temp.Time(~isnan(optitrack_temp.Yr)),optitrack_temp.Yr(~isnan(optitrack_temp.Yr)),optitrack_temp.Time);
+    optitrack_temp.Zr = interp1(optitrack_temp.Time(~isnan(optitrack_temp.Zr)),optitrack_temp.Zr(~isnan(optitrack_temp.Zr)),optitrack_temp.Time);
+    optitrack_temp.Wr = interp1(optitrack_temp.Time(~isnan(optitrack_temp.Wr)),optitrack_temp.Wr(~isnan(optitrack_temp.Wr)),optitrack_temp.Time);
+    optitrack_temp.X = interp1(optitrack_temp.Time(~isnan(optitrack_temp.X)),optitrack_temp.X(~isnan(optitrack_temp.X)),optitrack_temp.Time);
+    optitrack_temp.Y = interp1(optitrack_temp.Time(~isnan(optitrack_temp.Y)),optitrack_temp.Y(~isnan(optitrack_temp.Y)),optitrack_temp.Time);
+    optitrack_temp.Z = interp1(optitrack_temp.Time(~isnan(optitrack_temp.Z)),optitrack_temp.Z(~isnan(optitrack_temp.Z)),optitrack_temp.Time);
+end
+
 clear dataArray
 clearvars filename formatSpec fileID dataArray header_length;
 

calc_CellMetrics/loadSpikes.m

@@ -7,7 +7,8 @@
 %      KlustaViewa/Klustasuite
 %      MClust
 %      NWB
-%      Phy (default)
+%      CellExplorer (default)
+%      Phy (second default)
 %      Sebastien Royer's lab standard
 %      SpyKING Circus
 %      UltraMegaSort2000
@@ -70,7 +71,7 @@
 p = inputParser;
 addParameter(p,'basepath',pwd,@ischar); % basepath with dat file, used to extract the waveforms from the dat file
 addParameter(p,'clusteringpath',[],@ischar); % relativ clustering path to spike data (optional)
-addParameter(p,'format',[],@ischar); % clustering format: phy, klustakwik/neurosuite, KlustaViewa, NWB, Wave_clus, MClust, UltraMegaSort2000, ALF, AllenSDK
+addParameter(p,'format','CellExplorer',@ischar); % clustering format: phy, klustakwik/neurosuite, KlustaViewa, NWB, Wave_clus, MClust, UltraMegaSort2000, ALF, AllenSDK
                                                      % TODO: 'SpyKING CIRCUS', 'MountainSort', 'IronClust'
 addParameter(p,'basename','',@ischar); % The basename file naming convention
 addParameter(p,'electrodeGroups',nan,@isnumeric); % electrodeGroups: Loading only a subset of electrodeGroups from the spike format (only applicable to Klustakwik/neurosuite and KlustaViewa)
@@ -119,7 +120,7 @@
     basename = basenameFromBasepath(basepath);
 end
 
-if exist(fullfile(basepath,[basename,'.spikes.cellinfo.mat']),'file') && ~parameters.forceReload
+if strcmpi(format,'CellExplorer') && exist(fullfile(basepath,[basename,'.spikes.cellinfo.mat']),'file') && ~parameters.forceReload
     load(fullfile(basepath,[basename,'.spikes.cellinfo.mat']))
     if ~isfield(spikes,'processinginfo') || (isfield(spikes,'processinginfo') && spikes.processinginfo.version < 3 && strcmp(spikes.processinginfo.function,'loadSpikes') )
         parameters.forceReload = true;
@@ -129,6 +130,11 @@
     disp('loadSpikes: Using existing spikes file')
 % elseif exist(fullfile(basepath,[basename,'.spikes.cellinfo.mat']),'file') 
 %     load(fullfile(basepath,[basename,'.spikes.cellinfo.mat']))
+elseif strcmpi(format,'CellExplorer')
+    format = 'Phy';
+    parameters.forceReload = true;
+    spikes = [];
+    showGUI = true;
 else
     parameters.forceReload = true;
     spikes = [];
@@ -246,7 +252,7 @@
                 dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'HeaderLines' ,startRow-1, 'ReturnOnError', false);
                 fclose(fileID);
                 if isempty(dataArray{1})
-                    disp(['Noc clusters found in ', filename1,'. Will use the labels from KiloSort'])
+                    disp(['No clusters found in ', filename1,'. Will use the labels from KiloSort'])
                     filename = filename3;
                 else
                     filename = filename1;

calc_CellMetrics/mex/CCG.m

@@ -29,7 +29,7 @@
 %   ccg     [t x ngroups x ngroups] matrix where ccg(t,i,j) is the
 %           number (or rate) of events of group j at time lag t with  
 %           respect to reference events from group i
-%   t       time lag vector (units: seonds)
+%   t       time lag vector (units: seconds)
 %
 %  SEE
 %

calc_CellMetrics/preprocessOpenEphysData.m

@@ -57,10 +57,8 @@
         session.epochs{i}.startTime = startTime;
         if exist(fullfile(basepath,session.epochs{i}.name,'continuous','Neuropix-PXI-100.0','continuous.bin'),'file')
             inputFiles{i} = fullfile(basepath,session.epochs{i}.name,'continuous','Neuropix-PXI-100.0','continuous.bin');
-
         elseif exist(fullfile(basepath,session.epochs{i}.name,'continuous','Neuropix-PXI-100.0','continuous.dat'),'file')
             inputFiles{i} = fullfile(basepath,session.epochs{i}.name,'continuous','Neuropix-PXI-100.0','continuous.dat');
-
         else
             error(['Epoch duration could not be estimated as raw data file does not exist: ', inputFiles{i}]);
         end
@@ -108,13 +106,15 @@
 
 
         % 8. Merge digital timeseries
-        TTL_paths = {};
-        TTL_offsets = [];
-        for i = 1:numel(session.epochs)
-            TTL_paths{i} = fullfile(session.epochs{i}.name,'events','Neuropix-PXI-100.0','TTL_1');
-            TTL_offsets(i) = session.epochs{i}.startTime;
+        disp('Attempting to merge digital time series.')
+        digSeriesFilename = fullfile(basepath,'openephysDig.digitalseries.mat');
+        if exist(digSeriesFilename, 'file')
+            disp(['A merged digital time series file ' digSeriesFilename ...
+              'already exists. Skipping digital time series merger. Double check if the existing file is correct.'])
+        else
+            loadOpenEphysDigital(session);
+            disp('Successfully merged digital time series.')
         end
-        openephysDig = loadOpenEphysDigital(session,TTL_paths,TTL_offsets);
 
     end