From 226df8c54a021657c71b63f065cf1eb7d19b3d4f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jes=C3=BAs=20P=C3=A9rez-Ortega?= Date: Tue, 28 Jul 2020 15:42:35 -0400 Subject: [PATCH] Update --- .DS_Store | Bin 6148 -> 6148 bytes common/.DS_Store | Bin 0 -> 6148 bytes common/.MATLABDriveTag | 1 + common/Bin_Signal.m | 27 ++ common/Circle_Mask.m | 13 + common/Downsample_And_Crop_Stim.m | 12 + common/Find_Peaks_Or_Valleys.m | 284 ++++++++++++++++++ ...ind_Threshold_In_Cumulative_Distribution.m | 23 ++ .../Get_Indices_Before_And_After_Stimulus.m | 12 + common/Get_Peak_Or_Valley_Indices.m | 13 + common/Get_Raster_Correlation.m | 15 + common/Get_Raster_From_MEA_Units.m | 42 +++ common/Get_Raster_From_Time_Stamps.m | 51 ++++ common/Get_SEM.m | 8 + common/Hold_Axes.m | 8 + common/Hold_Figure.m | 16 + common/Hotelling_T2_Test.m | 113 +++++++ common/PlotR_JP.m | 52 ++++ common/Plot_Neurons_Measure.m | 18 ++ common/Plot_Raster.m | 84 ++++++ common/Plot_Raster_By_Stimulus.m | 20 ++ common/Plot_Raster_Motion_Stimulus.m | 58 ++++ common/Read_Colors.m | 77 +++++ common/Reshape_Raster.m | 21 ++ common/Save_Figure.m | 77 +++++ common/Scale.m | 36 +++ common/Set_Axes.m | 11 + common/Set_Figure.m | 36 +++ common/Set_Label_Time.m | 58 ++++ common/Shuffle_Raster_Test.m | 90 ++++++ common/Validate_Name.m | 21 ++ 31 files changed, 1297 insertions(+) create mode 100644 common/.DS_Store create mode 100644 common/.MATLABDriveTag create mode 100644 common/Bin_Signal.m create mode 100644 common/Circle_Mask.m create mode 100644 common/Downsample_And_Crop_Stim.m create mode 100644 common/Find_Peaks_Or_Valleys.m create mode 100644 common/Find_Threshold_In_Cumulative_Distribution.m create mode 100644 common/Get_Indices_Before_And_After_Stimulus.m create mode 100644 common/Get_Peak_Or_Valley_Indices.m create mode 100644 common/Get_Raster_Correlation.m create mode 100644 common/Get_Raster_From_MEA_Units.m create mode 100644 common/Get_Raster_From_Time_Stamps.m create mode 100644 common/Get_SEM.m create mode 100644 common/Hold_Axes.m create mode 100644 common/Hold_Figure.m create mode 100644 common/Hotelling_T2_Test.m create mode 100755 common/PlotR_JP.m create mode 100644 common/Plot_Neurons_Measure.m create mode 100644 common/Plot_Raster.m create mode 100644 common/Plot_Raster_By_Stimulus.m create mode 100644 common/Plot_Raster_Motion_Stimulus.m create mode 100644 common/Read_Colors.m create mode 100644 common/Reshape_Raster.m create mode 100644 common/Save_Figure.m create mode 100644 common/Scale.m create mode 100644 common/Set_Axes.m create mode 100644 common/Set_Figure.m create mode 100644 common/Set_Label_Time.m create mode 100644 common/Shuffle_Raster_Test.m create mode 100644 common/Validate_Name.m diff --git a/.DS_Store b/.DS_Store index e01012ecd916c37e9f816cb061e3eedb5e7980d6..0cecb9f629bcf1156ad48bd7c58bb4fffc36e940 100644 GIT binary patch literal 6148 zcmeHK-A)rh6h6c6mWnKY0veNzjfn|FA)?>~p%zpyAV>=q0n2VXlnvXRW_Qc4HGKtd zH734+4`AZ8@ddmxK7zi0iD!N)RK!amW=}KcJ9FloGdo|qGaUdRiBfI|KnDN~8^?wg zRDTdY&f77mQ=UddB5E*jtsTE;Im|j19Xdh=LI%zn1N_}>f-Ep7;HUQe{bUs#+Ybf@ zKa;xphQ_95P1B;$*63PtlkFq}Gbs59(_7)ac{A9wl$~PC`DL5dl1+Cqde;PV1Jib? zCz9H7dHKaFNh`r!NJyK)D!fLcW|ek!TkTThLv@v3b-lXINK5ow>-jdl zy`!_etG)X|cTaD3-^ELP8NJgovs-zOKk)=}n3r%$bS3ZExn$AI+HO7_(A^+u@6(Ll z@muHu%{@DGrS^&#({<<&%qdalUB{3Yu7h{y8NYmRAAXUPZN7-JIh@81oU>{z>8+Zrr z;Uj#3Z}5XO6P;Wk{bZPok#RCb9+L%{m6(p;hzEZYMq0VQwD% zJSPfM_^Wmx2N|4y7by5}9SX(3nko5#hY?p`5q%Uq)HmRvh$snr5`7GcV4>c{v01E; zf*I5a>>YV$u=lr_9;um5coxh(9Z{TW(bJheoHb-1WZ(=l!1o6O8^^f9ibSz>pb}RA z#0GQ=L7U$tc#fekuCO8zM^K1NMU<&TT``DEN57%s;tDGgWjYXb`5TEr>IDXa|*w(Re r;c*p-A_bK?j)lRG;$Cb*aLnKaVq9TGB6?85KLUb=P=yTqQwDwlEEw?3 delta 97 zcmZoMXfc=|#>B`mu~2NHo}wrd0|Nsi1A_nqLp*~cLj^-BLotKi#KPs14MbQrzh;wR t+^oRC&a$!LD&uB$4t@@x#?68p-H1@V-^m;4Wg<&0T*E43hX&L&p$$qDprKhvt+--jT7}7np#A3 zem<@ulZcFPQ@L2!n>{z**++&mCkOWA81W14cNZlEfg7;MkzE(HCqgga^y>{tEnwC%0;vJ&^%eQ zLs35+`xjp>T00)'; + elseif step<0 + binned = Find_Peaks_Or_Valleys(signal,0.5,true,true,0,step,true)>0; + binned = (Find_Peaks_Or_Valleys(double(binned),0.5,true,true,0,1,true)>0)'; + else + binned = Find_Peaks_Or_Valleys(signal,0.5,true,true,0,step,true)>0; + binned = (Find_Peaks_Or_Valleys(double(binned),0.5,true,false,0,-1,true)>0)'; + end + case 'after_onset' + binned = (Find_Peaks_Or_Valleys(signal,0.5,true,true,0,step,true)>0)'; + case 'after' + binned = (Find_Peaks_Or_Valleys(signal,0.5,true,false,0,step,true)>0)'; + case 'before' + binned = (Find_Peaks_Or_Valleys(signal,0.5,true,true,0,-step,true)>0)'; +end \ No newline at end of file diff --git a/common/Circle_Mask.m b/common/Circle_Mask.m new file mode 100644 index 0000000..129e0e0 --- /dev/null +++ b/common/Circle_Mask.m @@ -0,0 +1,13 @@ +function mask = Circle_Mask(image_size,xy_center,radius) +% Create a circle mask +% +% Jesus Perez-Ortega, jesus.epo@gmail.com +% March, 2019 + +% Get size of image +h = image_size(1); +w = image_size(2); +[x,y] = meshgrid(1:w,1:h); + +% Get mask +mask = sqrt((x-xy_center(1)).^2 + (y-xy_center(2)).^2) < radius; \ No newline at end of file diff --git a/common/Downsample_And_Crop_Stim.m b/common/Downsample_And_Crop_Stim.m new file mode 100644 index 0000000..d7e760c --- /dev/null +++ b/common/Downsample_And_Crop_Stim.m @@ -0,0 +1,12 @@ +function stim = Downsample_And_Crop_Stim(stimuli,period_ms,samples) +% Get the visual stimulus of the same length of the imaging +% +% The signal should be recorded in ms, and period should be a integer +% number of ms. +% +% stim = Downsample_And_Crop_Stim(stimuli,period_ms,samples) +% +% By Jesus Perez-Ortega, July 2019 + +stim = downsample(round(stimuli*2),period_ms); +stim = stim(1:samples); \ No newline at end of file diff --git a/common/Find_Peaks_Or_Valleys.m b/common/Find_Peaks_Or_Valleys.m new file mode 100644 index 0000000..10c11fc --- /dev/null +++ b/common/Find_Peaks_Or_Valleys.m @@ -0,0 +1,284 @@ +function [indices,widths,amplitudes,ini_fin] = Find_Peaks_Or_Valleys(data,threshold,join,... +detect_peaks,minimum_width,fixed_width,ignore_ini_fin) +% Find peaks from signal by a given threshold. +% +% [indices,widths,amplitudes,ini_fin] = Find_Peaks_Or_Valleys(data,threshold,join,... +% detect_peaks,minimum_width,fixed_width,ignore_ini_fin) +% +% default: threshold = 0; join = true; detect_peaks = true; +% minimum_width = 0; fixed_width = 0; ignore_ini_fin = false; +% +% Inputs +% data = data as vector Fx1 (F = #frames) +% threshold = threshold +% join = set mode to get peaks (0 = each vector above threshold is a peak; +% 1 = joining of adjacent vectors above the threshold is a peak) +% +% Outputs +% indices = Fx1 vector containing the peak indices +% +% 1. Get indices above threshold +% 2. Ignore initial and final peaks (or valleys) - Optional +% 3. Restriction to minimum width - Optional +% 4. Set a fixed width - Optional +% 5. Set the number at each peak (or valley) +% 6. Join peaks or valleys - Optional +% +% by Jesus E. Perez-Ortega, Feb-2012 +% last modification July-2019 + +switch nargin + case 6 + ignore_ini_fin = false; + case 5 + ignore_ini_fin = false; + fixed_width = 0; + case 4 + ignore_ini_fin = false; + fixed_width = 0; + minimum_width = 0; + case 3 + ignore_ini_fin = false; + fixed_width = 0; + minimum_width = 0; + detect_peaks = true; + case 2 + ignore_ini_fin = false; + fixed_width = 0; + minimum_width = 0; + detect_peaks = true; + join = true; + case 1 + ignore_ini_fin = false; + fixed_width = 0; + minimum_width = 0; + detect_peaks = true; + join = true; + threshold = 0; +end + +% 0. Correct signal data +if(size(data,1)==1) + data = data'; +end +original_data = data; + +% 1. Get peak or valley indices +[idx,count] = Get_Peak_Or_Valley_Indices(data,threshold,detect_peaks); +% Size of data +F = numel(data); +indices=zeros(F,1); +if ~count + if(detect_peaks) + disp('No peaks found!') + widths = []; + amplitudes = []; + ini_fin = []; + else + disp('No valleys found!') + widths = []; + amplitudes = []; + ini_fin = []; + end + return +end + +% 2. Ignore initial and final peak +if ignore_ini_fin + + % Delete if start above threshold + last=1; + idx=idx'; + for i=idx + if(last==i) + if detect_peaks + data(i)=threshold-1; + else + data(i)=threshold+1; + end + last=last+1; + else + break; + end + end + + % Delete if ends above threshold + last = F; + idx = fliplr(idx); + for i = idx + if last==i + if detect_peaks + data(i)=threshold-1; + else + data(i)=threshold+1; + end + last=last-1; + else + break; + end + end + + % Get peak or valley indices + [idx,count] = Get_Peak_Or_Valley_Indices(data,threshold,detect_peaks); + if ~count + if detect_peaks + disp('No peaks found!') + else + disp('No valleys found!') + end + return + end +end + +% 3. Minimum width (after join peaks or valleys) +if minimum_width + + % Join peaks or valleys + is = find(idx~=[0; idx(1:numel(idx)-1)+1]); % index of same peak + % number of total peaks or valleys + count = numel(is); + if count + for j = 1:count-1 + indices(idx(is(j)):idx(is(j+1)-1),1)=j; % set #peak + end + indices(idx(is(count)):max(idx),1)=count; + end + + % Get peaks or valleys width + widths=[]; + for i=1:count + widths(i)=length(find(indices==i)); + end + + % Evaluate peaks less than or equal to minimum width + idx_eval=find(widths<=minimum_width); + widths=widths(idx_eval); + + % number of peaks to eliminate + count_less=length(widths); + + % Detect initial and final times + if count_less>0 + for i=1:count_less + peak=find(indices==idx_eval(i)); + ini_peak=peak(1); + end_peak=peak(end); + if(detect_peaks) + data(ini_peak:end_peak)=threshold-1; + else + data(ini_peak:end_peak)=threshold+1; + end + end + end + + % Get peak or valley indices + [idx,count] = Get_Peak_Or_Valley_Indices(data,threshold,detect_peaks); + if ~count + if(detect_peaks) + disp('No peaks found!') + else + disp('No valleys found!') + end + return + end +end + +% 4. Set fixed width +if fixed_width + last_end = -1; + end_before = false; + for i = idx' + if i==last_end+1 + if detect_peaks + data(i) = threshold-1; + else + data(i) = threshold+1; + end + if end_before + end_before = false; + else + last_end = i; + end + else + if i>last_end + if fixed_width<0 + ini = i+fixed_width; + fin = i-1; + if ini<1 + ini = 1; + end + fixed_width_peak = ini:fin; + last_end = fin+1; + else + fin = i+fixed_width-1; + if(fin>F) + fin = F; + end + fixed_width_peak = i:fin; + last_end = fin; + end + + if detect_peaks + data(fixed_width_peak) = threshold+1; + if fixed_width<0 + fixed_width_peak = fixed_width_peak+length(fixed_width_peak); + data(fixed_width_peak) = threshold-1; + elseif sum(data(fixed_width_peak)threshold) + end_before = true; + end + end + end + end + end + + % Get peak or valley indices + [idx,count] = Get_Peak_Or_Valley_Indices(data,threshold,detect_peaks); +end + +% 5. Put numbers to peaks +indices=zeros(F,1); +for i=1:count + indices(idx(i))=i; +end + +% 6. Join peaks or valleys +if join + is = find(idx~=[0; idx(1:numel(idx)-1)+1]); % index of same peak + + % number of total peaks + count = numel(is); + if count + for j = 1:count-1 + indices(idx(is(j)):idx(is(j+1)-1),1)=j; % set #peak + end + indices(idx(is(count)):max(idx),1)=count; + end +end + +% Get peaks or valleys width +widths = zeros(count,1); +ini_fin = zeros(count,2); +for i = 1:count + id = find(indices==i); + ini_fin(i,1) = id(1); + ini_fin(i,2) = id(end); + widths(i) = length(id); +end + +% Get peaks or valleys amplitud +amplitudes = zeros(count,1); +for i = 1:count + if detect_peaks + value = max(original_data(indices==i)); + else + value = min(original_data(indices==i)); + end + amplitudes(i) = value; +end \ No newline at end of file diff --git a/common/Find_Threshold_In_Cumulative_Distribution.m b/common/Find_Threshold_In_Cumulative_Distribution.m new file mode 100644 index 0000000..3c6a6d4 --- /dev/null +++ b/common/Find_Threshold_In_Cumulative_Distribution.m @@ -0,0 +1,23 @@ +% Find the threshold upper a given alpha +% +% Jesús Pérez-Ortega nov 2018 + +function th = Find_Threshold_In_Cumulative_Distribution(data,alpha) + + if(min(data)==max(data)) + th = max(data)+1; + else + if(max(data)<=1) + x = 0:0.001:1; + elseif(max(data)<=10) + x = 0:0.01:10; + else + x = 0:max(data); + end + y = hist(data,x); + cdy = cumsum(y); + cdy = cdy/max(cdy); + id = find(cdy>(1-alpha),1,'first'); + th = x(id); + end +end \ No newline at end of file diff --git a/common/Get_Indices_Before_And_After_Stimulus.m b/common/Get_Indices_Before_And_After_Stimulus.m new file mode 100644 index 0000000..45542a9 --- /dev/null +++ b/common/Get_Indices_Before_And_After_Stimulus.m @@ -0,0 +1,12 @@ +function indices = Get_Indices_Before_And_After_Stimulus(stimulus,before_samples,after_samples) +% Get the indices befor and after and stimulus +% +% indices = Get_Indices_Before_And_After_Stimulus(stimulus,before_samples,after_samples) +% +% Jesus Perez-Ortega May 2019 + +during = stimulus>0; +before = Bin_Signal(during,before_samples,'before'); +after = Bin_Signal(during,after_samples,'after'); +indices = find(before | during' | after); + \ No newline at end of file diff --git a/common/Get_Peak_Or_Valley_Indices.m b/common/Get_Peak_Or_Valley_Indices.m new file mode 100644 index 0000000..8116d45 --- /dev/null +++ b/common/Get_Peak_Or_Valley_Indices.m @@ -0,0 +1,13 @@ +function [indices,count] = Get_Peak_Or_Valley_Indices(data,threshold,detect_peaks) +% Get peak or valley indices +% +% by Jesus E. Perez-Ortega - Feb 2012 +% Modified June 2019 + +if detect_peaks + indices = find(data>threshold); +else + % detect valleys + indices = find(datamax_time + max_time = max_i; + end +end +max_time = ceil(max_time)*1000; + +% Set raster size +raster = zeros(1,max_time,'logical'); + +% Build raster +k=1; +for i = 1:n_units + time_stamps_unit = evalin('base',[units_names{i} ';']); + spikes_id = round(time_stamps_unit*1000/window_ms)+1; + if ~isnan(spikes_id) + raster(k,spikes_id) = true; + end + k=k+1; +end \ No newline at end of file diff --git a/common/Get_Raster_From_Time_Stamps.m b/common/Get_Raster_From_Time_Stamps.m new file mode 100644 index 0000000..deda783 --- /dev/null +++ b/common/Get_Raster_From_Time_Stamps.m @@ -0,0 +1,51 @@ +function raster = Get_Raster_From_Time_Stamps(window,initial_time,final_time) +% Convert raw data (time stamps in seconds) to raster +% +% Time stamps have to be in the workspace, and window defines the bin of the raster in s. +% +% raster = Get_Raster_From_Time_Stamps(window,initial_time,final_time) +% +% By Jesús Pérez-Ortega July-2018 +% modified June 2019 + +switch nargin + case 1 + partial = false; + case 3 + partial = true; +end + +unit_names = evalin('base','who'); +n_cells = length(unit_names); + +if partial + max_time=ceil((final_time-initial_time)*1000/window); +else + % detect length of raster + max_time=0; + for cell=1:n_cells + data_cell=evalin('base',unit_names{cell}); + max_cell=max(data_cell); + if (max_cell>max_time) + max_time=max_cell; + end + end + max_time=ceil(max_time*1000/window); +end + +% Define raster size +raster = zeros(n_cells,max_time,'logical'); + +% Fill the raster +for cell = 1:n_cells + data_cell = evalin('base',unit_names{cell}); + spikes = round(data_cell*1000/window)+1; + if ~isnan(spikes) + if partial + spikes = spikes(spikes>=initial_time*1000/window); + spikes = spikes(spikes<=final_time*1000/window); + spikes = spikes-floor(initial_time*1000/window)+1; + end + raster(cell,spikes) = true; + end +end \ No newline at end of file diff --git a/common/Get_SEM.m b/common/Get_SEM.m new file mode 100644 index 0000000..ad47e32 --- /dev/null +++ b/common/Get_SEM.m @@ -0,0 +1,8 @@ +function SEM = Get_SEM(data) +% get the standar error of the mean (SEM) +% +% SEM = Get_SEM(data) +% +% By Jesus Perez-Ortega, MArch 2020 + +SEM = std(data)/sqrt(numel(data)); diff --git a/common/Hold_Axes.m b/common/Hold_Axes.m new file mode 100644 index 0000000..5c64fc5 --- /dev/null +++ b/common/Hold_Axes.m @@ -0,0 +1,8 @@ +function h = Hold_Axes(axes_tag) +%% Set the axes specified by its tag +% +% Jesus Perez-Ortega March-18 +h = findobj('Tag',axes_tag); +if(~isempty(h)) + axes(h); +end \ No newline at end of file diff --git a/common/Hold_Figure.m b/common/Hold_Figure.m new file mode 100644 index 0000000..0c497ae --- /dev/null +++ b/common/Hold_Figure.m @@ -0,0 +1,16 @@ +function existFigure = Hold_Figure(titleName) +%% Hold on a Figure with an specific name +% +% existFigure = Hold_Figure(titleName) +% +% Jesus Perez-Ortega March-18 +% Modified Sep 2019 + +h = findobj('name',titleName); + +if isempty(h) + existFigure = false; +else + figure(h); + existFigure = true; +end diff --git a/common/Hotelling_T2_Test.m b/common/Hotelling_T2_Test.m new file mode 100644 index 0000000..aa42994 --- /dev/null +++ b/common/Hotelling_T2_Test.m @@ -0,0 +1,113 @@ +function [p, h1]= Hotelling_T2_Test(X,alpha,mu) +%Hotelling's T-Squared test for one multivariate sample. +% +% Syntax: function [T2Hot1] = T2Hot1(X,alpha) +% +% Inputs: +% X - multivariate data matrix. +% alpha - significance level (default = 0.05). +% +% Output: +% n - sample-size. +% p - variables. +% T2 - Hotelling's T-Squared statistic. +% Chi-sqr. or F - the approximation statistic test. +% df's - degrees' of freedom of the approximation statistic test. +% P - probability that null Ho: is true. +% +% +% Example: For the example given by Johnson and Wichern (1992, p. 183), +% with 20 cases (n = 20) and three variables (p = 3). We are interested +% to test if there is any difference against the expected mean vector +% [4 50 10] at a significance level = 0.10. +% -------------- -------------- +% x1 x2 x3 x1 x2 x3 +% -------------- -------------- +% 3.7 48.5 9.3 3.9 36.9 12.7 +% 5.7 65.1 8.0 4.5 58.8 12.3 +% 3.8 47.2 10.9 3.5 27.8 9.8 +% 3.2 53.2 12.0 4.5 40.2 8.4 +% 3.1 55.5 9.7 1.5 13.5 10.1 +% 4.6 36.1 7.9 8.5 56.4 7.1 +% 2.4 24.8 14.0 4.5 71.6 8.2 +% 7.2 33.1 7.6 6.5 52.8 10.9 +% 6.7 47.4 8.5 4.1 44.1 11.2 +% 5.4 54.1 11.3 5.5 40.9 9.4 +% -------------- -------------- +% +% Total data matrix must be: +% X=[3.7 48.5 9.3;5.7 65.1 8.0;3.8 47.2 10.9;3.2 53.2 12.0;3.1 55.5 9.7; +% 4.6 36.1 7.9;2.4 24.8 14.0;7.2 33.1 7.6;6.7 47.4 8.5;5.4 54.1 11.3; +% 3.9 36.9 12.7;4.5 58.8 12.3;3.5 27.8 9.8;4.5 40.2 8.4;1.5 13.5 10.1; +% 8.5 56.4 7.1;4.5 71.6 8.2;6.5 52.8 10.9;4.1 44.1 11.2;5.5 40.9 9.4]; +% +% Calling on Matlab the function: +% T2Hot1(X) +% Immediately it ask: +% -Do you have an expected mean vector? (y/n): +% For this example we must to put: +% y (meaning 'yes') +% Then it ask: +% -Give me the expected mean vector: +% Giving the mean vector: +% [4 50 10] +% Otherwise (n; meaning 'no') it consider a zero mean vector. +% +% Answer is: +% ------------------------------------------------------------------------------------ +% Sample-size Variables T2 F df1 df2 P +% ------------------------------------------------------------------------------------ +% 20 3 9.7388 2.9045 3 17 0.0649 +% ------------------------------------------------------------------------------------ +% Mean vectors results significant. +% +% +% Created by A. Trujillo-Ortiz and R. Hernandez-Walls +% Facultad de Ciencias Marinas +% Universidad Autonoma de Baja California +% Apdo. Postal 453 +% Ensenada, Baja California +% Mexico. +% atrujo@uabc.mx +% And the special collaboration of the post-graduate students of the 2002:2 +% Multivariate Statistics Course: Karel Castro-Morales, Alejandro Espinoza-Tenorio, +% Andrea Guia-Ramirez. +% +% Copyright (C) December 2002 +% +% References: +% +% Johnson, R. A. and Wichern, D. W. (1992), Applied Multivariate Statistical Analysis. +% 3rd. ed. New-Jersey:Prentice Hall. pp. 180-181,199-200. +% +% Modified by Jesus Perez-Ortega May 2019 + +[n,p]=size(X); + +switch nargin + case 1 + alpha = 0.05; + mu = zeros([1,p]); + case 2 + mu = zeros([1,p]); +end + +if n <= p + error('Requires that sample-size (n) must be greater than the number of variables (p).'); +else + m = mean(X); % Mean vector from data matrix X. + s = cov(X); % Covariance matrix from data matrix X. + T2 = n*(m-mu)*(s\(m-mu)'); % Hotelling's T-Squared statistic. + %T2 = n*(m-mu)*inv(s)*(m-mu)'; % Hotelling's T-Squared statistic. + if n >= 50 + % Chi-square approximation. + p = 1-chi2cdf(T2,p); % Probability that null Ho: is true. + h1 = p < alpha; + else + % F approximation + F = (n-p)/((n-1)*p)*T2; + v2 = n-p; % Denominator degrees of freedom. + p = 1-fcdf(F,p,v2); % Probability that null Ho: is true. + h1 = p < alpha; + end +end diff --git a/common/PlotR_JP.m b/common/PlotR_JP.m new file mode 100755 index 0000000..7b4b475 --- /dev/null +++ b/common/PlotR_JP.m @@ -0,0 +1,52 @@ +% Plot Raster +% +% Plot raster in figure, plot coactive cells and plot cells activity. +% +% [Co Ac] = PlotR_JP(X, figure, titlePlot, color) +% +% Inputs +% X = binary data as matrix FxC (F = #frames, C = #cells) +% numFig = number of the figure to plot. +% titlePlot = title of the plot as string. +% color = color of the plot as vector 1x3 [R G B]. +% +% Outputs +% Figure with plots. +% Co = Coactive cells along the frames. +% Ac = Accumulated activity along the cells. +% +% ..:: by Jesús E. Pérez-Ortega ::.. Mar-2012 + +function [Co Ac] = PlotR_JP(X, numFig, titlePlot, color) + +[F C]=size(X); +figure(numFig); clf + +% Plot Raster +subplot(5,5,[1:4 6:9 11:14 16:19]) +title(titlePlot) +axis([1 F 1 C]) +hold on +for i=1:C + idx=find(X(:,i)); + plot(idx,X(idx,i)*i,'.','color',color) +end +hold off +set(gca,'XTicklabel','') +set(gca,'XTick',0) +box + +% Plot Coactive cells +subplot(5,5,[21:24]) +Co=sum(X,2); +plot(Co,'color',color) +xlim([1 F]) + +% Plot Accumulated activity +subplot(5,5,[5;10;15;20]) +Ac=sum(X,1); +plot(Ac,'*-','color',color) +set(gca,'XTicklabel','') +set(gca,'XTick',0) +xlim([1 C]) +view([90 -90]) \ No newline at end of file diff --git a/common/Plot_Neurons_Measure.m b/common/Plot_Neurons_Measure.m new file mode 100644 index 0000000..9b8229e --- /dev/null +++ b/common/Plot_Neurons_Measure.m @@ -0,0 +1,18 @@ +function Plot_Neurons_Measure(data,name,new_fig) +% Plot measure for neurons in vertical manner +% +% Plot_Neurons_Measure(data,name,new_fig) +% +% By Jesus Perez-Ortega, June 2019 + +if nargin == 2 + new_fig = false; +end + +if new_fig + Set_Figure(name,[0 0 100 400]); +end + +plot(data,'.-k') +view([90 -90]);xlim([1 length(data)]) +ylabel(name) \ No newline at end of file diff --git a/common/Plot_Raster.m b/common/Plot_Raster.m new file mode 100644 index 0000000..9fd6358 --- /dev/null +++ b/common/Plot_Raster.m @@ -0,0 +1,84 @@ +function Plot_Raster(data,name,spikes,reshape,create_figure) +% Plot Raster +% +% Plot_Raster(data,name,spikes,reshape,create_figure) +% +% default: name=''; spikes = true; reshape = true; create_figure = true; +% +% Pérez-Ortega Jesús 2018 +% modified May 2019 +% Modified Sep 2019 + +switch nargin + case 1 + name=''; + spikes = true; + reshape = true; + create_figure = true; + case 2 + spikes = true; + reshape = true; + create_figure = true; + case 3 + reshape = true; + create_figure = true; + case 4 + create_figure = true; +end + +% Get size of the raster +[cells,frames] = size(data); + +% Create a new figure +if create_figure + Set_Figure(['Raster (' name ')'],[0 0 1220 460]); + Set_Axes(['RasterAxes' name],[0 0.34 1 0.66]); +end +axis([0.5 frames 0.5 cells+0.5]); hold on + +% Plot raster +if spikes + if reshape>1 + win = reshape; + elseif reshape + if frames<2000 + win = 1; + elseif frames<4000 + win = 2; + elseif frames<10000 + win = 5; + elseif frames<20000 + win = 10; + elseif frames<500000 + win = 20; + else + win = 50; + end + else + win = 1; + end + + if (nnz(data(:)==0)+nnz(data(:)==1)) ~= numel(data) + imagesc(data,[0 max(data(:))]); + colormap(gca,flipud(gray)) + else + data = Reshape_Raster(data,win); + axis([1 frames/win 1 cells]) + +% hold on +% for i = 1:cells +% id = find(data(i,:)); +% plot(id,data(i,id)*i,'.k') +% end +% set(gca,'XTick',[]) + + imagesc(data,[0,1]); + colormap(gca,[1 1 1; 0 0 0]) + end +else + imagesc(data,[-4 4]); Set_Colormap_Blue_White_Red(); +end +ylim([0.5 cells+0.5]) +box on +set(gca,'XTicklabel','','XTick',[0 frames]) +title(strrep(name,'_','-')) \ No newline at end of file diff --git a/common/Plot_Raster_By_Stimulus.m b/common/Plot_Raster_By_Stimulus.m new file mode 100644 index 0000000..e481cae --- /dev/null +++ b/common/Plot_Raster_By_Stimulus.m @@ -0,0 +1,20 @@ +function Plot_Raster_By_Stimulus(raster,stimulus,before_samples,after_samples,locomotion,fps,save) + +% Get number of stimulus +n_stims = max(stimulus); + +% Plot raster for each kind of stimulus +for i = 1:n_stims + title_name = ['stimulus ' num2str(i)]; + + % Get indices from each stimulus + id = Get_Indices_Before_And_After_Stimulus(stimulus==i,before_samples,after_samples); + + % Plot raster + Plot_Raster_Motion_Stimulus(raster(:,id),title_name,locomotion(id),stimulus(id),fps,max_data) + + % Save Figure + if save + Save_Figure(title_name) + end +end \ No newline at end of file diff --git a/common/Plot_Raster_Motion_Stimulus.m b/common/Plot_Raster_Motion_Stimulus.m new file mode 100644 index 0000000..af45d38 --- /dev/null +++ b/common/Plot_Raster_Motion_Stimulus.m @@ -0,0 +1,58 @@ +function Plot_Raster_Motion_Stimulus(raster,name,locomotion,stimulus,fps,max_data) +% Plot raster, coactivity, moyion and the stimulus +% +% Plot_Raster_Motion_Stimulus(data,name,locomotion,stimulus,fps) +% +% Perez-Ortega Jesus - May 2018 + +if nargin == 5 + max_data = max(raster(:)); +end + +[n_cells,n_frames] = size(raster); + +if n_cells == 1 + Set_Figure(name,[0 0 1200 300]); + subplot(3,1,1) + Plot_Raster(raster,name,true,false,false) + set(gca,'xtick',[],'clim',[0 max_data]) + ylabel('neuron #') + + subplot(3,1,2) + plot(stimulus,'k') + set(gca,'xtick',[],'ytick',[]) + ylabel('stimulus') + xlim([0 n_frames]) + + subplot(3,1,3) + plot(locomotion,'k') + set(gca,'xtick',[],'ytick',[]) + ylabel('motion') + xlim([0 n_frames]) + Set_Label_Time(n_frames,fps) +else + Set_Figure(name,[0 0 1200 700]); + subplot(8,1,1:5) + Plot_Raster(raster,name,true,false,false) + set(gca,'xtick',[],'clim',[0 max_data]) + ylabel('neuron #') + + subplot(8,1,6) + plot(sum(raster),'k') + set(gca,'xtick',[]) + xlim([0 n_frames]) + ylabel('coactivity') + + subplot(8,1,7) + plot(stimulus,'k') + set(gca,'xtick',[],'ytick',[]) + ylabel('stimulus') + xlim([0 n_frames]) + + subplot(8,1,8) + plot(locomotion,'k') + set(gca,'xtick',[],'ytick',[]) + ylabel('motion') + xlim([0 n_frames]) + Set_Label_Time(n_frames,fps) +end \ No newline at end of file diff --git a/common/Read_Colors.m b/common/Read_Colors.m new file mode 100644 index 0000000..6aa0bf3 --- /dev/null +++ b/common/Read_Colors.m @@ -0,0 +1,77 @@ +function colors = Read_Colors(num_colors) +% Return the numbers of colors specified +% +% colors = Read_Colors(num_colors) +% +% colors=[0.9 0.0 0.0;... % red +% 0.5 0.8 0.0;... % green +% 0.0 0.5 1.0;... % sky blue +% 0.9 0.9 0.0;... % yellow +% 0.9 0.0 0.9;... % magenta +% 0.0 0.9 0.9;... % cyan +% 1.0 0.5 0.0;... % orange +% 0.5 0.0 1.0;... % purple +% 1.0 0.7 0.7;... % rose +% 0.6 0.3 0.0;... % brown +% 0.95 0.5 0.5;... % red light +% 0.75 0.9 0.5;... % green light +% 0.5 0.75 1.0;... % sky blue light +% 0.95 0.95 0.5;... % yellow light +% 0.95 0.5 0.95;... % magenta light +% 0.5 0.95 0.95;... % cyan light +% 1.0 0.75 0.5;... % orange light +% 0.75 0.5 1.0;... % purple light +% 1.0 0.85 0.85;... % rose light +% 0.8 0.65 0.5;... % brown light +% 0.4 0.0 0.0;... % red dark +% 0.0 0.3 0.0;... % green dark +% 0.0 0.0 0.5;... % sky blue dark +% 0.4 0.4 0.0;... % yellow dark +% 0.4 0.0 0.4;... % magenta dark +% 0.0 0.4 0.4;... % cyan dark +% 0.5 0.0 0.0;... % orange dark +% 0.0 0.0 0.5;... % purple dark +% 0.5 0.2 0.2;... % rose dark +% 0.2 0.1 0.0]; % brown dark +% +% Perez-Ortega Jesus - May 2018 + +if nargin==0 + num_colors=10; +end + +if(num_colors<=10) + colors=[0.9 0.0 0.0;... % red + 0.5 0.8 0.0;... % green + 0.0 0.5 1.0;... % sky blue + 0.9 0.9 0.0;... % yellow + 0.9 0.0 0.9;... % magenta + 0.0 0.9 0.9;... % cyan + 1.0 0.5 0.0;... % orange + 0.5 0.0 1.0;... % purple + 1.0 0.7 0.7;... % rose + 0.6 0.3 0.0;... % brown + 0.95 0.5 0.5;... % red light + 0.75 0.9 0.5;... % green light + 0.5 0.75 1.0;... % sky blue light + 0.95 0.95 0.5;... % yellow light + 0.95 0.5 0.95;... % magenta light + 0.5 0.95 0.95;... % cyan light + 1.0 0.75 0.5;... % orange light + 0.75 0.5 1.0;... % purple light + 1.0 0.85 0.85;... % rose light + 0.8 0.65 0.5;... % brown light + 0.4 0.0 0.0;... % red dark + 0.0 0.3 0.0;... % green dark + 0.0 0.0 0.5;... % sky blue dark + 0.4 0.4 0.0;... % yellow dark + 0.4 0.0 0.4;... % magenta dark + 0.0 0.4 0.4;... % cyan dark + 0.5 0.0 0.0;... % orange dark + 0.0 0.0 0.5;... % purple dark + 0.5 0.2 0.2;... % rose dark + 0.2 0.1 0.0]; % brown dark + colors=colors(1:num_colors,:); +else + colors=hsv(num_colors); +end \ No newline at end of file diff --git a/common/Reshape_Raster.m b/common/Reshape_Raster.m new file mode 100644 index 0000000..9674454 --- /dev/null +++ b/common/Reshape_Raster.m @@ -0,0 +1,21 @@ +% Reshape raster +% +% Binarize the raster by a given window +% +% Pérez-Ortega Jesús - May 2018 + +function new_raster=Reshape_Raster(raster,window) + [c,n]=size(raster); + + if(window==1) + new_raster = raster; + else + new_n=floor(n/window); + new_raster=zeros(c,new_n); + for i=1:new_n + ini=(i-1)*window+1; + fin=i*window; + new_raster(:,i)=logical(sum(raster(:,ini:fin),2)); + end + end +end \ No newline at end of file diff --git a/common/Save_Figure.m b/common/Save_Figure.m new file mode 100644 index 0000000..fd5e814 --- /dev/null +++ b/common/Save_Figure.m @@ -0,0 +1,77 @@ +function Save_Figure(name,format,path,background,lineWidth) +% Save current figure with default style +% +% Save_Figure(name,format,path,background,lineWidth) +% +% default: format = ''; path = ''; background = ''; lineWidth = ''; +% +% Perez-Ortega Jesus - June 2018 +% modified Feb 2019 +% modified May 2019 +% modified Sep 2019 +% modified May 2020 + +switch nargin + case 1 + format = ''; + path = ''; + background = ''; + lineWidth = ''; + case 2 + s.Format = format; + path = ''; + background = ''; + lineWidth = ''; + case 3 + s.Format = format; + background = ''; + lineWidth = ''; + case 4 + s.Format = format; + s.Background = background; + lineWidth = ''; + case 5 + s.Format = format; + s.Background = background; + s.Background = background; + s.FixedLineWidth = lineWidth; +end + +if isempty(format), s.Format = 'png'; end +if isempty(background), s.Background = 'w'; end +if isempty(lineWidth), s.FixedLineWidth='1'; end + +s.Version='1'; +s.Preview='none'; +s.Width='auto'; +s.Height='auto'; +s.Units='inches'; +s.Color='rgb'; +s.FixedFontSize='14'; +s.ScaledFontSize='auto'; +s.FontMode='fixed'; +s.FontSizeMin='8'; +s.ScaledLineWidth='auto'; +s.LineMode='fixed'; % 'none' or 'fixed' +s.LineWidthMin='0.5'; +s.FontName='auto'; +s.FontWeight='auto'; +s.FontAngle='auto'; +s.FontEncoding='latin1'; +s.PSLevel='3'; +s.Renderer='auto'; +s.Resolution='300'; +s.LineStyleMap='none'; +s.ApplyStyle='0'; +s.Bounds='loose'; +s.LockAxes='on'; +s.LockAxesTicks='off'; +s.ShowUI='on'; +s.SeparateText='off'; + +if strcmp(s.Format,'eps') + set(gcf,'renderer','Painters') + print('-depsc','-tiff','-r300', '-painters', [path name '.eps']) +else + hgexport(gcf,[path name '.' s.Format],s); +end \ No newline at end of file diff --git a/common/Scale.m b/common/Scale.m new file mode 100644 index 0000000..3cfab2b --- /dev/null +++ b/common/Scale.m @@ -0,0 +1,36 @@ +function data_out = Scale(data,min_value,max_value) +% Scale data between 0 and 1 +% +% data_out = Scale(data,min_value,max_value) +% +% Jesus Perez-Ortega April-19 + +if nargin==1 + min_value = 0; + max_value = 1; +end + +% Change to double +data_out = double(data); + +% Substract minimum +min_data = min(data_out(:)); +data_out = data_out-min_data; + +% Divide by maximum +max_data = max(data_out(:)); +if max_data + data_out = data_out./max_data; + + % Multiply by factor + factor = max_value-min_value; + data_out = data_out.*factor; + + % Add minimum + data_out = data_out+min_value; +else + data_out = min_value*ones(size(data)); +end + + + diff --git a/common/Set_Axes.m b/common/Set_Axes.m new file mode 100644 index 0000000..6875a59 --- /dev/null +++ b/common/Set_Axes.m @@ -0,0 +1,11 @@ +function h = Set_Axes(AxesName,Position) +% Set Axes +% +% h = Set_Axes(AxesName,Position) + +ax = axes('outerposition',Position); +set(ax,'Tag',AxesName) + +if nargout + h = ax; +end \ No newline at end of file diff --git a/common/Set_Figure.m b/common/Set_Figure.m new file mode 100644 index 0000000..e94dc29 --- /dev/null +++ b/common/Set_Figure.m @@ -0,0 +1,36 @@ +function fig = Set_Figure(title_name,position) +%% Set Figure +% +% Create a figure with basic properties: title, size and position. If the +% figure already exist, it is only cleaned. +% +% fig = Set_Figure(title_name,position) +% +% position = [x y width height]) +% +% Jesus Perez-Ortega modififed April-19 +% Modified Oct 2019 + +if isempty(title_name) + if nargin == 1 + h = figure('name','unnamed','numbertitle','off'); + else + h = figure('name','unnamed','numbertitle','off','position',position); + end +else + h = findobj('name',title_name); + if isempty(h) + if nargin == 1 + h = figure('name',title_name,'numbertitle','off'); + else + h = figure('name',title_name,'numbertitle','off','position',position); + end + else + figure(h); + end +end +clf + +if nargout + fig = h; +end diff --git a/common/Set_Label_Time.m b/common/Set_Label_Time.m new file mode 100644 index 0000000..40f56b6 --- /dev/null +++ b/common/Set_Label_Time.m @@ -0,0 +1,58 @@ +function Set_Label_Time(samples,sampleFrequency,sampleShift) +% Set convenient label of time +% +% Set_Label_Time(samples,sampleFrequency) +% +% Jesus Perez-Ortega April-19 +% Modified May 2019 +% Modified Oct 2019 +% Modified Nov 2019 +% Modified Jul 2020 + +if nargin==2 + sampleShift = 0; +end + +% Plot depending on the duration +if(samples/sampleFrequency<=1) % Less than 1 s + step = 0.1; + xlabel('Time (ms)'); + factor = 0.001; +elseif(samples/sampleFrequency<=3) % Less than 2 s + step = 0.5; + xlabel('Time (ms)'); + factor = 0.001; +elseif(samples/sampleFrequency<=15) % Less than 10 s + step = 1; + xlabel('Time (s)'); + factor = 1; +elseif(samples/sampleFrequency<=30) % Less than 30 s + step = 5; + xlabel('Time (s)'); + factor = 1; +elseif(samples/sampleFrequency/60<3) % Less than 3 min + step = 20; + xlabel('Time (s)'); + factor = 1; +elseif(samples/sampleFrequency/60<30) % Less than 30 min + step = 60; + xlabel('Time (min)'); + factor = 60; +elseif(samples/sampleFrequency/60<60) % Less than 60 min + step = 60*5; + xlabel('Time (min)'); + factor = 60; +elseif(samples/sampleFrequency/60<180) % Less than 180 min + step = 60*30; + xlabel('Time (min)'); + factor = 60; +else + step = 60*60; + xlabel('Time (h)'); + factor = 3600; +end + +set(gca,'box','off','xtick',(1:step*sampleFrequency:(samples+1))+sampleShift,... + 'xticklabel',(0:samples/sampleFrequency/step)*step/factor) +xlim([1 samples+1]) + diff --git a/common/Shuffle_Raster_Test.m b/common/Shuffle_Raster_Test.m new file mode 100644 index 0000000..441b912 --- /dev/null +++ b/common/Shuffle_Raster_Test.m @@ -0,0 +1,90 @@ +function [threshold,cdf] = Shuffle_Raster_Test(raster,shuffle_method,iterations,alpha,plot_figure) +% Test for significant number of coactive neurons +% +% function [threshold,cpd] = Shuffle_Raster_Test(raster,shuffle_method,iterations,alpha, +% plot_figure) +% +% default: shuffled_method = 'time_shift'; iterations = 1000; alpha = 0.05; +% plot_figure = false +% +% Perez-Ortega Jesus - June 2019 + +switch nargin + case 1 + shuffle_method='time_shift'; + iterations = 1000; + alpha=0.05; + plot_figure = false; + case 2 + iterations = 1000; + alpha=0.05; + plot_figure = false; + case 3 + alpha=0.05; + plot_figure = false; + case 4 + plot_figure = false; +end + +% Get raster +coactivity = sum(raster); + +% Initialize variables +samples = size(raster,2); +coactivity_shuffled = zeros(iterations,samples); + +% Plot only if necessary +if plot_figure + Set_Figure('Raster | real VS shuffled',[0 0 1000 800]); + subplot(5,1,1) + imagesc(raster) + title('Observed') + colormap([1 1 1; 0 0 0]) + Set_Figure('Coactivity | real VS shuffled',[0 0 1000 800]); + subplot(5,1,1) + plot(coactivity) + title('Observed') + y_limits=get(gca,'ylim'); +end + +% SHUFFLED data +for i = 1:iterations + % Get raster shuffled + raster_shuffled = shuffle(raster,shuffle_method); + + % Get coactivity + coactivity_shuffled(i,:) = sum(raster_shuffled); + + % Plot only if necessary + if(plot_figure && i<5) + Hold_Figure('Raster | real VS shuffled'); + subplot(5,1,i+1) + imagesc(raster_shuffled) + title(['Shuffled ' num2str(i)]) + Hold_Figure('Coactivity | real VS shuffled') + subplot(5,1,i+1) + plot(coactivity_shuffled(i,:)) + title(['Shuffled ' num2str(i)]) + ylim(y_limits) + end +end + +% cumulative distribution function of coactive cells +interval = 0:max(coactivity_shuffled(:)); +cdf = histcounts(coactivity_shuffled(:),interval,'Normalization','cdf'); + +% Set threshold from Probability distribution < alpha +threshold = find(cdf>(1-alpha),1,'first'); + +% Plot distribution shuffled +if plot_figure + Set_Figure('Shuffled test',[0 0 600 300]); + plot(interval(1:end-1),cdf,'o-b'); hold on + title({['Cumulative distribution of raster shuffled (' num2str(iterations) ' iterations; method '... + strrep(shuffle_method,'_','-') ')'];['th = ' num2str(threshold) ' at alpha=' num2str(alpha)]}) + ylabel('P(x)') + xlabel('x') + if(~isnan(threshold)) + xlim([0 threshold*3]) + end +end \ No newline at end of file diff --git a/common/Validate_Name.m b/common/Validate_Name.m new file mode 100644 index 0000000..e72cf12 --- /dev/null +++ b/common/Validate_Name.m @@ -0,0 +1,21 @@ +function validated = Validate_Name(name) +% Validate string to write a variable name in worskpace +% +% validated = Validate_Name(name) +% +% Jesus Perez-Ortega May-19 + +no_alphanum = ~isstrprop(name,'alphanum'); + +ini = find(no_alphanum==0,1,'first'); +fin = find(no_alphanum==0,1,'last'); + +validated = name(ini:fin); +no_alphanum = no_alphanum(ini:fin); + +validated(no_alphanum) = '_'; +n = sum(no_alphanum); + +for i = 1:n + validated = strrep(validated,'__','_'); +end \ No newline at end of file