diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000..e01012e Binary files /dev/null and b/.DS_Store differ diff --git a/CCidx2_JP.m b/CCidx2_JP.m new file mode 100644 index 0000000..62dc96c --- /dev/null +++ b/CCidx2_JP.m @@ -0,0 +1,39 @@ +% Cluster Connectivity index 2 +% Get cluster connectivity index of neuronal groups of peaks. +% Mean of the groups' connectivity index +% (This is an idea of Jesús E. Pérez-Ortega, José Bargas, et al.) +% +% [CCI] = CCidx2_JP(g,Xp,idx) +% +% Inputs +% g = number of groups +% Xp = binary data as matrix PxC (P = #peaks, C = #cells) +% idx = indexes of group to which each data point belongs +% +% Outputs +% CCI = Connectivity index +% +% ..:: by Jesús E. Pérez-Ortega ::.. Mar-2012 + +function [CCI] = CCidx2_JP(g,Xp,idx) + +P=size(Xp,1); +C_intra=zeros(g,1); % connectivity inter-group +single=0; % single neuronal vector +for i=1:g + clust = find(idx==i); + PM=Xp(clust,:); + idxPM=find(sum(PM,1)); + PM=PM(:,idxPM); + if size(PM,1)==1 % if single not be taken into account + C_intra(i,1)=0; % unlike with the first version: 1 instead of 0 + single=single+1; + else + C_intra(i,1)=CI_JP(PM); + end +end +% disp(C_intra) +CCI=sum(C_intra)/g; % unlike with the first version: g instead of (g-single) + + + diff --git a/Change_Names_MEA_Channels.m b/Change_Names_MEA_Channels.m new file mode 100644 index 0000000..90e15e0 --- /dev/null +++ b/Change_Names_MEA_Channels.m @@ -0,0 +1,15 @@ +% Exclusivo para MEAs de 128 canales + +function new_channels = Change_Names_MEA_Channels(channels) + + ch=channels(:,1); + ends=[find(diff(ch)~=0); length(ch)]; + inis=[1; ends(1:end-1)+1]; + + for i=1:length(inis) + new(inis(i):ends(i))=i; + end + + channels(:,1)=new; + new_channels=channels; +end \ No newline at end of file diff --git a/Circle_Mask.m b/Circle_Mask.m new file mode 100644 index 0000000..129e0e0 --- /dev/null +++ b/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/ClustIdx_JP.m b/ClustIdx_JP.m new file mode 100644 index 0000000..5131f4c --- /dev/null +++ b/ClustIdx_JP.m @@ -0,0 +1,62 @@ +% Clustering indexes +% Get indexes for evaluating clustering from hierarchical cluster tree +% +% [ClustIdx g] = ClustIdx_JP(Tree, SimOrXpeaks, Metric, numFig) +% +% Inputs +% Tree = hierarchical cluster tree +% SimOrXpeaks = Sim, similarity as matrix PxP (P=#peaks) for metrics Dunn & +% Contrast; Xpeaks, peaks vectors as matrix PxC for metrics +% Connectivity & Davies +% (P = #peaks; C=#cells) +% Metric = index to compute ('Dunn','Connectivity','Davies','Contrast') +% numFig = number of the figure to plot +% +% Output +% ClustIdx = clustering indexes of 'Metric' from 2 to 10 groups +% g = best number of groups according to the index selected +% +% ..:: by Jesús E. Pérez-Ortega ::.. Jun-2012 +% modified March-2018 + +function [ClustIdx g] = ClustIdx_JP(tree_or_data,similitud,method,num_fig, clust_method,groups) + +if(nargin==5) + groups=2:30; +end + +dist=1-similitud; +j=1; +for i=groups + switch(clust_method) + case 'hierarchical' + T = cluster(tree_or_data,'maxclust',i); + case 'kmeans' + T = kmeans(tree_or_data,i); + end + g=max(T); + + switch(method) + case 'Dunn' + ClustIdx(j)=DunnIdx_JP(g,dist,T); + case 'Connectivity' + ClustIdx(j)=CCidx2_JP(g,similitud,T); + case 'Davies' + ClustIdx(j)=DaviesIdx_JP(g,similitud,T); + case 'Contrast' + ClustIdx(j)=ContrastIdx_JP(g,similitud,T); + otherwise + method='Dunn'; + ClustIdx(j)=DunnIdx_JP(g,dist,T); + end + j=j+1; +end +figure(num_fig) +plot(groups,ClustIdx) +hold on + +[gIdx,g]=max(ClustIdx); + +plot(groups(g),gIdx,'*r') +hold off +title([method '''s index (' num2str(groups(g)) ' groups recommended)']) diff --git a/Compare_Times.m b/Compare_Times.m new file mode 100644 index 0000000..c2b3e9b --- /dev/null +++ b/Compare_Times.m @@ -0,0 +1,89 @@ +% Times comparison Manual VS Auto + +function [coincidence,diffs] = Compare_Times(times_1,times_2,coactivity) + if(nargin==2) + coactivity=[]; + end + + n_1=size(times_1,1); + n_2=size(times_2,1); + Times_1=times_1(:); + Times_2=times_2(:); + + Y_manual=1.5*ones(length(Times_1),1); + Y_auto=1.5*ones(length(Times_2),1); + + % Comparison + n_coincidence=0; + auto_extra=0; + manual_extra=0; + auto_actual=1; + coincidence=[]; + for i=1:n_1 + i_m=times_1(i,1); + f_m=times_1(i,2); + + for j=auto_actual:n_2 + i_a=times_2(j,1); + f_a=times_2(j,2); + + % If data is in the same place + %if(i_m>=i_a && i_m<=f_a || f_m>=i_a && f_m<=f_a || i_m<=i_a && f_m>=f_a) + if abs(i_m-i_a)<100 + diff_i=i_a-i_m; + diff_f=f_a-f_m; + n_coincidence=n_coincidence+1; + coincidence(n_coincidence,:)=[i_m i_a f_m f_a]; + diffs(n_coincidence,:)=[diff_i diff_f]; + auto_actual=j+1; + break; + % If are differents + elseif(i_m>f_a) + auto_extra=auto_extra+1; + elseif(f_m 0); + case 2 + id = find(raster_analysis.Peaks.Indices(range) > 0)+range(1)-1; + case 3 + if(peaks) + id = find(raster_analysis.Peaks.Indices(range) > 0)+range(1)-1; + else + id = find(raster_analysis.Peaks.Indices(range) == 0)+range(1)-1; + end + end + raster_peaks = raster_analysis.Raster.Raster(:,id); +end \ No newline at end of file diff --git a/Filter_MEA_Channel.m b/Filter_MEA_Channel.m new file mode 100644 index 0000000..9dcb3d0 --- /dev/null +++ b/Filter_MEA_Channel.m @@ -0,0 +1,28 @@ +% Filter MEA +% +% Lab Peńa-Ortega +% 250 - 7,500 Hz band pass filter +% butterworth order 2 +% +% By Pérez-Ortega Jesús, jan-2019 + +function filtered = Filter_MEA_Channel(signal) + + fs = 25000; % sampling frequency + + % Filter design + order=2; + % low pass + fc = 7500; % cutoff frequency + [b,a] = butter(order,fc/(fs/2),'low'); + + % Filter 1 + y = filter(b,a,signal); + + % high pass + fc = 250; % cutoff frequency + [b,a] = butter(order,fc/(fs/2),'high'); + + % Filter 2 + filtered = filter(b,a,y); +end \ No newline at end of file diff --git a/Find_Peaks_Or_Valleys.m b/Find_Peaks_Or_Valleys.m new file mode 100644 index 0000000..474d8f1 --- /dev/null +++ b/Find_Peaks_Or_Valleys.m @@ -0,0 +1,269 @@ +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) +% +% 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 +% JP debug 1-oct-2012 +% JP add 3rd input to join or not +% modified March-2018 +% modified April-2018 +% modified May-2018 +% modified Jul-2018 +% modified Nov-2018 +% modified Dec-2018 +% modified Jan-2019 + +if(nargin==6) + ignore_ini_fin=false; +elseif(nargin==5) + fixed_width = 0; + ignore_ini_fin=false; +elseif(nargin==4) + minimum_width=0; + fixed_width=0; + ignore_ini_fin=false; +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 = flipud(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<0) + ini=0; + 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 + data(fixed_width_peak-fixed_width)=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/Find_Threshold_In_Cumulative_Distribution.m b/Find_Threshold_In_Cumulative_Distribution.m new file mode 100644 index 0000000..3c6a6d4 --- /dev/null +++ b/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/Fit_Power_Law.m b/Fit_Power_Law.m new file mode 100644 index 0000000..b211554 --- /dev/null +++ b/Fit_Power_Law.m @@ -0,0 +1,18 @@ +% Fit_Power_Law +% +% Jesús E. Pérez-Ortega - june 2018 + +function [slope, intercept, R2] = Fit_Power_Law(x,y) + + % Set power law function + power_law = fittype('a*x^b'); + + % Compute fit + [fit1,gof] = fit(x,y,power_law,'StartPoint',[1 -1]); + coeffs = coeffvalues(fit1); + + % Get coefficients + slope=coeffs(2); + intercept=coeffs(1); + R2=gof.rsquare; +end diff --git a/Get_Adjacency_From_Raster.m b/Get_Adjacency_From_Raster.m new file mode 100644 index 0000000..64eaef0 --- /dev/null +++ b/Get_Adjacency_From_Raster.m @@ -0,0 +1,27 @@ +% Get adjacency matrix from raster peaks +% +% Negative numbers and NaNs are set to 0s +% +% Pérez-Ortega Jesús - march 2018 +% Modified - may 2018 +% Modified - june 2018 + +function adjacency=Get_Adjacency_From_Raster(raster,connectivity_method) + cells=size(raster,1); + switch(connectivity_method) + case 'coactivity' + %warning('Data is set to binary to obtain the adjacency matrix.') + raster=double(raster>0); + adjacency=raster*raster'.*(1-eye(cells)); + case 'jaccard' + %warning('Data is set to binary to obtain the adjacency matrix.') + raster=double(raster>0); + adjacency=squareform(1-pdist(raster,'jaccard')); + adjacency(isnan(adjacency))=0; + otherwise + adjacency=corr(raster','type',connectivity_method); + adjacency(adjacency<0)=0; + adjacency(isnan(adjacency))=0; + adjacency=adjacency.*(ones(cells)-eye(cells)); + end +end \ No newline at end of file diff --git a/Get_Adjacency_From_Sequence.m b/Get_Adjacency_From_Sequence.m new file mode 100644 index 0000000..960dde2 --- /dev/null +++ b/Get_Adjacency_From_Sequence.m @@ -0,0 +1,18 @@ +% Get a directed adjacency matrix from a sequence +% +% Get a directed adjacency matrix including loops. +% +% Pérez-Ortega Jesús - march 2018 + +function Adjacency = Get_Adjacency_From_Sequence(sequence) + n=length(sequence); + max_state=max(sequence); + Adjacency=zeros(max_state); + + for i=1:n-1 + j=i+1; + state_i=sequence(i); + state_j=sequence(j); + Adjacency(state_i,state_j)=Adjacency(state_i,state_j)+1; + end +end \ No newline at end of file diff --git a/Get_And_Filter_Coactivity.m b/Get_And_Filter_Coactivity.m new file mode 100644 index 0000000..dd00026 --- /dev/null +++ b/Get_And_Filter_Coactivity.m @@ -0,0 +1,18 @@ +% Get Coactivity +% +% double smoothing filter +% +% By Jesús Pérez-Ortega jan-2018 +% modified feb-2018 +% modified april-2018 + +function smooth_coactivity = Get_And_Filter_Coactivity(raster,bin) + coactivity=sum(raster)'; + if(bin>1) + % double smoothing + smooth_coactivity=smooth(coactivity,bin); + %smooth_coactivity=smooth(smooth_coactivity,bin); + else + smooth_coactivity=coactivity; + end +end \ No newline at end of file diff --git a/Get_Circular_XY.m b/Get_Circular_XY.m new file mode 100644 index 0000000..882a424 --- /dev/null +++ b/Get_Circular_XY.m @@ -0,0 +1,27 @@ +function [xy,step] = Get_Circular_XY(n,radio,offset) +% Get circular coordinates for n elements given a radio and offset in +% degrees +% +% [xy,step] = Get_Circular_XY(n,radio,offset) +% +% Pérez-Ortega Jesús - April 2018 +% modified April 2019 + +switch(nargin) + case 1 + radio = 1; + offset = 0; + case 2 + offset = 0; +end + +% convert to radians +offset = offset*pi/180; + +if (n==1) + xy=[0 0]; + step = 0; +else + xy=[cos(offset+2*pi*(1:n)'/n) sin(offset+2*pi*(1:n)'/n)].*radio; + step = 360/n; +end \ No newline at end of file diff --git a/Get_Efficiency.m b/Get_Efficiency.m new file mode 100644 index 0000000..e290694 --- /dev/null +++ b/Get_Efficiency.m @@ -0,0 +1,24 @@ +% Get efficiency between a group of vectors +% +% Jesus Perez-Ortega - Sep 2018 + +function [Efficiency,Error] = Get_Efficiency(vectors) + [n,length_vector]=size(vectors); + if(n>1) + % Compute efficiency + Error=zeros(n); + for i=1:n-1 + for j=i+1:n + error_ij=length(find(abs(vectors(i,:)-vectors(j,:)))); + Error(i,j)=error_ij; + Error(j,i)=error_ij; + end + end + Error=Error./length_vector; + Efficiency=1-Error; + else + warning('There are no data to compare!') + Efficiency=[]; + Error=[]; + end +end diff --git a/Get_Entropy.m b/Get_Entropy.m new file mode 100644 index 0000000..f22b169 --- /dev/null +++ b/Get_Entropy.m @@ -0,0 +1,22 @@ +function entropy = Get_Entropy(frequencies) +% Entropy from a list of frequencies +% +% By Pérez-Ortega Jesús, Feb 2019 + +n = length(frequencies); +total = sum(frequencies); + +if(~total) + entropy = nan; + return +end + +p = frequencies/total; + +entropy = 0; +for i = 1:n + if(p(i)) + entropy = entropy + p(i)*log2(p(i)); + end +end +entropy = -entropy; diff --git a/Get_Force_XY.m b/Get_Force_XY.m new file mode 100644 index 0000000..4b73d3a --- /dev/null +++ b/Get_Force_XY.m @@ -0,0 +1,21 @@ +function xy = Get_Force_XY(adjacency) +% Get coordinates to plot network in 2D with force layout +% +% xy = Get_Force_XY(adjacency) +% +% Jesús Pérez-Ortega nov 2018 + +if isempty(adjacency) + xy = []; + return +end + +% Get coordinates +G=graph(adjacency); +figure(); +%G_plot=plot(G,'Layout','force','usegravity',true); +%G_plot=plot(G,'Layout','force','weighteffect','direct'); +%G_plot=plot(G,'Layout','force','usegravity',true,'iterations',250); +G_plot=plot(G,'Layout','force','iterations',250); +xy = [G_plot.XData' G_plot.YData']; +close; \ No newline at end of file diff --git a/Get_Groups_Of_Neurons.m b/Get_Groups_Of_Neurons.m new file mode 100644 index 0000000..93b2243 --- /dev/null +++ b/Get_Groups_Of_Neurons.m @@ -0,0 +1,51 @@ +% Get groups of neurons +% +% By Jesús Pérez-Ortega jan-2018 + +function [CellsGroups Indexes IndexesComb] = Get_Groups_Of_Neurons(Freqs,Groups,idx_vector_state,th_freq_group) + + % Get states of each neuron + C=size(Freqs,1); + CellsGroups=zeros(Groups,C); + CellsAcGroups=zeros(Groups,C); + for i=1:Groups + PeaksGroupIdx= idx_vector_state==i; + n_peaks_group=length(PeaksGroupIdx); + CellsAcGroups(i,:)=sum(Freqs(:,PeaksGroupIdx),2)/n_peaks_group; + CellsGroup=CellsAcGroups(i,:)>th_freq_group; + CellsGroups(i,CellsGroup)=i; + end + + % Get states combinations + [Groups C]=size(CellsGroups); + States=[]; Labels={}; + Idx=[]; Indexes=[]; + Combinations=2^Groups-1; + for i=1:Combinations + bin_number=dec2bin(i,Groups); + GroupsToCompare=[]; + Label=[]; + for j=1:Groups + Compare=str2num(bin_number(j)); + if Compare + GroupsToCompare=[GroupsToCompare; CellsGroups(j,:)>=1]; + if Label + Label=[Label ',' num2str(j)]; + else + Label=num2str(j); + end + else + GroupsToCompare=[GroupsToCompare; CellsGroups(j,:)==0]; + end + end + Idx{i}=find(sum(GroupsToCompare,1)==Groups); + States(i)=length(Idx{i})/C; + Labels(i)={Label}; + end + [Labels idxSort]=sort(Labels); + States=States(idxSort); + for i=1:Combinations + Indexes=[Indexes Idx{idxSort(i)}]; + IndexesComb(Idx{i})=i; + end +end \ No newline at end of file diff --git a/Get_ISI.m b/Get_ISI.m new file mode 100644 index 0000000..e29a1b1 --- /dev/null +++ b/Get_ISI.m @@ -0,0 +1,11 @@ +% Get Inter spike interval (ISI) +% +% inputs: data = vector with spikes (0: no spike, 1: spike) +% f = samples per second (sampling frequency) +% +% By Jesús E. Pérez-Ortega - jan - 2019 + +function isi = Get_ISI(data,f) + id = find(data); + isi = diff(id)/f; +end \ No newline at end of file diff --git a/Get_Instant_Freq.m b/Get_Instant_Freq.m new file mode 100644 index 0000000..7d5f647 --- /dev/null +++ b/Get_Instant_Freq.m @@ -0,0 +1,34 @@ +% Plot instantaneous frequency (inverse of ISI) +% +% inputs: data = vector with spikes (0: no spike, 1: spike) +% sps = samples per second (sampling frequency) +% +% By Jesús E. Pérez-Ortega - may 2017 +% modified jan-2018 + +function frequencies = Get_Instant_Freq(data,sps) + + N=length(data); + frequencies=zeros(1,N); + + idx=find(data); + + % if there are spikes + if (~isempty(idx)) + % If first data is not a spike + if (idx(1)~=1) + frequencies(1:idx(1))=sps/(idx(1)-1); + end + + % All ISI inverse + n=length(idx)-1; + for j=1:n + frequencies(idx(j):idx(j+1))=sps/(idx(j+1)-idx(j)); + end + + % If last data is not a spike + if(idx(n+1)~=N) + frequencies(idx(n+1):N)=sps/(N-idx(n+1)); + end + end +end \ No newline at end of file diff --git a/Get_Network_Activation_Sequence.m b/Get_Network_Activation_Sequence.m new file mode 100644 index 0000000..604b6e1 --- /dev/null +++ b/Get_Network_Activation_Sequence.m @@ -0,0 +1,74 @@ +% Get network activation sequence +% +% Jesús Pérez-Ortega - Dic 2018 +% modified Jan 2019 +function [sequence,singles] = Get_Network_Activation_Sequence(raster,network,link_sequence) + + if(nargin==2) + link_sequence = false; + end + + %[shorted,short_times]=Delete_Consecutive_Coactivation(raster); + n_length=size(raster,2); + adjacency=Get_Adjacency_From_Raster(raster,'coactivity'); + coincidence=double(adjacency>0).*network; + + if(link_sequence) + coincidence = squareform(coincidence,'tovector'); + end + + %total_links=sum(adjacency_core_inspiration(:)); + %n_co=1; + %single_links=0; + + sequence=[]; + for i=1:n_length + single=Get_Adjacency_From_Raster(raster(:,i),'coactivity'); + if(link_sequence) + single=squareform(single,'tovector'); + end + single_coincidence=single.*coincidence; + sequence=[sequence setdiff(find(sum(single_coincidence)),sequence)]; + singles{i}=single_coincidence; + + % Plot network + %{ + core_coactivation=single.*adjacency_core_inspiration; + single_links=single_links+sum(coactivation_remaining(:)); + + if(sum(coactivation_remaining(:))~=0) + %Set_Figure('Network',[0 0 300 300]); + net_color=colors(short_times(l),:); + node_size=10; + subplot(5,11,n_co+11*(m-1)) + %Plot_WU_Network(core_coactivation,[],net_color,mean([net_color; 0.5 0.5 0.5]),node_size); hold on + Plot_WU_Network(coactivation_remaining,[],net_color,mean([net_color; 0.5 0.5 0.5]),node_size) + title([num2str(short_times(l)) ' ms (' num2str(round(single_links/total_links*100)) '%)']) + if(n_co==1) + ylabel(['inspiration #' num2str(k)]) + end + subplot(5,11,11*m) + Plot_WU_Network(coactivation_remaining,[],net_color,mean([net_color; 0.5 0.5 0.5]),node_size);hold on + % frame = getframe(gcf); + % writeVideo(v,frame); + n_co=n_co+1; + if(n_co>10) + break; + end + end + remaining=remaining-coactivation_remaining; + if(remaining==0) + break; + end + %input('press any key to continue...') + %} + end + %{ + m=m+1; + if(m>5) + Save_Figure(['Networks_' num2str(j) '_' num2str(ceil(k/5))]); + close + m=1; + end + %} +end \ No newline at end of file diff --git a/Get_Normalized_Instant_Freq.m b/Get_Normalized_Instant_Freq.m new file mode 100644 index 0000000..f399602 --- /dev/null +++ b/Get_Normalized_Instant_Freq.m @@ -0,0 +1,39 @@ +% Get Normalized Instantaneous Frequency from raster +% +% Instantaneous frequency = inverse of ISI +% +% By Jesús Pérez-Ortega jan-2018 +% modified feb-2018 +% modified jan-2019 + +function freqs = Get_Normalized_Instant_Freq(data,fps,norm_type,bin) + + if(nargin==3) + bin=0; + end + + [c,f]=size(data); + freqs=zeros([c f]); + for i=1:c + % Get frequencies + freqs(i,:)=Get_Instant_Freq(data(i,:),fps); + + if(~freqs(i,:)) + freqs(i,:)=zeros(1,f); + else + % Normalize the firing frequency + if (bin) + freqs(i,:)=smooth(freqs(i,:),bin); + end + if (~strcmp(norm_type,'none')) + freqs(i,:)=freqs(i,:)-mean(freqs(i,:)); + switch(norm_type) + case 'norm' + freqs(i,:)=freqs(i,:)/max(freqs(i,:)); % Normalize + case 'zscore' + freqs(i,:)=freqs(i,:)/std(freqs(i,:)); % Z-score + end + end + end + end +end \ No newline at end of file diff --git a/Get_Peak_Or_Valley_Indices.m b/Get_Peak_Or_Valley_Indices.m new file mode 100644 index 0000000..ad6a5cd --- /dev/null +++ b/Get_Peak_Or_Valley_Indices.m @@ -0,0 +1,13 @@ +% Get peak or valley indices +% +% by Jesús E. Pérez-Ortega - Feb 2012 + +function [indices,count] = Get_Peak_Or_Valley_Indices(data,threshold,detect_peaks) + if(detect_peaks) + indices = find(data>=threshold); + else + % detect valleys + indices = find(data0; + end + case 'average' + vectors=zeros(peaks,C); + for i=1:peaks + DataPeak_i=data(:,peak_indices==i); + vectors(i,:)=mean(DataPeak_i,2); + end + case 'network' + vectors=zeros(peaks,C*(C-1)/2); + for i=1:peaks + DataPeak_i=data(:,peak_indices==i); + A=Get_Adjacency_From_Raster(Reshape_Raster(DataPeak_i,bin_network),... + connectivity_method); + + % Get significant network + %{ +% extra=250; +% id = find(peak_indices==i); +% id = id(1)-extra:id(end)+extra; +% DataPeak_i = data(:,id); + iterations = 20; + alpha = 0.05; + single_th = true; + shuffle_method = 'time_shift'; + A = Get_Significant_Network_From_Raster(DataPeak_i,bin_network,iterations,... + alpha,connectivity_method,shuffle_method,single_th); + %} + vectors(i,:)=squareform(A,'tovector'); + end + end + else + disp('There are no data peaks!') + end +end + +%{ + if (n_peaks) + switch(Sim_method) + case 'RV' % 2h + peaks=max(PeaksIdx); + + C=size(Data,1); + DataPeaks3=zeros(C,C,peaks); + DataPeaks3_2=zeros(C,C,peaks); + + for i=1:peaks + I=Data(:,PeaksIdx==i); + DataPeaks3(:,:,i)=I*I'; + end + for i=1:peaks + Trace_DataPeaks3(i)=trace(abs(DataPeaks3(:,:,i)^2)); + end + Sim=eye(peaks); + for i=1:(peaks-1) + for j=(i+1):peaks + Sim(i,j) = trace(DataPeaks3(:,:,i)*DataPeaks3(:,:,j))/... + sqrt(Trace_DataPeaks3(i)*Trace_DataPeaks3(j)); + Sim(j,i) = Sim(i,j); + end + end + Sim=(Sim+1)/2; % Normalization + case 'Network Hamming Coactivity' % 5min + peaks=max(PeaksIdx); + for i=1:peaks + P=Data(:,PeaksIdx==i); + A=Get_Adjacency_From_Raster(P>0,'Coactivity'); + A=A/size(P,2); + A(A<1)=0; + As(i,:)=squareform(A,'tovector'); + end + Distance=squareform(pdist(As,'Hamming')); + Sim=1-Distance; + case 'Network Pearson Euclidean' %5min + peaks=max(PeaksIdx); + for i=1:peaks + P=Data(:,PeaksIdx==i); + A=Get_Adjacency_From_Raster(P,'Pearson'); + As(i,:)=squareform(A,'tovector'); + end + Distance=squareform(pdist(As,'Euclidean')); + Sim=1-Distance/max(Distance(:)); + case 'Network Pearson-Size Euclidean' % 5min + peaks=max(PeaksIdx); + for i=1:peaks + P=Data(:,PeaksIdx==i); + A=Get_Adjacency_From_Raster(P,'Pearson')*size(P,2); + A(A<0)=0; + As(i,:)=squareform(A,'tovector'); + end + Distance=squareform(pdist(As,'Euclidean')); + Sim=1-Distance/max(Distance(:)); + case 'Network Pearson-Size Euclidean' % 5min + % Convert to binary if needed + DataPeaks=PeaksVectors_JP(Data,PeaksIdx,Vector_method); + + % Similarity + Distance=squareform(pdist(DataPeaks,Sim_method)); + Sim=1-Distance/max(Distance(:)); % Normalization + otherwise + % Convert to binary if needed + DataPeaks=PeaksVectors_JP(Data,PeaksIdx,Vector_method); + + % Similarity + Distance=squareform(pdist(DataPeaks,Sim_method)); + Sim=1-Distance/max(Distance(:)); % Normalization + end + end + %} diff --git a/Get_Peaks_Sequence_Sorted.m b/Get_Peaks_Sequence_Sorted.m new file mode 100644 index 0000000..a613f25 --- /dev/null +++ b/Get_Peaks_Sequence_Sorted.m @@ -0,0 +1,23 @@ +% Get the states sequence in order of appearence +% +% By Jesús Pérez-Ortega jan-2018 +% modified march 2018 +% modified may 2018 + +function [peaks_sequence_sorted,idx_sort] = Get_Peaks_Sequence_Sorted(peaks_sequence) + % Get peak sequence + idx_group=[]; + groups=max(peaks_sequence); + for i=1:groups + idx_group=[idx_group find(peaks_sequence==i,1,'first')]; + end + [~,idx_sort]=sort(idx_group); + + % Get sequence sorted + peaks_sequence_sorted=peaks_sequence; + j=1; + for i=idx_sort + peaks_sequence_sorted(peaks_sequence==i)=j; + j=j+1; + end +end \ No newline at end of file diff --git a/Get_Peaks_Similarity.m b/Get_Peaks_Similarity.m new file mode 100644 index 0000000..7bc800d --- /dev/null +++ b/Get_Peaks_Similarity.m @@ -0,0 +1,23 @@ +% Similarity Index +% +% Compute similarity between peak vectors +% +% By Jes?s P?rez-Ortega Jan-2018 +% modified Sep-2018 + +function [Sim,Distance] = Get_Peaks_Similarity(vectors,Sim_method) + if(size(vectors,1)>1) + % Compute similarity + Distance=squareform(pdist(vectors,Sim_method)); + if(max(Distance(:))==0) + Sim=ones(length(Distance)); + else + Sim=1-Distance/max(Distance(:)); % Normalization + %Sim=1./(1+Distance); + end + else + warning('There are no data to compare!') + Sim=[]; + Distance=[]; + end +end \ No newline at end of file diff --git a/Get_Peaks_Similarity_OLD.m b/Get_Peaks_Similarity_OLD.m new file mode 100644 index 0000000..da238ba --- /dev/null +++ b/Get_Peaks_Similarity_OLD.m @@ -0,0 +1,84 @@ +% Similarity Index +% +% Get the peaks from coactivity and compute similarity between peaks +% +% By Jesús Pérez-Ortega jan-2018 +% modified march-2018 + +function [Sim PeaksIdx DataPeaks] = Get_Peaks_Similarity_OLD(Data,Smooth_Co,Th,Join,Binary,Sim_method) + + % Find peaks + PeaksIdx=FindPeaks_JP(Smooth_Co,Th,Join); + n_peaks=max(PeaksIdx); + + if (n_peaks) + switch(Sim_method) + case 'RV' % 2h + peaks=max(PeaksIdx); + + C=size(Data,1); + DataPeaks3=zeros(C,C,peaks); + DataPeaks3_2=zeros(C,C,peaks); + + for i=1:peaks + I=Data(:,PeaksIdx==i); + DataPeaks3(:,:,i)=I*I'; + end + for i=1:peaks + Trace_DataPeaks3(i)=trace(abs(DataPeaks3(:,:,i)^2)); + end + Sim=eye(peaks); + for i=1:(peaks-1) + for j=(i+1):peaks + Sim(i,j) = trace(DataPeaks3(:,:,i)*DataPeaks3(:,:,j))/sqrt(Trace_DataPeaks3(i)*Trace_DataPeaks3(j)); + Sim(j,i) = Sim(i,j); + end + end + Sim=(Sim+1)/2; % Normalization + case 'Network Hamming Coactivity' % 5min + peaks=max(PeaksIdx); + for i=1:peaks + P=Data(:,PeaksIdx==i); + A=Get_Adjacency_From_Raster(P>0,'Coactivity'); + A=A/size(P,2); + A(A<1)=0; + As(i,:)=squareform(A,'tovector'); + end + Distance=squareform(pdist(As,'Hamming')); + Sim=1-Distance; + case 'Network Pearson Euclidean' %5min + peaks=max(PeaksIdx); + for i=1:peaks + P=Data(:,PeaksIdx==i); + A=Get_Adjacency_From_Raster(P,'Pearson'); + As(i,:)=squareform(A,'tovector'); + end + Distance=squareform(pdist(As,'Euclidean')); + Sim=1-Distance/max(Distance(:)); + case 'Network Pearson-Size Euclidean' % 5min + peaks=max(PeaksIdx); + for i=1:peaks + P=Data(:,PeaksIdx==i); + A=Get_Adjacency_From_Raster(P,'Pearson')*size(P,2); + A(A<0)=0; + As(i,:)=squareform(A,'tovector'); + end + Distance=squareform(pdist(As,'Euclidean')); + Sim=1-Distance/max(Distance(:)); + case 'Network Pearson-Size Euclidean' % 5min + % Convert to binary if needed + DataPeaks=PeaksVectors_JP(Data,PeaksIdx,Binary); + + % Similarity + Distance=squareform(pdist(DataPeaks,Sim_method)); + Sim=1-Distance/max(Distance(:)); % Normalization + otherwise + % Convert to binary if needed + DataPeaks=PeaksVectors_JP(Data,PeaksIdx,Binary); + + % Similarity + Distance=squareform(pdist(DataPeaks,Sim_method)); + Sim=1-Distance/max(Distance(:)); % Normalization + end + end +end \ No newline at end of file diff --git a/Get_Raster_From_MEA.m b/Get_Raster_From_MEA.m new file mode 100644 index 0000000..dc47d07 --- /dev/null +++ b/Get_Raster_From_MEA.m @@ -0,0 +1,41 @@ +% Build a raster from Multi-Electrode Array (MEA) +% +% Each variable is a matrix with 79 +% columns. The 2nd column indicate the unit and the 3rd column indicate the +% time stamp. +% +% Pérez-Ortega Jesús - Sep 2018 + +function raster = Get_Raster_From_MEA(channel_data,window_ms,omit) + if(nargin==1) + window_ms=1; + omit=[]; + elseif(nargin==2) + omit=[]; + end + + n_channels=unique(channel_data(:,1))'; + + % Get the time length of the raster + max_time=ceil(max(channel_data(:,3)))*1000; + + % Set raster size + raster=zeros(1,max_time); + + % Build raster + k=1; + for i=n_channels + if(~ismember(i,omit)) + channel=channel_data(channel_data(:,1)==i,:); + n_units=max(channel(:,2)); + for j=1:n_units + time_stamps_unit=channel(channel(:,2)==j,3); + spikes_id=round(time_stamps_unit*1000/window_ms)+1; + if ~isnan(spikes_id) + raster(k,spikes_id)=1; + end + k=k+1; + end + end + end +end \ No newline at end of file diff --git a/Get_Raster_From_MEA_Channels.m b/Get_Raster_From_MEA_Channels.m new file mode 100644 index 0000000..c14cc38 --- /dev/null +++ b/Get_Raster_From_MEA_Channels.m @@ -0,0 +1,44 @@ +% Build a raster from Multi-Electrode Array (MEA) +% +% Read channels from the workspace. Each variable is a matrix with 79 +% columns. The 2nd column indicate the unit and the 3rd column indicate the +% time stamp. +% +% Pérez-Ortega Jesús - March 2018 + +function raster = Get_Raster_From_MEA_Channels(window_ms) + if(nargin==0) + window_ms=1; + end + + channel_names=evalin('base','who'); + n_channels=length(channel_names); + + % Get the time length of the raster + max_time=0; + for i=1:n_channels + max_i=evalin('base',['max(' channel_names{i} '(:,3));']); + if(max_i>max_time) + max_time=max_i; + end + end + max_time=ceil(max_time)*1000; + + % Set raster size + raster=zeros(1,max_time); + + % Build raster + k=1; + for i=1:n_channels + channel=evalin('base',[channel_names{i} ';']); + n_units=max(channel(:,2)); + for j=1:n_units + time_stamps_unit=channel(channel(:,2)==j,3); + spikes_id=round(time_stamps_unit*1000/window_ms)+1; + if ~isnan(spikes_id) + raster(k,spikes_id)=1; + end + k=k+1; + end + end +end \ No newline at end of file diff --git a/Get_Raster_From_Time_Stamps.m b/Get_Raster_From_Time_Stamps.m new file mode 100644 index 0000000..31be580 --- /dev/null +++ b/Get_Raster_From_Time_Stamps.m @@ -0,0 +1,49 @@ +% 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. +% +% By Jesús Pérez-Ortega July-2018 + +function raster=Get_Raster_From_Time_Stamps(window,initial_time,final_time) + + if(nargin==1) + partial=false; + elseif(nargin==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); + + % 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)=1; + end + end +end \ No newline at end of file diff --git a/Get_Significant_Network_From_Raster.m b/Get_Significant_Network_From_Raster.m new file mode 100644 index 0000000..daddfce --- /dev/null +++ b/Get_Significant_Network_From_Raster.m @@ -0,0 +1,65 @@ +% Get significant network from raster +% +% Jesús Pérez-Ortega, nov 2018 + +function [A_significant,A_raw,th,As] = Get_Significant_Network_From_Raster(raster,bin,... + iterations,alpha,network_method,shuffle_method,single_th) + + if(nargin == 6) + single_th = false; + elseif(nargin == 5) + single_th = false; + network_method = 'coactivity'; + elseif(nargin == 4) + single_th = false; + shuffle_method = 'time_shift'; + elseif(nargin == 3) + single_th = false; + shuffle_method = 'time_shift'; + network_method = 'coactivity'; + alpha = 0.05; + elseif(nargin == 2) + single_th = false; + shuffle_method = 'time_shift'; + network_method = 'coactivity'; + alpha = 0.05; + iterations = 1000; + elseif(nargin == 1) + single_th = false; + shuffle_method = 'time_shift'; + network_method = 'coactivity'; + alpha = 0.05; + iterations = 1000; + bin=1; + end + + % Reduce raster in bin + raster_bin = Reshape_Raster(raster,bin); + + % Get original adjacency network + A_raw = Get_Adjacency_From_Raster(raster_bin,network_method); + + % Random versions + As = []; + for i = 1:iterations + shuffled = shuffle(raster,shuffle_method); + shuffled_bin = Reshape_Raster(shuffled,bin); + As(i,:) = squareform(Get_Adjacency_From_Raster(shuffled_bin,network_method),... + 'tovector'); + end + + if(single_th) + n_edges = size(As,2); + th = zeros(1,n_edges); + for i = 1:n_edges + th_i = Find_Threshold_In_Cumulative_Distribution(As(:,i),alpha); + th(i) = th_i; + end + th = squareform(th); + else + th = Find_Threshold_In_Cumulative_Distribution(As(:),alpha); + end + + % Get significant adjacency + A_significant = A_raw>th; +end diff --git a/Get_Significant_Network_From_Raster_And_Correlation.m b/Get_Significant_Network_From_Raster_And_Correlation.m new file mode 100644 index 0000000..ee43c0b --- /dev/null +++ b/Get_Significant_Network_From_Raster_And_Correlation.m @@ -0,0 +1,82 @@ +% Get significant network from raster an coactivity +% +% Jesús Pérez-Ortega, nov 2018 + +function [A_significant,A_raw,th] = Get_Significant_Network_From_Raster_And_Correlation(raster,... + correlation,bin,iterations,alpha,network_method,shuffle_method,single_th) + + if(nargin == 6) + single_th = false; + elseif(nargin == 5) + single_th = false; + network_method = 'coactivity'; + elseif(nargin == 4) + single_th = false; + shuffle_method = 'time_shift'; + elseif(nargin == 3) + single_th = false; + shuffle_method = 'time_shift'; + network_method = 'coactivity'; + alpha = 0.05; + elseif(nargin == 2) + single_th = false; + shuffle_method = 'time_shift'; + network_method = 'coactivity'; + alpha = 0.05; + iterations = 1000; + elseif(nargin == 1) + single_th = false; + shuffle_method = 'time_shift'; + network_method = 'coactivity'; + alpha = 0.05; + iterations = 1000; + bin=1; + end + + % Reduce raster in bin + raster_bin = Reshape_Raster(raster,bin); + neurons = size(raster_bin,1); + + % Get original adjacency network + A_raw = Get_Adjacency_From_Raster(raster_bin,network_method); + + % Get adjacency with correlation weights + co(1,:) = correlation; + ws = co'*co; + + A_raw_co = A_raw.*ws; + + % Random versions + As_shuffled = []; + for i = 1:iterations + shuffled = shuffle(raster,shuffle_method); + shuffled_bin = Reshape_Raster(shuffled,bin); + + % Get adjacency + A_shuffled = Get_Adjacency_From_Raster(shuffled_bin,network_method); + + % Get correlation with the coactivity + %ws = reshape(weights(randperm(neurons*neurons)),neurons,neurons); + co = co(randperm(neurons)); + ws = co'*co; + + + A_shuffled_co = A_shuffled.*ws; + As_shuffled(i,:) = squareform(A_shuffled_co,'tovector'); + end + + if(single_th) + n_edges = size(As_shuffled,2); + th = zeros(1,n_edges); + for i = 1:n_edges + th_i = Find_Threshold_In_Cumulative_Distribution(As_shuffled(:,i),alpha); + th(i) = th_i; + end + th = squareform(th); + else + th = Find_Threshold_In_Cumulative_Distribution(As_shuffled(:),alpha); + end + + % Get significant adjacency + A_significant = A_raw_co>th; +end diff --git a/Get_SmallWorld_Properties.m b/Get_SmallWorld_Properties.m new file mode 100644 index 0000000..215a5a3 --- /dev/null +++ b/Get_SmallWorld_Properties.m @@ -0,0 +1,70 @@ +function properties = Get_SmallWorld_Properties(adjacency,only_connected) +% Get the Small-World properties +% +% Pérez-Ortega Jesús E. +% jul 2018 +% modified Sep 18 +% modified Jan 19 + + if(nargin==1) + only_connected=false; + end + + if(only_connected) + nodes_connected=sum(adjacency)>0; + adjacency=adjacency(nodes_connected,nodes_connected); + end + + % Total nodes and total edges + N = length(adjacency); + Links = sum(adjacency); + K = sum(Links)/2; + K_mean = 2*K/N; + Kmax= N*(N-1)/2; + Rho = K/Kmax; + + + % Real + D=distance_bin(adjacency); % distance + [L,E]=charpath(D); % characteristic path length and efficiency + Clocal=clustering_coef_bu(adjacency); % local clustering coefficient + C=mean(Clocal); % clustering coefficient + + % Regular + Reg=Make_Regular_Ring_Network(N,K); + Dreg=distance_bin(Reg); % distance + [Lreg,Ereg]=charpath(Dreg); % characteristic path length and efficiency + Clocalreg=clustering_coef_bu(Reg); % local clustering coefficient + Creg=mean(Clocalreg); % clustering coefficient + + % Red aleatoria + Rand=makerandCIJ_und(N,K); + Drand=distance_bin(Rand); % distance + [Lrand,Erand]=charpath(Drand); % characteristic path length and efficiency + Clocalrand=clustering_coef_bu(Rand); % local clustering coefficient + Crand=mean(Clocalrand); % clustering coefficient + + % Medida de small-world + Omega=Lrand/L-C/Creg; + + properties.N=N; + properties.K=K; + properties.K_mean=K_mean; + properties.Rho=Rho; + properties.Links=Links; + + properties.Clocal=Clocal; + properties.C=C; + properties.L=L; + properties.E=E; + + properties.Creg=Creg; + properties.Lreg=Lreg; + properties.Ereg=Ereg; + + properties.Crand=Crand; + properties.Lrand=Lrand; + properties.Erand=Erand; + + properties.Omega=Omega; +end \ No newline at end of file diff --git a/Get_Spike_Rate.m b/Get_Spike_Rate.m new file mode 100644 index 0000000..f1a147c --- /dev/null +++ b/Get_Spike_Rate.m @@ -0,0 +1,17 @@ +% Get z-score of spike rate of a single neuron +% +% By Jesús Pérez-Ortega aug-2018 + +function spike_rate = Get_Spike_Rate(spikes,bin,step) + n=length(spikes); + n_sum=floor((n-bin)/step+1); + + spike_rate=zeros(1,n); + for i=1:n_sum + % Get spike rate + ini=(i-1)*step+1; + fin=(i-1)*step+bin; + spike_rate(ini:(ini+step-1))=sum(spikes(ini:fin)); + end + spike_rate = spike_rate(1:n); +end \ No newline at end of file diff --git a/Get_Spikes_Shape_From_MEA_Channels.m b/Get_Spikes_Shape_From_MEA_Channels.m new file mode 100644 index 0000000..14b2c4a --- /dev/null +++ b/Get_Spikes_Shape_From_MEA_Channels.m @@ -0,0 +1,52 @@ +% Get the spikes shape from Multi-Electrode Array (MEA) +% +% Read channels from the workspace. Each variable is a matrix with 79 +% columns. The shape of the spikes is from 4-79 to thecolumn indicate the unit and the 3rd column indicate the +% time stamp. +% +% Pérez-Ortega Jesús - March 2018 + +function [mean_shapes min_shapes max_shapes] = Get_Spikes_Shape_From_MEA_Channels(plot_save) + if(nargin==0) + plot_save=1; + end + + channel_names=evalin('base','who'); + n_channels=length(channel_names); + + % Build raster + if(plot_save) + Set_Figure('Spikes shape',[0 0 1000 800]); + end + + k=1; + for i=1:n_channels + channel=evalin('base',[channel_names{i} ';']); + n_units=max(channel(:,2)); + for j=1:n_units + spikes_shape=channel(channel(:,2)==j,4:end); + max_shapes(k,:)=max(spikes_shape); + min_shapes(k,:)=min(spikes_shape); + mean_shapes(k,:)=mean(spikes_shape); + + if(plot_save) + subplot(5,6,k) + mins=min_shapes(k,:); + maxs=max_shapes(k,:)-mins; + + h=area([mins;maxs]',-10000,'LineStyle',':'); hold on + set(h(2),'FaceColor', [0.9 0.9 1]) + set(h(1),'FaceColor', [1 1 1]) + + plot(mean_shapes(k,:),'b') + ylim([min(mins) -min(mins)]) + title(['Channel ' num2str(i) ' - unit ' num2str(j)]) + + if(~mod(k+5,6)) + ylabel('Amplitude (mV)') + end + end + k=k+1; + end + end +end \ No newline at end of file diff --git a/Get_XY_Ensembles.m b/Get_XY_Ensembles.m new file mode 100644 index 0000000..e81c963 --- /dev/null +++ b/Get_XY_Ensembles.m @@ -0,0 +1,77 @@ +function [xy, xy_colors, id, structure] = Get_XY_Ensembles(networks) +% Get coordinate of ensembles +% +% [xy, xy_colors, id, structure] = Get_XY_Ensembles(networks) +% +% Jesus Perez-Ortega April-19 + +% Get number of networks +n_networks = length(networks); +n_neurons = length(networks{1}); +colors = Read_Colors(n_networks); + +% Get identity of neurons +structure = zeros(n_networks,n_neurons); +network = zeros(n_neurons); +for i = 1:n_networks + id_i = sum(networks{i})>0; + if ~isempty(id_i) + structure(i,id_i) = i; + network = network | networks{i}; + end +end + +% Set XY for neurons without connections +id_inactive = find(sum(structure>0)==0); +xy_inactive = Get_Circular_XY(length(id_inactive),3); + +% Set XY for neurons of specificity +alone = sum(structure>0)==1; +xy_base = Get_Circular_XY(n_networks,2); +xy_specific = []; +id_specific = []; +for i = 1:n_networks + % get neurons id + id_i = find(alone & structure(i,:)); + + + if length(id_i)==1 + % Add to the existing + xy_specific = [xy_specific; [0 0]]; + id_specific = [id_specific id_i]; + elseif length(id_i)>1 + % Get network specific-state neurons + net = networks{i}(id_i,id_i); + + % Get coordinates + xy_i = Scale(Get_Force_XY(net),-1,1) + xy_base(i,:); + + % Add to the existing + xy_specific = [xy_specific; xy_i]; + id_specific = [id_specific id_i]; + end +end +% Set XY for hub neurons +id_hubs = find(sum(structure>0)>1); +[~,id_new] = Sort_Raster(flipud(structure(:,id_hubs))','ascend'); +id_hubs = id_hubs(id_new); +xy_hubs = Scale(Get_Force_XY(network(id_hubs,id_hubs)),-1,1); + +% Join +id = [id_specific id_hubs id_inactive]; +xy = [xy_specific; xy_hubs; xy_inactive]; +structure = structure(:,id); + +% Colors +xy_colors = zeros(n_neurons,3); +for i = 1:n_neurons + id_color = find(structure(:,i)); + if isempty(id_color) + xy_colors(i,:) = [0 0 0]; + elseif length(id_color)==1 + xy_colors(i,:) = colors(id_color,:); + else + xy_colors(i,:) = [0.5 0.5 0.5]; + %xy_colors(i,:) = mean(colors(id_color,:)); + end +end \ No newline at end of file diff --git a/Get_Zscore_From_Raster.m b/Get_Zscore_From_Raster.m new file mode 100644 index 0000000..dc541c9 --- /dev/null +++ b/Get_Zscore_From_Raster.m @@ -0,0 +1,20 @@ +% Get z-score of spike rate from raster +% +% By Jesús Pérez-Ortega aug-2018 + +function spike_rate = Get_Zscore_From_Raster(raster,bin,step) + + [c,f]=size(raster); + + spike_rate=zeros([c f]); + for i=1:c + % Get frequencies + spike_rate(i,:)=Get_Spike_Rate(raster(i,:),bin,step); + + % z-score + spike_rate(i,:)=spike_rate(i,:)-mean(spike_rate(i,:)); + if(std(spike_rate(i,:))) + spike_rate(i,:)=spike_rate(i,:)/std(spike_rate(i,:)); + end + end +end \ No newline at end of file diff --git a/Hold_Axes.m b/Hold_Axes.m new file mode 100644 index 0000000..5c64fc5 --- /dev/null +++ b/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/Hold_Figure.m b/Hold_Figure.m new file mode 100644 index 0000000..8230ba8 --- /dev/null +++ b/Hold_Figure.m @@ -0,0 +1,8 @@ +%% Hold on Figure +% +% Jesús Pérez-Ortega March-18 + +function Hold_Figure(TitleName) + h=findobj('name',TitleName); + figure(h); +end \ No newline at end of file diff --git a/Jitter_Raster.m b/Jitter_Raster.m new file mode 100644 index 0000000..aa09161 --- /dev/null +++ b/Jitter_Raster.m @@ -0,0 +1,42 @@ +% Jitter a raster with an specific bin +% +% Jes?s P?rez-Ortega - Sep 2018 + +function jittered = Jitter_Raster(raster,bin,exact_bin) + if(nargin==2) + exact_bin=false; + end + + [n_cells,n_samples]=size(raster); + + jittered=zeros(n_cells,n_samples); + for i=1:n_cells + cell=raster(i,:); + + t=find(cell); + + jit=rand(1,length(t)); + if(exact_bin) + jit(jit>0.5)=bin; + jit(jit<=0.5)=-bin; + else + jit=round(jit*2*bin)-bin; + end + + t_jit=t+jit; + t_jit=sort(t_jit); + + if(max(t_jit)>n_samples) + n_extra=length(find(t_jit>n_samples)); + t_jit=t_jit(1:end-n_extra); + end + + if(min(t_jit)<1) + n_extra=length(find(t_jit<1)); + t_jit=t_jit(1+n_extra:end); + end + + jittered(i,t_jit)=1; + + end +end \ No newline at end of file diff --git a/Join_Peaks.m b/Join_Peaks.m new file mode 100644 index 0000000..a5b3696 --- /dev/null +++ b/Join_Peaks.m @@ -0,0 +1,24 @@ +%% Join Peaks +% Modified jan-2018 +function Peaks = Join_Peaks(StatesIdx,PeaksIdx) + + NS=length(StatesIdx); + NP=length(PeaksIdx); + if NS~=NP + for i=1:NS + idx(i)=find(PeaksIdx==i,1,'first'); + end + else + idx=find(PeaksIdx); + end + + F=length(idx); + Peaks(1)=StatesIdx(1); + j=2; + for i=2:F + if idx(i-1)~=idx(i)-1 + Peaks(j)=StatesIdx(i); + j=j+1; + end + end +end \ No newline at end of file diff --git a/Make_Regular_Ring_Network.m b/Make_Regular_Ring_Network.m new file mode 100644 index 0000000..38dfaa3 --- /dev/null +++ b/Make_Regular_Ring_Network.m @@ -0,0 +1,50 @@ +% Generate Ring Regular Network +% +% This function generates binary undirected network with ring regular +% connections. +% +% regular_network = Make_Regular_Ring_Network(N,K) +% +% Inputs: +% N = number of nodes +% K = number of links +% +% Outputs: +% Reg = binary undirected adjacent matrix of ring regular network +% +% ..:: by Jesús E. Pérez-Ortega ::.. March-2013 +% Modified July-2018 + +function regular_network = Make_Regular_Ring_Network(N,K) + k=2*K/N; + regular_network=zeros(N); + km=round(ceil(k)/2); + for i=1:N + n_b=i-[km:-1:1]; + n_a=i+[1:km]; + + initial=find(n_b<1); + if initial + n_init=length(initial); + n_b(1:n_init)=[N:-1:(N-n_init+1)]; + end + + final=find(n_a>N); + if final + n_final=length(final); + n_a(end:-1:(end-n_final+1))=[1:n_final]; + end + + regular_network(i,[n_b, n_a])=1; + regular_network([n_b, n_a],i)=1; + end + + % Remove excess connections + n_remove=(sum(sum(regular_network))-ceil(k*N))/2; + [i,j]=find(triu(regular_network,2)); + idx=round(rand(1,n_remove)*length(i)); + for l=1:n_remove + regular_network(i(idx(l)),j(idx(l)))=0; + regular_network(j(idx(l)),i(idx(l)))=0; + end +end \ No newline at end of file diff --git a/Neural_Ensemble_Analysis.fig b/Neural_Ensemble_Analysis.fig new file mode 100644 index 0000000..5f240c4 Binary files /dev/null and b/Neural_Ensemble_Analysis.fig differ diff --git a/Neural_Ensemble_Analysis.m b/Neural_Ensemble_Analysis.m new file mode 100644 index 0000000..cd0e1a4 --- /dev/null +++ b/Neural_Ensemble_Analysis.m @@ -0,0 +1,2906 @@ +function varargout = Neural_Ensemble_Analysis(varargin) +% NEURAL_ENSEMBLE_ANALYSIS MATLAB code for Neural_Ensemble_Analysis.fig +% NEURAL_ENSEMBLE_ANALYSIS, by itself, creates a new NEURAL_ENSEMBLE_ANALYSIS +% or raises the existing singleton*. +% +% H = NEURAL_ENSEMBLE_ANALYSIS returns the handle to a new NEURAL_ENSEMBLE_ANALYSIS +% or the handle to the existing singleton*. +% +% NEURAL_ENSEMBLE_ANALYSIS('CALLBACK',hObject,eventData,handles,...) calls the local +% function named CALLBACK in NEURAL_ENSEMBLE_ANALYSIS.M with the given input arguments. +% +% NEURAL_ENSEMBLE_ANALYSIS('Property','Value',...) creates a new NEURAL_ENSEMBLE_ANALYSIS +% or raises the existing singleton*. +% Starting from the left, property value pairs are +% applied to the GUI before Neural_Ensemble_Analysis_OpeningFcn gets called. An +% unrecognized property name or invalid value makes property application +% stop. All inputs are passed to Neural_Ensemble_Analysis_OpeningFcn via varargin. +% +% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one +% instance to run (singleton)". +% +% See also: GUIDE, GUIDATA, GUIHANDLES + +% Edit the above text to modify the response to help Neural_Ensemble_Analysis + +% Last Modified by GUIDE v2.5 26-Apr-2019 13:04:50 + +% Begin initialization code - DO NOT EDIT +gui_Singleton = 1; +gui_State = struct('gui_Name', mfilename, ... + 'gui_Singleton', gui_Singleton, ... + 'gui_OpeningFcn', @Neural_Ensemble_Analysis_OpeningFcn, ... + 'gui_OutputFcn', @Neural_Ensemble_Analysis_OutputFcn, ... + 'gui_LayoutFcn', [] , ... + 'gui_Callback', []); +if nargin && ischar(varargin{1}) + gui_State.gui_Callback = str2func(varargin{1}); +end + +if nargout + [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); +else + gui_mainfcn(gui_State, varargin{:}); +end +% End initialization code - DO NOT EDIT +end + +function Neural_Ensemble_Analysis_OpeningFcn(hObject, eventdata, handles, varargin) + handles.output = hObject; + guidata(hObject, handles); + btnRefreshWorkspace_Callback(hObject, eventdata, handles); +end + +function varargout = Neural_Ensemble_Analysis_OutputFcn(~, ~, handles) + varargout{1} = handles.output; +end + +%% --- Functions --- + +%% Read analysis variable +function Read_Analysis(handles,experiment) + % Raster + % samples per second + samples_per_second=experiment.Raster.SamplesPerSecond; + set(handles.txtSamplesPerSecond,'string',num2str(samples_per_second)); + + % Coactivity + % Smooth filter + set(handles.txtSmoothFilter,'string',num2str(experiment.Coactivity.SmoothFilter)); + % Coactivity Threshold + set(handles.txtCoactivityThreshold,'string',num2str(... + experiment.Coactivity.CoactivityThreshold)); + % Z-score + set(handles.chkZscore,'value',experiment.Coactivity.ZScore); + % Remove oscillations + set(handles.chkRemoveOscillations,'value',experiment.Coactivity.RemoveOscillations); + + % Peaks + set(handles.chkJoin,'value',experiment.Peaks.Join); + set(handles.chkFixedWidth,'value',experiment.Peaks.FixedWidth); + set(handles.rdbSpikes,'value',experiment.Peaks.UseSpikes); + set(handles.rdbFrequencies,'value',~experiment.Peaks.UseSpikes); + set(handles.chkDividePeaks,'value',experiment.Peaks.DividePeaks); + set(handles.txtFixedWidth,'string',num2str(... + experiment.Peaks.FixedWidthWindow*1000/samples_per_second)); + set(handles.txtDividePeaks,'string',num2str(... + experiment.Peaks.DivisionWindow*1000/samples_per_second)); + set(handles.chkDetectPeaksOrValleys,'value',experiment.Peaks.DetectPeaksOrValleys) + if(experiment.Peaks.DetectPeaks) + set(handles.popPeaksValleys,'value',1); + elseif(experiment.Peaks.DetectValleys) + set(handles.popPeaksValleys,'value',2); + else + set(handles.popPeaksValleys,'value',3); + set(handles.popStimuli,'enable','on'); + end + + % Evaluate analysis computed to enable buttons + if(experiment.Peaks.DividePeaks) + set(handles.txtDividePeaks,'enable','on') + set(handles.btnPlotSequenceByPeak,'enable','on') + end + + % Method for create vectors + set(handles.popNetworkMethod,'enable','off'); + switch(experiment.Peaks.VectorMethod) + case 'sum' + num=1; + case 'binary' + num=2; + case 'average' + num=3; + case 'network' + num=4; + set(handles.popNetworkMethod,'enable','on'); + otherwise + num=1; + end + set(handles.popVectorMethod,'Value',num); + + % Method for create networks + switch(experiment.Peaks.NetworkMethod) + case 'coactivity' + num=1; + case 'jaccard' + num=2; + case 'pearson' + num=3; + otherwise + num=1; + end + set(handles.popNetworkMethod,'Value',num); + + % Similarity method + switch(experiment.Clustering.SimilarityMethod) + case 'euclidean' + num=1; + case 'jaccard' + num=2; + case 'cosine' + num=3; + case 'correlation' + num=4; + case 'hamming' + num=5; + case 'seuclidean' + num=6; + case 'cityblock' + num=7; + case 'minkowski' + num=8; + case 'chebychev' + num=9; + case 'mahalanobis' + num=10; + case 'spearman' + num=11; + otherwise + num=1; + end + set(handles.popSimilarityMethod,'Value',num); + + % Linkage method + switch(experiment.Clustering.LinkageMethod) + case 'ward' + num=1; + case 'centroid' + num=2; + case 'median' + num=3; + case 'single' + num=4; + case 'complete' + num=5; + case 'average' + num=6; + case 'weighted' + num=7; + otherwise + num=1; + end + set(handles.popLinkageMethod,'Value',num); + + % Clustering + groups=experiment.Clustering.Groups; + set(handles.txtGroups,'string',num2str(groups)); + if(isfield(experiment.Clustering,'GroupsToPlot')) + groups_to_plot_double = experiment.Clustering.GroupsToPlot; + groups_to_plot=[]; + for i=1:length(groups_to_plot_double) + groups_to_plot=[groups_to_plot num2str(groups_to_plot_double(i)) ',']; + end + else + groups_to_plot=[]; + for i=1:(groups-1) + groups_to_plot=[groups_to_plot num2str(i) ',']; + end + groups_to_plot=[groups_to_plot num2str(groups)]; + end + set(handles.txtGroupsToPlot,'string',groups_to_plot); + + % Plots + % Similarity method + try + switch(experiment.Plot.CurrentSorting) + case 'no sorting' + num=1; + case 'activity' + num=2; + case 'activation' + num=3; + case 'jaccard correlation' + num=4; + case 'linear correlation' + num=5; + otherwise + by_group = strsplit(experiment.Plot.CurrentSorting,' '); + group = str2num(by_group{2}); + num = 5 + group; + end + catch + num = 1; + end + set(handles.popSortingNeurons,'Value',num); + + % Evaluate analysis computed to enable buttons + set(handles.btnDunnTest,'enable','on') + set(handles.btnSaveBySeconds,'enable','on') + + % Vectors + % Normalized Instant Frequencies + if(isfield(experiment.Raster,'Frequencies')) + set(handles.btnGetFrequencies,'ForeGroundColor',[0 0.5 0]) + set(handles.btnPlotFrequencies,'enable','on') + set(handles.rdbFrequencies,'enable','on') + end + if(experiment.Peaks.DetectPeaksOrValleys) + % Enable + set(handles.popPeaksValleys,'enable','on') + set(handles.txtCoactivityThreshold,'enable','on') + set(handles.btnShuffleTest,'enable','on') + set(handles.chkFixedWidth,'enable','on') + if(get(handles.chkFixedWidth,'value')) + set(handles.txtFixedWidth,'enable','on') + end + else + % Disable + set(handles.popPeaksValleys,'enable','off') + set(handles.txtCoactivityThreshold,'enable','off') + set(handles.btnShuffleTest,'enable','off') + set(handles.chkFixedWidth,'enable','off') + % Set values + set(handles.chkFixedWidth,'value',true) + end + + % Coactivity + if(isfield(experiment.Coactivity,'Coactivity')) + set(handles.btnGetCoactivity,'ForeGroundColor',[0 0.5 0]) + set(handles.btnPlotCoactivity,'enable','on') + set(handles.btnGetVectors,'enable','on') + set(handles.popSortingNeurons,'string',{'no sorting','activity','activation',... + 'jaccard correlation','linear correlation'}) + end + + % Peaks + if(isfield(experiment.Peaks,'Vectors')) + set(handles.btnGetVectors,'ForeGroundColor',[0 0.5 0]) + set(handles.btnGetSimilarity,'enable','on') + if(experiment.Peaks.DetectStimuli) + set(handles.btnGetNetworksByStimuli,'enable','on') + end + end + + % Similarity + if(isfield(experiment.Clustering,'Similarity')) + set(handles.btnGetSimilarity,'ForeGroundColor',[0 0.5 0]) + set(handles.btnPlotSimilarity,'enable','on') + set(handles.btnGetClusters,'enable','on') + set(handles.btnDunnTest,'enable','on') + end + + % Clustering + if(isfield(experiment.Clustering,'GroupIndices')) + set(handles.btnGetClusters,'ForeGroundColor',[0 0.5 0]) + set(handles.btnPlotPeaks,'enable','on') + set(handles.btnPlotSequence,'enable','on') + set(handles.btnGetNetworks,'enable','on') + % Enable options on sorting raster/frequencies + labels={'no sorting','activity','activation','jaccard correlation',... + 'linear correlation'}; + for i=1:groups + labels{end+1}=['group ' num2str(i)]; + end + set(handles.popSortingNeurons,'string',labels); + + if(get(handles.popNetworkMethod,'value')) + set(handles.btnGetNetworks,'enable','on') + end + end + + % Network + if(isfield(experiment,'Network')) + if(isfield(experiment.Network,'Network')) + set(handles.btnGetNetworks,'ForeGroundColor',[0 0.5 0]) + set(handles.btnPlotNetworks,'enable','on') + end + if(isfield(experiment.Network,'Stimulus')) + set(handles.btnGetNetworksByStimuli,'ForeGroundColor',[0 0.5 0]) + set(handles.btnPlotStimulusNetworks,'enable','on') + end + if(isfield(experiment.Network,'WholeNetwork')) + set(handles.btnGetNetworkRaster,'ForeGroundColor',[0 0.5 0]) + set(handles.btnPlotWholeNetwork,'enable','on') + end + + end + +end + +%% Read experiment +function experiment=Read_Experiment(handles) + data_num=get(handles.popWorkspace,'Value'); + data_strings=get(handles.popWorkspace,'String'); + name=data_strings{data_num}; + if evalin('base',['exist(''' name '_analysis'')']) + experiment=evalin('base',[name '_analysis']); + else + experiment=[]; + end +end + +%% Write experiment data +function Write_Experiment(handles,experiment) + data_num=get(handles.popWorkspace,'Value'); + data_strings=get(handles.popWorkspace,'String'); + name=data_strings{data_num}; + assignin('base',[name '_analysis'],experiment); +end + +%% Refresh workspace data +function btnRefreshWorkspace_Callback(~, ~, handles) + data_strings = evalin('base','who'); + set(handles.popWorkspace,'String',[{'Select raster'};data_strings]) + set(handles.popWorkspace,'Value',1) + set(handles.popStimuli,'String',[{'select stimuli'};data_strings]) + set(handles.popStimuli,'Value',1) + set(handles.lblRaster,'String','...') + set(handles.lblRaster,'ForegroundColor','black') + Disable_Buttons(handles,'Raster') +end + +%% User number entry control +function user_entry = User_Number_Entry(hObject,default,minimum,maximum,multiple,negative) + if(nargin==4) + multiple=false; + negative=false; + elseif(nargin==5) + negative=false; + end + + user_entry = str2double(get(hObject,'string')); + neg=false; + % verify if the data is a number m?ltiplo + if ~isnan(user_entry) + user_entry=floor(user_entry); + if(user_entry<0) + neg=true; + user_entry=-user_entry; + end + % verify if the data is between the limits + if(user_entry>maximum) + user_entry=maximum; + elseif(user_entrymaximum) + user_entry=maximum; + elseif(user_entrymaximum) + user_entry=maximum; + elseif(user_entry1 && alpha<0) + alpha=0.05; + end + + shuffle_method=str2num(answer{3}); + if(isempty(shuffle_method)) + shuffle_method=1; + end + switch(shuffle_method) + case 1 + shuffle_method='time_shift'; + case 2 + shuffle_method='frames'; + case 3 + shuffle_method='time'; + case 4 + shuffle_method='isi'; + case 5 + shuffle_method='cell'; + case 6 + shuffle_method='exchange'; + otherwise + shuffle_method='time_shift'; + end + + interval_by_cell=str2num(answer{4}); + if(isempty(interval_by_cell)) + interval_by_cell=1; + end + + % Read experiment + experiment=Read_Experiment(handles); + raster=experiment.Raster.Raster; + name=experiment.Raster.Name; + smooth_window=experiment.Coactivity.SmoothWindow; + z_score=experiment.Coactivity.ZScore; + zscore_window=experiment.Coactivity.ZScoreWindow; + only_mean=experiment.Coactivity.ZScoreOnlyMean; + remove_oscillations=experiment.Coactivity.RemoveOscillations; + + % Process + plot=false; + coactivity_threshold_shuffeld_test = Shuffle_Test(raster,smooth_window,... + iterations,shuffle_method,alpha,interval_by_cell,remove_oscillations,... + z_score,only_mean,zscore_window,plot); + + % Save + if(get(handles.chkSavePlot,'value')) + Save_Figure(['Shuffle test (' shuffle_method ') - ' name]); + end + + % Write experiment + experiment.Coactivity.ShuffleMethod=shuffle_method; + experiment.Coactivity.CoactivityThresholdShuffleTest=... + coactivity_threshold_shuffeld_test; + experiment.Coactivity.AlphaShuffleTest=alpha; + Write_Experiment(handles,experiment); + end + % Color black + set(hObject,'ForeGroundColor',[0 0 0]) +end + +%% Get coactivity +function btnGetCoactivity_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + tic + + % Read experiment + experiment=Read_Experiment(handles); + raster=experiment.Raster.Raster; + smooth_window=experiment.Coactivity.SmoothWindow; + z_score=experiment.Coactivity.ZScore; + zscore_window=experiment.Coactivity.ZScoreWindow; + only_mean=experiment.Coactivity.ZScoreOnlyMean; + remove_oscillations=experiment.Coactivity.RemoveOscillations; + + % Get coactivity + coactivity = Get_And_Filter_Coactivity(raster,smooth_window); + + % Remove Oscillations + if(remove_oscillations) + percentage = 0.1; + [coactivity,removed]= Remove_Oscillations(coactivity,percentage); + experiment.Coactivity.CoactivityRemoved=removed; + end + + % Z-score + if(z_score) + coactivity=Z_Score_Coactivity(coactivity,zscore_window,only_mean); + end + + % Write experiment + experiment.Coactivity.Coactivity=coactivity; + experiment.Plot.Correlation = []; + experiment.Plot.JaccardCorrelation = []; + experiment.Plot.CorrelationState = []; + + Write_Experiment(handles,experiment); + + % Refresh coactivity threshold + txtCoactivityThreshold_Callback(handles.txtCoactivityThreshold,[],handles); + + % Enable options on sorting raster/frequencies + set(handles.popSortingNeurons,'string',{'no sorting','activity','activation',... + 'jaccard correlation','linear correlation'}) + set(handles.popSortingNeurons,'value',1) + + + % Disable/Enable buttons + Disable_Buttons(handles,'Peaks') + Enable_Buttons(handles,'Peaks') + + % Color green + set(hObject,'ForeGroundColor',[0 0.5 0]) + time=toc; disp(['Done in ' num2str(time) ' s.']) +end + +%% Get neural vectors +function btnGetVectors_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + tic + + % Read experiment + experiment=Read_Experiment(handles); + coactivity=experiment.Coactivity.Coactivity; + coactivity_threshold=experiment.Coactivity.CoactivityThreshold; + vector_method=experiment.Peaks.VectorMethod; + network_method=experiment.Peaks.NetworkMethod; + fixed_width=experiment.Peaks.FixedWidth; + join=experiment.Peaks.Join; + divide_peaks=experiment.Peaks.DividePeaks; + use_spikes=experiment.Peaks.UseSpikes; + fixed_width_window=experiment.Peaks.FixedWidthWindow; + division_window=experiment.Peaks.DivisionWindow; + samples_per_second=experiment.Raster.SamplesPerSecond; + detect_peaks_or_valleys=experiment.Peaks.DetectPeaksOrValleys; + + if join + % Ask for minimum + prompt = {['Enter the minimum time width in ms to consider peak'... + ' (0 = without restriction):']}; + title = 'Enter parameters'; + dims = [1 50]; + default_input = {'0'}; + answer = inputdlg(prompt,title,dims,default_input); + if(~isempty(answer)) + minimum_width=str2num(answer{1})/1000*samples_per_second; + if(isempty(minimum_width)) + minimum_width=0; + elseif(minimum_width<1) + minimum_width=0; + end + else + minimum_width=0; + end + %} + else + minimum_width=0; + end + + if(use_spikes) + data=experiment.Raster.Raster; + else + data=experiment.Raster.Frequencies; + end + + if(fixed_width) + width=fixed_width_window; + else + width=0; + end + + if detect_peaks_or_valleys + if(get(handles.popPeaksValleys,'value')==3) + stimuli = experiment.Peaks.Stimuli; + detect_peaks = true; + ignore_ini_fin = true; + [vector_indices,vector_widths,vector_amplitudes,vector_ini_fin] = ... + Find_Peaks_Or_Valleys(stimuli,coactivity_threshold,join,detect_peaks,... + minimum_width,width,ignore_ini_fin); + experiment.Peaks.DetectPeaks = false; + experiment.Peaks.DetectValleys = false; + experiment.Peaks.DetectStimuli = true; + else + % Find peaks indexes + detect_peaks = get(handles.popPeaksValleys,'value')==1; + ignore_ini_fin = true; + [vector_indices,vector_widths,vector_amplitudes,vector_ini_fin] = ... + Find_Peaks_Or_Valleys(coactivity,coactivity_threshold,join,detect_peaks,... + minimum_width,width,ignore_ini_fin); + experiment.Peaks.DetectPeaks = true; + experiment.Peaks.DetectValleys = false; + experiment.Peaks.DetectStimuli = false; + end + % Write experiment + experiment.Peaks.Widths = vector_widths; + experiment.Peaks.Amplitudes = vector_amplitudes; + experiment.Peaks.IniFin = vector_ini_fin; + else + % Create vector indices + samples=experiment.Raster.Samples; + n=ceil(samples/width); + vector_indices=zeros(samples,1); + for i=1:n + ini=(i-1)*width+1; + fin=i*width; + if(fin>samples) + fin=samples; + end + vector_indices(ini:fin)=i; + end + + % Write experiment + experiment.Peaks.DetectPeaks=false; + experiment.Peaks.DetectValleys = true; + experiment.Peaks.DetectStimuli = false; + end + + if divide_peaks + vector_indices = Divide_Peaks(vector_indices, division_window); + end + + if max(vector_indices)>0 + % Create a matrix with all vector peaks + bin_network = 1; % bin_network = 25; + vectors=Get_Peak_Vectors(data,vector_indices,vector_method,network_method,bin_network); + count=size(vectors,1); + + % Create raster with only peaks + raster_peaks = data(:,vector_indices>0); + raster_no_peaks = data(:,~vector_indices); + + % Write experiment + experiment.Peaks.Indices=vector_indices; + experiment.Peaks.Vectors=vectors; + experiment.Peaks.Count=count; + experiment.Peaks.RasterPeaks=raster_peaks; + experiment.Peaks.RasterNoPeaks=raster_no_peaks; + Write_Experiment(handles,experiment); + + % Disable/Enable buttons + Disable_Buttons(handles,'Similarity') + Enable_Buttons(handles,'Similarity') + if(get(handles.popPeaksValleys,'value')==3) + set(handles.btnGetNetworksByStimuli,'enable','on') + else + set(handles.btnGetNetworksByStimuli,'enable','off') + end + + % Color green + set(hObject,'ForeGroundColor',[0 0.5 0]) + else + warning('No vectors created') + % Color red + set(hObject,'ForeGroundColor',[0.5 0 0]) + end + time=toc; disp(['Done in ' num2str(time) ' s.']) +end + +%% Get similarity +function btnGetSimilarity_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + tic + + % Read experiment + experiment=Read_Experiment(handles); + vectors=experiment.Peaks.Vectors; + similarity_method=experiment.Clustering.SimilarityMethod; + linkage_method=experiment.Clustering.LinkageMethod; + + % Get similarity + similarity = Get_Peaks_Similarity(vectors,similarity_method); + if(~isempty(similarity)) + tree=linkage(squareform(1-similarity,'tovector'),linkage_method); + else + tree=[]; + end + + % Write experiment + experiment.Clustering.Similarity=similarity; + experiment.Clustering.Tree=tree; + Write_Experiment(handles,experiment); + + % Enable buttons + Enable_Buttons(handles,'Clustering') + + % Color green + set(hObject,'ForeGroundColor',[0 0.5 0]) + time=toc; disp(['Done in ' num2str(time) ' s.']) +end + +%% Groups test +function btnDunnTest_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + + prompt = {'Enter test type (1=Contrast; 2=Dunn; 3=Davies)'}; + title = 'Enter parameters'; + dims = [1 50]; + default_input = {'1'}; + answer = inputdlg(prompt,title,dims,default_input); + if(~isempty(answer)) + test_method=str2num(answer{1}); + if(isempty(test_method)) + test_method=1; + end + + switch(test_method) + case 1 + test_method='Contrast'; + case 2 + test_method='Dunn'; + case 3 + test_method='Davies'; + otherwise + test_method='Contrast'; + end + + % Read experiment + experiment=Read_Experiment(handles); + name=experiment.Raster.Name; + tree=experiment.Clustering.Tree; + similarity=experiment.Clustering.Similarity; + groups=experiment.Clustering.Groups; + + if(groups>30) + groups=2:ceil(groups/10):groups; + else + groups=2:30; + end + + + % Process + obj=Set_Figure(['Clustering test (' name ')'],[0 0 600 300]); + [~,groups_test]=ClustIdx_JP(tree, similarity,test_method,obj,'hierarchical',... + groups); % 'Connectivity', 'Davies', 'Contrast' + + % Save + if(get(handles.chkSavePlot,'value')) + Save_Figure(['Clustering test (' test_method ') - ' name]); + end + + % Write experiment + experiment.Clustering.GroupsTest=groups_test; + Write_Experiment(handles,experiment); + end + + % Color black + set(hObject,'ForeGroundColor',[0 0 0]) +end + +%% Get clusters +function btnGetClusters_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + tic + + % Read experiment + experiment=Read_Experiment(handles); + raster = experiment.Raster.Raster; + indices = experiment.Peaks.Indices; + tree=experiment.Clustering.Tree; + groups=experiment.Clustering.Groups; + + % Process + % Get clusters + group_indices=cluster(tree,'maxclust',groups); + % Get Sequences + [sequence_sorted,sorted_indices] = Get_Peaks_Sequence_Sorted(group_indices); + + % Get raster from states + for i = 1:groups + peaks = find(sequence_sorted==i); + n_peaks = length(peaks); + + raster_state = []; + for j = 1:n_peaks + peak = raster(:,indices==peaks(j)); + raster_state = [raster_state peak]; + end + raster_states{i}=raster_state; + end + + % Write experiment + experiment.Clustering.GroupIndices = group_indices; + experiment.Clustering.SortedIndices = sorted_indices; + experiment.Clustering.SequenceSorted = sequence_sorted; + experiment.Clustering.RasterStates = raster_states; + experiment.Plot.CorrelationState = []; + Write_Experiment(handles,experiment); + + % Enable options on sorting raster/frequencies + labels={'no sorting','activity','activation','jaccard correlation','linear correlation'}; + for i=1:groups + labels{end+1}=['group ' num2str(i)]; + end + set(handles.popSortingNeurons,'string',labels); + set(handles.popSortingNeurons,'value',1); + + % Enable buttons + Enable_Buttons(handles,'PlotStates') + Enable_Buttons(handles,'Network') + + if(experiment.Peaks.DividePeaks) + set(handles.btnPlotSequenceByPeak,'enable','on') + end + + % Color green + set(hObject,'ForeGroundColor',[0 0.5 0]) + time=toc; disp(['Done in ' num2str(time) ' s.']) +end + +%% Get network +function btnGetNetworkRaster_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + tic + + % Read experiment + experiment = Read_Experiment(handles); + raster = experiment.Raster.Raster; + bin = 1; + iterations = 1000; + alpha = 0.05; + network_method = experiment.Peaks.NetworkMethod; + shuffle_method = 'time_shift'; + single_th = true; + + % Get significant network + [whole_network,whole_coactivity] = Get_Significant_Network_From_Raster(... + raster,bin,iterations,alpha,network_method,shuffle_method,single_th); + + % Write experiment + experiment.Network.WholeCoactivity = whole_coactivity; + experiment.Network.WholeNetwork = whole_network; + Write_Experiment(handles,experiment); + + % Enable buttons + set(handles.btnPlotWholeNetwork,'enable','on') + + % Color green + set(hObject,'ForeGroundColor',[0 0.5 0]) + time=toc; disp(['Done in ' num2str(time) ' s.']) +end + +%% Get networks from states +function btnGetNetworks_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + + % Read experiment + experiment = Read_Experiment(handles); + n = experiment.Raster.Neurons; + groups = experiment.Clustering.Groups; + + % clear current networks + experiment.Network.State = {}; + experiment.Network.StateSignificant = {}; + + with_th = true; + if(with_th) + bin = 1; %bin = 25; + iterations = 1000; + alpha = 0.05; + network_method = experiment.Peaks.NetworkMethod; + shuffle_method = 'time_shift'; + single_th = true; + + groups_to_plot = experiment.Clustering.GroupsToPlot; + % + %correlation = experiment.Plot.Correlation; + +% sequence_sorted = experiment.Clustering.SequenceSorted; +% indices = experiment.Peaks.Indices; +% raster = experiment.Raster.Raster; + % Get significant networks of each state + for i=1:groups + if(ismember(i,groups_to_plot)) + + % Get state raster + raster_state = experiment.Clustering.RasterStates{i}; + + % compute correlation + try + correlation = experiment.Plot.CorrelationState(i,:); + catch + correlation = Compute_Correlation_by_State(experiment,i); + experiment.Plot.CorrelationState(i,:)=correlation; + end + + tic + % Add raster with no peks to the raster peaks + %{ + raster_noise = experiment.Peaks.RasterNoPeaks; + n_state=size(raster_state,2); + n_noise=size(raster_noise,2); + if(n_noise>n_state) + id = randperm(n_noise); + id = id(1:n_state); + raster_noise = raster_noise(:,id); + end + raster_state = [raster_state raster_noise]; + %} + % Get raster from states + %{ + extra=500; + peaks = find(sequence_sorted==i); + n_peaks = length(peaks); + raster_state = []; + for j = 1:n_peaks + id = find(indices==peaks(j)); + id = id(1)-extra:id(end)+extra; + peak = raster(:,id); + raster_state = [raster_state peak]; + end + %} + + % Get significant network +% [state_network_th,state_network] = Get_Significant_Network_From_Raster(... +% raster_state,bin,iterations,alpha,network_method,shuffle_method,single_th); + [state_network_th,state_network] = ... + Get_Significant_Network_From_Raster_And_Correlation(raster_state,... + correlation,bin,iterations,alpha,network_method,shuffle_method,single_th); + + % Minimum two coactivations + %state_network_th(state_network < 2) = 0; + + % Write in experiment + experiment.Network.State{i}=state_network; + experiment.Network.StateSignificant{i}=state_network_th; + toc + + end + end + %} + % Get core network +% raster_peaks = experiment.Peaks.RasterPeaks; +% [core_network_th,core_network] = Get_Significant_Network_From_Raster(raster_peaks,... +% bin,iterations,alpha,network_method,shuffle_method,single_th); + + % Network of all networks = union of states and core network = intersection + core_network = ones(n); + network = zeros(n); + for i=1:groups + if(ismember(i,groups_to_plot)) + state_network_th = experiment.Network.StateSignificant{i}; + network = network + state_network_th; + core_network = core_network .* state_network_th; + end + end + core_network_th = core_network; + network_th = logical(network); + else + networks = experiment.Peaks.Vectors; + sequence = experiment.Clustering.SequenceSorted; + % Get significant networks of each state + for i=1:groups + if(ismember(i,groups_to_plot)) + % Get significant network + if(sum(sequence==i)==1) + state_network = squareform(networks(sequence==i,:)>0); + else + state_network = squareform(mean(networks(sequence==i,:)>0)); + end + state_network_th = state_network==1; + + % Write in experiment + experiment.Network.State{i}=state_network; + experiment.Network.StateSignificant{i}=state_network_th; + end + end + + % Network of whole raster = intersection of states + core_network = squareform(mean(networks>0)); + core_network_th = core_network==1; + + % Network of all networks = union of states + network = zeros(n); + for i=1:groups + if(ismember(i,groups_to_plot)) + state_network_th = experiment.Network.StateSignificant{i}; + network = network + state_network_th; + end + end + network_th = logical(network); + end + + % Write experiment + experiment.Network.Network=network; + experiment.Network.Significant=network_th; + experiment.Network.CoreNetwork=core_network; + experiment.Network.CoreSignificant=core_network_th; + Write_Experiment(handles,experiment); + + % Enable buttons + Enable_Buttons(handles,'PlotNetwork') + + % Color green + set(hObject,'ForeGroundColor',[0 0.5 0]) +end + +%% Get networks by stimuli +function btnGetNetworksByStimuli_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + + % Read experiment + experiment = Read_Experiment(handles); + n = experiment.Raster.Neurons; + stimuli = experiment.Peaks.Stimuli; + raster = experiment.Raster.Raster; + + groups = sum(unique(stimuli)>0); + stimuli_sequence = zeros(size(stimuli)); + j=1; + for i = unique(stimuli)' + if(i) + stimuli_sequence(stimuli==i)=j; + j=j+1; + end + end + stimuli_sequence = Delete_Consecutive_Coactivation(stimuli_sequence'); + stimuli_sequence = Get_Peaks_Sequence_Sorted(stimuli_sequence); + stimuli_sequence = stimuli_sequence(stimuli_sequence>0); + + idx = Get_Peak_Or_Valley_Indices(stimuli,0.0001,true); + 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 + + % clear current networks + experiment.Network.Stimulus = {}; + experiment.Network.StimulusSignificant = {}; + + bin = 1; %bin = 25; + iterations = 100; + alpha = 0.05; + network_method = experiment.Peaks.NetworkMethod; + shuffle_method = 'time_shift'; + single_th = true; + + % Get significant networks of each stimulus + for i=1:groups + % Get state raster + peaks = find(stimuli_sequence==i); + n_peaks = length(peaks); + + raster_state = []; + for j = 1:n_peaks + peak = raster(:,indices==peaks(j)); + raster_state = [raster_state peak]; + end + raster_state = experiment.Clustering.RasterStates{i}; + + tic + % Get significant network + [state_network_th,state_network] = Get_Significant_Network_From_Raster(... + raster_state,bin,iterations,alpha,network_method,shuffle_method,single_th); +% % Minimum two coactivations + %state_network_th(state_network < 2) = 0; + + % Write in experiment + experiment.Network.Stimulus{i}=state_network; + experiment.Network.StimulusSignificant{i}=state_network_th; + toc + end + %} + % Network of all networks = union of states and core network = intersection + core_network = ones(n); + network = zeros(n); + for i=1:groups + state_network_th = experiment.Network.StimulusSignificant{i}; + network = network + state_network_th; + core_network = core_network .* state_network_th; + end + core_network_th = core_network; + network_th = logical(network); + + % Write experiment + experiment.Network.UnionStimulus=network; + experiment.Network.UnionSignificantStimulus=network_th; + experiment.Network.IntersectionStimulus=core_network; + experiment.Network.IntersectionSignificantStimulus=core_network_th; + Write_Experiment(handles,experiment); + + % Enable buttons + Enable_Buttons(handles,'PlotNetworkByStimulus') + + % Color green + set(hObject,'ForeGroundColor',[0 0.5 0]) +end + +%% --- Callbacks --- +%% Raster +function popWorkspace_Callback(hObject,~,handles) + % Get data + data_num=get(hObject,'Value'); + set(handles.lblRaster,'String','...') + Disable_Buttons(handles,'Raster'); + if data_num>1 + data_strings=get(hObject,'String'); + name=data_strings{data_num}; + % Check if exist + if evalin('base',['exist(''' name ''',''var'')']) + raster=evalin('base',name); + % Evaluate data + set(handles.lblRaster,'ForegroundColor',[0 0.5 0]) + [cells,samples]=size(raster); + if (samples>1 && cells>1 && length(size(raster))==2) + % Data with less than 500 samples or more than 500 cells + if (samplessamples_per_second*samples) + zscore_window=0; + else + zscore_window=round(zscore_window*samples_per_second); + end + experiment.Coactivity.ZScoreWindow=zscore_window; + experiment.Coactivity.ZScoreOnlyMean=only_mean; + end + end + % Write experiment + experiment.Coactivity.ZScore=z_score; + experiment.Coactivity.CoactivityThreshold=coactivity_threshold; + Write_Experiment(handles,experiment); + + % Disable buttons + Disable_Buttons(handles,'Peaks') + set(handles.btnGetCoactivity,'ForegroundColor','black') +end + +function chkRemoveOscillations_Callback(hObject,~,handles) + % Read experiment + experiment=Read_Experiment(handles); + + % Process + remove_oscillations = get(hObject,'value'); + + % Write experiment + experiment.Coactivity.RemoveOscillations=remove_oscillations; + Write_Experiment(handles,experiment); + + % Disable buttons + Disable_Buttons(handles,'Peaks') + set(handles.btnGetCoactivity,'ForegroundColor','black') +end + +%% Neural vectors + +% Selection +function chkDetectPeaksOrValleys_Callback(hObject,~,handles) + % Read experiment + experiment=Read_Experiment(handles); + + % Process + detect_peaks_or_valleys=get(hObject,'value'); + + % Enable/Disable + if(detect_peaks_or_valleys) + % Enable + set(handles.popPeaksValleys,'enable','on') + set(handles.txtCoactivityThreshold,'enable','on') + set(handles.btnShuffleTest,'enable','on') + set(handles.chkFixedWidth,'enable','on') + if(get(handles.chkFixedWidth,'value')) + set(handles.txtFixedWidth,'enable','on') + end + % Set values + coactivity_threshold=str2double(get(handles.txtCoactivityThreshold,'string')); + experiment.Coactivity.CoactivityThreshold=coactivity_threshold; + else + % Disable + set(handles.popPeaksValleys,'enable','off') + set(handles.txtCoactivityThreshold,'enable','off') + set(handles.btnShuffleTest,'enable','off') + set(handles.chkFixedWidth,'enable','off') + set(handles.txtFixedWidth,'enable','on') + % Set values + set(handles.chkFixedWidth,'value',true) + experiment.Peaks.FixedWidth=true; + experiment.Coactivity.CoactivityThreshold=[]; + end + + % Write experiment + experiment.Peaks.DetectPeaksOrValleys=detect_peaks_or_valleys; + Write_Experiment(handles,experiment); + + % Disable buttons + Disable_Buttons(handles,'Similarity') + set(handles.btnGetVectors,'ForegroundColor','black') +end + +% Peaks or valleys +function popPeaksValleys_Callback(hObject,~,handles) + + % Enable Stimuli signal if selected + if(get(hObject,'value')==3) + set(handles.popStimuli,'enable','on') + else + set(handles.popStimuli,'enable','off') + end + + % Disable buttons + Disable_Buttons(handles,'Similarity') + set(handles.btnGetVectors,'ForegroundColor','black') +end + +% Stimuli +function popStimuli_Callback(hObject,~,handles) + set(hObject,'ForegroundColor','black') + + % Read experiment + experiment=Read_Experiment(handles); + + % Get data + data_num=get(hObject,'Value'); + if data_num>1 + data_strings=get(hObject,'String'); + name=data_strings{data_num}; + % Check if exist + if evalin('base',['exist(''' name ''',''var'')']) + stimuli=evalin('base',name); + % Evaluate data + set(handles.lblRaster,'ForegroundColor',[0 0.5 0]) + if(isvector(stimuli)) + if(length(stimuli)~=experiment.Raster.Samples) + set(hObject,'ForegroundColor','red') + end + else + set(hObject,'ForegroundColor','red') + end + else + btnRefreshWorkspace_Callback([],[],handles) + Disable_Buttons(handles,'Raster'); + end + end + + % Write experiment + experiment.Peaks.Stimuli = stimuli; + Write_Experiment(handles,experiment); + + % Disable buttons + Disable_Buttons(handles,'Similarity') + set(handles.btnGetVectors,'ForegroundColor','black') +end + +% Coactivity Threshold +function txtCoactivityThreshold_Callback(hObject,~,handles) + % Read experiment + experiment=Read_Experiment(handles); + z_score=experiment.Coactivity.ZScore; + coactivity=experiment.Coactivity.Coactivity; + + % Process + if(z_score) + maximum=max(floor(coactivity*10)/10); + else + maximum=experiment.Raster.Neurons; + end + coactivity_threshold = User_Double_Entry(hObject,0,0,maximum); + + % Write experiment + experiment.Coactivity.CoactivityThreshold=coactivity_threshold; + Write_Experiment(handles,experiment); + + % Disable buttons + Disable_Buttons(handles,'Similarity') + set(handles.btnGetVectors,'ForegroundColor','black') +end + +% Use spikes +function rdbSpikes_Callback(~,~,handles) + % Read experiment + experiment=Read_Experiment(handles); + + % Process + if(~get(handles.rdbSpikes,'value')) + set(handles.rdbSpikes,'value',true) + end + set(handles.rdbFrequencies,'value',false) + use_spikes=true; + + % Write experiment + experiment.Peaks.UseSpikes=use_spikes; + Write_Experiment(handles,experiment); + + % Disable buttons + Disable_Buttons(handles,'Similarity') + set(handles.btnGetVectors,'ForegroundColor','black') +end + +% Use frequencies +function rdbFrequencies_Callback(~,~,handles) + % Read experiment + experiment=Read_Experiment(handles); + + % Process + if(~get(handles.rdbFrequencies,'value')) + set(handles.rdbFrequencies,'value',true) + end + set(handles.rdbSpikes,'value',false) + use_spikes=false; + + % Write experiment + experiment.Peaks.UseSpikes=use_spikes; + Write_Experiment(handles,experiment); + + % Disable buttons + Disable_Buttons(handles,'Similarity') + set(handles.btnGetVectors,'ForegroundColor','black') +end + +% Fixed width +function chkFixedWidth_Callback(hObject,~,handles) + % Read experiment + experiment=Read_Experiment(handles); + + % Process + fixed_width=get(hObject,'value'); + + if(fixed_width) + set(handles.txtFixedWidth,'enable','on') + else + set(handles.txtFixedWidth,'enable','off') + end + + % Write experiment + experiment.Peaks.FixedWidth=fixed_width; + Write_Experiment(handles,experiment); + + % Disable buttons + Disable_Buttons(handles,'Similarity') + set(handles.btnGetVectors,'ForegroundColor','black') +end + +% Define fixed with +function txtFixedWidth_Callback(hObject,~,handles) + % Read experiment + experiment = Read_Experiment(handles); + samples = experiment.Raster.Samples; + samples_per_second = experiment.Raster.SamplesPerSecond; + detect_peaks_or_valleys = experiment.Peaks.DetectPeaksOrValleys; + + % Process + minimum = 1000/samples_per_second; + maximum = samples*1000/samples_per_second; + multiple = true; + fixed_width_window = round(User_Number_Entry(hObject,minimum,minimum,maximum,multiple,... + detect_peaks_or_valleys)*samples_per_second/1000); + + % Write experiment + experiment.Peaks.FixedWidthWindow=fixed_width_window; + Write_Experiment(handles,experiment); + + % Disable buttons + Disable_Buttons(handles,'Similarity') + set(handles.btnGetVectors,'ForegroundColor','black') +end + +% Join peaks +function chkJoin_Callback(hObject,~,handles) + % Read experiment + experiment=Read_Experiment(handles); + + % Process + join=get(hObject,'value'); + + if(join) + set(handles.chkDividePeaks,'enable','on') + if(get(handles.chkDividePeaks,'value')) + set(handles.txtDividePeaks,'enable','on') + end + else + set(handles.chkDividePeaks,'enable','off') + set(handles.txtDividePeaks,'enable','off') + end + + % Write experiment + experiment.Peaks.Join=join; + Write_Experiment(handles,experiment); + + % Disable buttons + Disable_Buttons(handles,'Similarity') + set(handles.btnGetVectors,'ForegroundColor','black') +end + +% Divide Peaks +function chkDividePeaks_Callback(hObject,~,handles) + % Read experiment + experiment=Read_Experiment(handles); + + % Process + divide_peaks=get(hObject,'value'); + + if(divide_peaks) + set(handles.txtDividePeaks,'enable','on') + else + set(handles.txtDividePeaks,'enable','off') + set(handles.btnPlotSequenceByPeak,'enable','off') + end + + % Write experiment + experiment.Peaks.DividePeaks=divide_peaks; + Write_Experiment(handles,experiment); + + % Disable buttons + Disable_Buttons(handles,'Similarity') + set(handles.btnGetVectors,'ForegroundColor','black') +end + +% Define division time +function txtDividePeaks_Callback(hObject,~,handles) + % Read experiment + experiment=Read_Experiment(handles); + samples_per_second=experiment.Raster.SamplesPerSecond; + fixed_width_window=experiment.Peaks.FixedWidthWindow; + + % Process + minimum=1000/samples_per_second; + maximum=abs(fixed_width_window)*1000/samples_per_second; + multiple=true; + division_window=User_Number_Entry(hObject,minimum,minimum,maximum,... + multiple)*samples_per_second/1000; + + % Write experiment + experiment.Peaks.DivisionWindow=division_window; + Write_Experiment(handles,experiment); + + % Disable buttons + Disable_Buttons(handles,'Similarity') + set(handles.btnGetVectors,'ForegroundColor','black') +end + + +% Vector method +function popVectorMethod_Callback(hObject,~,handles) + % Get data + data_num=get(hObject,'Value'); + data_strings=get(hObject,'String'); + vector_method=data_strings{data_num}; + + % Disable/enable button + if(data_num<4) + set(handles.popNetworkMethod,'enable','off') + else + set(handles.popNetworkMethod,'enable','on') + end + + % Read experiment + experiment=Read_Experiment(handles); + + % Write experiment + experiment.Peaks.VectorMethod=vector_method; + Write_Experiment(handles,experiment); + + % Disable buttons + Disable_Buttons(handles,'Similarity') + set(handles.btnGetVectors,'ForegroundColor','black') +end +% Network method +function popNetworkMethod_Callback(hObject,~,handles) + % Get data + data_num=get(hObject,'Value'); + data_strings=get(hObject,'String'); + network_method=data_strings{data_num}; + + % Read experiment + experiment=Read_Experiment(handles); + + % Write experiment + experiment.Peaks.NetworkMethod=network_method; + Write_Experiment(handles,experiment); + + % Disable buttons + Disable_Buttons(handles,'Similarity') + set(handles.btnGetVectors,'ForegroundColor','black') +end + +%% Similarity +% Similarity method +function popSimilarityMethod_Callback(hObject,~,handles) + % Get data + data_num=get(hObject,'Value'); + data_strings=get(hObject,'String'); + similarity_method=data_strings{data_num}; + + % Read experiment + experiment=Read_Experiment(handles); + + % Write experiment + experiment.Clustering.SimilarityMethod=similarity_method; + Write_Experiment(handles,experiment); + + % Disable buttons + Disable_Buttons(handles,'Clustering') + set(handles.btnGetSimilarity,'ForegroundColor','black') +end +% Linkage method +function popLinkageMethod_Callback(hObject,~,handles) + % Get data + data_num=get(hObject,'Value'); + data_strings=get(hObject,'String'); + linkage_method=data_strings{data_num}; + + % Read experiment + experiment=Read_Experiment(handles); + + % Write experiment + experiment.Clustering.LinkageMethod=linkage_method; + Write_Experiment(handles,experiment); + + % Disable buttons + Disable_Buttons(handles,'Clustering') + set(handles.btnGetSimilarity,'ForegroundColor','black') +end + +%% Clustering +% Groups +function txtGroups_Callback(hObject,~,handles) + % Read experiment + experiment = Read_Experiment(handles); + count = experiment.Peaks.Count; + + % Process + groups = User_Number_Entry(hObject,5,2,count); + groups_string = []; + for i=1:(groups-1) + groups_string = [groups_string num2str(i) ',']; + end + groups_string = [groups_string num2str(groups)]; + set(handles.txtGroupsToPlot,'string',groups_string) + groups_to_plot = 1:groups; + + % Write experiment + experiment.Clustering.Groups = groups; + experiment.Clustering.GroupsToPlot = groups_to_plot; + Write_Experiment(handles,experiment); + + % Disable buttons + set(handles.btnGetClusters,'ForegroundColor','black') +end + +%% Plots +% Sort neurons +function popSortingNeurons_Callback(hObject, ~, handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + + % Read experiment + experiment=Read_Experiment(handles); + raster = experiment.Raster.Raster; + neurons = experiment.Raster.Neurons; + + % Process + data_num=get(hObject,'Value'); + data_strings=get(hObject,'String'); + sorting_method=data_strings{data_num}; + + % Select sorting method + switch(sorting_method) + case 'no sorting' + cell_indices = 1:neurons; + case 'activity' + activity=experiment.Raster.Activity; + [~,cell_indices]=sort(activity,'descend'); + case 'activation' + [~,cell_indices] = Sort_Raster(raster); + case 'jaccard correlation' + coactivity=experiment.Coactivity.Coactivity; + try + correlation = experiment.Plot.JaccardCorrelation(1,:); + [~, cell_indices]=sort(correlation,'descend'); + catch + disp('Computing correlation...') + tic + for i=1:neurons + D(i)=pdist([coactivity'; raster(i,:)],'jaccard'); + end + correlation=1-D; + [~, cell_indices]=sort(correlation,'descend'); + experiment.Plot.JaccardCorrelation=correlation; + toc + end + case 'linear correlation' + if(get(handles.popPeaksValleys,'value')==3) + stimuli = experiment.Peaks.Stimuli; + + disp('Computing correlation...') + tic + for i=1:neurons + D(i)=pdist([stimuli';raster(i,:)],'correlation'); + end + correlation=1-D; + correlation(isnan(correlation))=0; + correlation(correlation<0)=0; + [~, cell_indices]=sort(correlation,'descend'); + toc + else + try + correlation = experiment.Plot.Correlation(1,:); + [~, cell_indices]=sort(correlation,'descend'); + catch + disp('Computing correlation...') + tic + coactivity=experiment.Coactivity.Coactivity; + for i=1:neurons + D(i)=pdist([coactivity';raster(i,:)],'correlation'); + end + correlation=1-D; + correlation(isnan(correlation))=0; + correlation(correlation<0)=0; + [~, cell_indices]=sort(correlation,'descend'); + toc + end + end + experiment.Plot.Correlation=correlation; + otherwise + by_group=strsplit(sorting_method,' '); + group=str2num(by_group{2}); + try + correlation = experiment.Plot.CorrelationState(group,:); + [~, cell_indices]=sort(correlation,'descend'); + catch + [correlation,cell_indices] = Compute_Correlation_by_State(experiment,group); + experiment.Plot.CorrelationState(group,:)=correlation; + end + end + + % Write experiment + experiment.Plot.CurrentIndices = cell_indices; + experiment.Plot.CurrentSorting = sorting_method; + Write_Experiment(handles,experiment); + + % Color green + set(hObject,'ForeGroundColor',[0 0 0]) +end + +function [correlation,cell_indices] = Compute_Correlation_by_State(experiment,group) + + disp('Computing correlation...') + tic + + coactivity = experiment.Coactivity.Coactivity; + indices = experiment.Peaks.Indices; + group_sequence = experiment.Clustering.SequenceSorted; + neurons = experiment.Raster.Neurons; + raster = experiment.Raster.Raster; + + % find indices of group peaks + idx_group=find(group_sequence==group)'; + idx_indices=zeros(size(indices)); + for i=idx_group + idx_indices(indices==i)=1; + end + coactivity=coactivity.*idx_indices; + + % compute correlation + for i=1:neurons + D(i) = pdist([coactivity';raster(i,:)],'correlation'); + end + correlation = 1-D; + correlation(isnan(correlation)) = 0; + correlation(correlation<0) = 0; + % sort + [~, cell_indices] = sort(correlation,'descend'); + + toc +end + +% Groups to plot +function txtGroupsToPlot_Callback(hObject,~, handles) + % Read experiment + experiment=Read_Experiment(handles); + groups=experiment.Clustering.Groups; + + % Process + % default + default_entry=[]; + for i=1:(groups-1) + default_entry=[default_entry num2str(i) ',']; + end + default_entry=[default_entry num2str(groups)]; + default_groups_to_plot=[1:groups]; + + % validate user entry + user_entry = get(hObject,'string'); + if ~isnan(user_entry) + cell_groups=strsplit(user_entry,','); + n=length(cell_groups); + groups_to_plot=[]; + for i=1:n + group_num=str2num(cell_groups{i}); + if(isempty(group_num)) + user_entry=default_entry; + groups_to_plot=default_groups_to_plot; + break; + end + groups_to_plot=[groups_to_plot group_num]; + end + else + user_entry=default_entry; + groups_to_plot=default_groups_to_plot; + end + set(hObject,'string',num2str(user_entry)); + + % Write experiment + experiment.Clustering.GroupsToPlot=groups_to_plot; + Write_Experiment(handles,experiment); +end + +%% --- Plots --- +%% Plot raster +function btnPlotRaster_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + + % Read experiment + experiment = Read_Experiment(handles); + raster = experiment.Raster.Raster; + name = experiment.Raster.Name; + + try + cell_indices = experiment.Plot.CurrentIndices; + current_sorting = experiment.Plot.CurrentSorting; + catch + neurons = experiment.Raster.Neurons; + cell_indices = 1:neurons; + current_sorting = 'no sorting'; + + % Write experiment + experiment.Plot.CurrentIndices = cell_indices; + experiment.Plot.CurrentSorting = current_sorting; + Write_Experiment(handles,experiment); + end + + % Plot raster + plot_spikes=true; + Plot_Raster(raster(cell_indices,:),name,plot_spikes) + if (strcmp(current_sorting,'no sorting')) + ylabel('neurons') + else + ylabel({'neuorns sorted by'; current_sorting}) + end + + % Save + if(get(handles.chkSavePlot,'value')) + Save_Figure(['Raster - ' name]); + end + + % Color black + set(hObject,'ForeGroundColor',[0 0 0]) +end + +%% Plot frequencies +function btnPlotFrequencies_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + + % Read experiment + experiment = Read_Experiment(handles); + frequencies = experiment.Raster.Frequencies; + name = experiment.Raster.Name; + + try + cell_indices = experiment.Plot.CurrentIndices; + current_sorting = experiment.Plot.CurrentSorting; + catch + neurons = experiment.Raster.Neurons; + cell_indices = 1:neurons; + current_sorting = 'no sorting'; + + % Write experiment + experiment.Plot.CurrentIndices = cell_indices; + experiment.Plot.CurrentSorting = current_sorting; + Write_Experiment(handles,experiment); + end + + % Plot frequencies + plot_spikes=false; + Plot_Raster(frequencies(cell_indices,:),name,plot_spikes) + if (strcmp(current_sorting,'no sorting')) + ylabel('neurons') + else + ylabel({'neuorns sorted by'; current_sorting}) + end + + % Save + if(get(handles.chkSavePlot,'value')) + Save_Figure(['Frequencies - ' name]); + end + + % Color black + set(hObject,'ForeGroundColor',[0 0 0]) +end + +%% Plot coactivity +function btnPlotCoactivity_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + + % Read experiment + experiment=Read_Experiment(handles); + coactivity_threshold=experiment.Coactivity.CoactivityThreshold; + name=experiment.Raster.Name; + samples_per_second=experiment.Raster.SamplesPerSecond; + smooth_coactivity=experiment.Coactivity.Coactivity; + + + if(get(handles.popPeaksValleys,'value')==3) + coactivity_threshold = []; + end + + % Plot coactivity + Plot_Coactivity(smooth_coactivity,name,coactivity_threshold,samples_per_second) + if(experiment.Coactivity.ZScore) + ylabel({'coactivity';'(z-score)'}) + end + if (experiment.Coactivity.SmoothFilter>1000/samples_per_second) + if(experiment.Coactivity.RemoveOscillations) + title(['smooth filter (' num2str(experiment.Coactivity.SmoothFilter)... + ' ms) - slow oscillations removed']); + else + title(['smooth filter (' num2str(experiment.Coactivity.SmoothFilter) ' ms)']); + end + elseif(experiment.Coactivity.RemoveOscillations) + title('slow oscillations removed'); + else + title(''); + end + + % Save + if(get(handles.chkSavePlot,'value')) + Save_Figure(['Coactivity - ' name]); + end + + % Color black + set(hObject,'ForeGroundColor',[0 0 0]) +end + +%% Plot similarity +function btnPlotSimilarity_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + + % Read experiment + experiment=Read_Experiment(handles); + similarity=experiment.Clustering.Similarity; + name=experiment.Raster.Name; + tree=experiment.Clustering.Tree; + similarity_method=experiment.Clustering.SimilarityMethod; + linkage_method=experiment.Clustering.LinkageMethod; + + % Plot Similarity + Plot_Similarity(similarity,name,tree,similarity_method,linkage_method) + + % Save + if(get(handles.chkSavePlot,'value')) + Save_Figure(['Similarity - ' name]); + end + + % Color black + set(hObject,'ForeGroundColor',[0 0 0]) +end + +%% Plot peaks +function btnPlotPeaks_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + + % Read experiment + experiment=Read_Experiment(handles); + name=experiment.Raster.Name; + samples_per_second=experiment.Raster.SamplesPerSecond; + coactivity=experiment.Coactivity.Coactivity; + coactivity_threshold=experiment.Coactivity.CoactivityThreshold; + peak_indices=experiment.Peaks.Indices; + indices=experiment.Clustering.SequenceSorted; + groups_to_plot=experiment.Clustering.GroupsToPlot; + groups=experiment.Clustering.Groups; + + if(get(handles.popPeaksValleys,'value')==3) + coactivity_threshold = []; + end + + % noise_group + noise_group=setdiff(1:groups,groups_to_plot); + + % Plot Peaks + peak_number=false; + Plot_Coactivity(coactivity,name,coactivity_threshold,samples_per_second); + Plot_Peaks_On_Coactivity(coactivity,peak_indices,indices,... + noise_group,samples_per_second,peak_number) + if(experiment.Coactivity.ZScore) + ylabel({'coactivity';'(z-score)'}) + end + if (experiment.Coactivity.SmoothFilter>1000/samples_per_second) + if(experiment.Coactivity.RemoveOscillations) + title(['smooth filter (' num2str(experiment.Coactivity.SmoothFilter)... + ' ms) - slow oscillations removed']); + else + title(['smooth filter (' num2str(experiment.Coactivity.SmoothFilter) ' ms)']); + end + elseif(experiment.Coactivity.RemoveOscillations) + title('slow oscillations removed'); + end + + % Save + if(get(handles.chkSavePlot,'value')) + Save_Figure(['Coactivity Peaks - ' name]); + end + + % Color black + set(hObject,'ForeGroundColor',[0 0 0]) +end + +%% Plot sequence +function btnPlotSequence_Callback(hObject,~,handles) + + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + + % Ask for division + prompt = {'Enter the seconds of windows to divide the raster: (0=all sequence)'}; + title = 'Enter parameters'; + dims = [1 50]; + default_input = {'0'}; + answer = inputdlg(prompt,title,dims,default_input); + if(~isempty(answer)) + window_sec=str2num(answer{1}); + if(isempty(window_sec)) + window_sec=1; + end + + % Read experiment + experiment=Read_Experiment(handles); + name=experiment.Raster.Name; + groups_to_plot=experiment.Clustering.GroupsToPlot; + sequence=experiment.Clustering.SequenceSorted; + groups=experiment.Clustering.Groups; + + % noise_group + noise_group=setdiff(1:groups,groups_to_plot); + + if(window_sec==0) + % Plot Sequence + Plot_States_Sequence(sequence,noise_group,name); + + % Save + if(get(handles.chkSavePlot,'value')) + Save_Figure(['Sequence - ' name]); + end + else + % read indices + indices=experiment.Peaks.Indices; + samples=experiment.Raster.Samples; + samples_per_second=experiment.Raster.SamplesPerSecond; + window=window_sec*samples_per_second; + + n_seqs=ceil(samples/window); + indices_seq=zeros(length(sequence),1); + for i=1:n_seqs + if(i*window>samples) + idx=find(indices(((i-1)*window+1):end)>0,1,'first'); + ini=indices((i-1)*window+idx); + fin=max(indices(((i-1)*window+1):end)); + else + idx=find(indices(((i-1)*window+1):i*window)>0,1,'first'); + ini=indices((i-1)*window+idx); + fin=max(indices(((i-1)*window+1):i*window)); + end + + % Plot Sequence + Plot_States_Sequence(sequence(ini:fin),noise_group,... + [name ' (' num2str(i) '-' num2str(n_seqs) ')'],... + Read_Colors(max(sequence))); + + if(i==1) + Hold_Axes('States count'); + ylims = get(gca,'ylim'); + else + Hold_Axes('States count'); + ylim(ylims) + end + + % Save + if(get(handles.chkSavePlot,'value')) + Save_Figure(['Sequence - ' name '_' num2str(i) '-' num2str(n_seqs)]); + end + indices_seq(ini:fin)=i; + end + % Write experiment + experiment.Plot.SecondsForSequences=window_sec; + experiment.Plot.IndicesForSequences=indices_seq; + Write_Experiment(handles,experiment); + end + end + + % Color black + set(hObject,'ForeGroundColor',[0 0 0]) +end + +%% Plot sequence by peak +function btnPlotSequenceByPeak_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + + % Read experiment + experiment=Read_Experiment(handles); + name=experiment.Raster.Name; + division_window=experiment.Peaks.DivisionWindow; + fixed_width_window=abs(experiment.Peaks.FixedWidthWindow); + sequence=experiment.Clustering.SequenceSorted; + samples_per_second=experiment.Raster.SamplesPerSecond; + + % Peak width + width=ceil(fixed_width_window/division_window); + + % Plot save + save_plot=get(handles.chkSavePlot,'value'); + + % Plot states by width + division_ms=division_window*1000/samples_per_second; + + % Plot state sequences + sequences=Plot_States_By_Width(sequence,width,name,save_plot); + Plot_Sequences(sequences,division_ms,save_plot,'Sequences'); + + % Plot graphs from sequencies + %adjacency_vectors=Plot_Sequence_Graph_By_Width(sequence,width,name,save_plot); + %Plot_Sequences(double(adjacency_vectors>0),division_ms,save_plot,'Adjacencies',... + %flipud(gray(max(adjacency_vectors(:))))); + %Plot_Sequences(adjacency_vectors,division_ms,save_plot,'Adjacencies',... + %flipud(gray(max(adjacency_vectors(:))))); + + % Write experiment + experiment.Plot.Sequences=sequences; + %experiment.Plot.AdjacencyVectors=adjacency_vectors; + Write_Experiment(handles,experiment); + + % Color black + set(hObject,'ForeGroundColor',[0 0 0]) +end + +%% Plot whole network +function btnPlotWholeNetwork_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + + % Read experiment + experiment = Read_Experiment(handles); + name = strrep(experiment.Raster.Name,'_','-'); + n = experiment.Raster.Neurons; + coactivity = experiment.Network.WholeCoactivity; + network = experiment.Network.WholeNetwork; + + try + cell_indices = experiment.Plot.CurrentIndices; + catch + neurons = experiment.Raster.Neurons; + cell_indices = 1:neurons; + current_sorting = 'no sorting'; + + % Write experiment + experiment.Plot.CurrentIndices = cell_indices; + experiment.Plot.CurrentSorting = current_sorting; + Write_Experiment(handles,experiment); + end + + % Get coordinates of network +% xy = Get_Force_XY(coactivity); + xy = Get_Circular_XY(n); + + % Plot network + save_plot=get(handles.chkSavePlot,'value'); + + edge_color = [0.5 0.5 0.5]; + node_color = [0.8 0.8 0.8]; + network_plot = coactivity.*network; + Plot_Adjacencies_And_Network(coactivity(cell_indices,cell_indices),... + network_plot(cell_indices,cell_indices),['All networks (union) - ' name],... + xy,node_color,edge_color,save_plot) + + % Color black + set(hObject,'ForeGroundColor',[0 0 0]) +end + + +%% Plot networks +function btnPlotNetworks_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + + % Read experiment + experiment = Read_Experiment(handles); + name = strrep(experiment.Raster.Name,'_','-'); + n = experiment.Raster.Neurons; + groups = experiment.Clustering.Groups; + groups_to_plot=experiment.Clustering.GroupsToPlot; + network = experiment.Network.Network; + network_th = experiment.Network.Significant; + core_network = experiment.Network.CoreNetwork; + core_network_th = experiment.Network.CoreSignificant; + + try + cell_indices = experiment.Plot.CurrentIndices; + catch + neurons = experiment.Raster.Neurons; + cell_indices = 1:neurons; + current_sorting = 'no sorting'; + + % Write experiment + experiment.Plot.CurrentIndices = cell_indices; + experiment.Plot.CurrentSorting = current_sorting; + Write_Experiment(handles,experiment); + end + + % Get coordinates of network +% xy = Get_Force_XY(network); + xy = Get_Circular_XY(n); + + % Plot network + save_plot=get(handles.chkSavePlot,'value'); + + edge_color = [0.5 0.5 0.5]; + node_color = [0.8 0.8 0.8]; + network_plot = network.*network_th; + Plot_Adjacencies_And_Network(network(cell_indices,cell_indices),... + network_plot(cell_indices,cell_indices),['All networks (union) - ' name],... + xy,node_color,edge_color,save_plot) + lims_x = get(gca,'xlim'); + lims_y = get(gca,'ylim'); + + % Plot core network + network_plot = core_network.*core_network_th; + Plot_Adjacencies_And_Network(core_network(cell_indices,cell_indices),... + network_plot(cell_indices,cell_indices),... + ['Core network (intersection) - ' name],... + xy,node_color,edge_color,save_plot) + xlim(lims_x) + ylim(lims_y) + + % Plot network of each state + colors = Read_Colors(groups); + for i=1:groups + if(ismember(i,groups_to_plot)) + state_network = experiment.Network.State{i}; + state_network_th = experiment.Network.StateSignificant{i}; + + network_plot = state_network.*state_network_th; + Plot_Adjacencies_And_Network(state_network(cell_indices,cell_indices),... + network_plot(cell_indices,cell_indices),... + ['State ' num2str(i) ' network - ' name],... + xy,colors(i,:),edge_color,save_plot) + + xlim(lims_x) + ylim(lims_y) + end + end + %btnPlotRandomNetworks_Callback(hObject,[], handles) + + [~, xy_colors, id, structure] = Get_XY_Ensembles(experiment.Network.StateSignificant); + Plot_Ensembles(network_th(id,id),[],xy_colors,structure,name,save_plot); + + % Color black + set(hObject,'ForeGroundColor',[0 0 0]) +end + +%% Plot networks by stimulus +function btnPlotStimulusNetworks_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + + % Read experiment + experiment = Read_Experiment(handles); + name = strrep(experiment.Raster.Name,'_','-'); + n = experiment.Raster.Neurons; + stimuli = experiment.Peaks.Stimuli; + groups = sum(unique(stimuli)>0); + + network = experiment.Network.UnionStimulus; + network_th = experiment.Network.UnionSignificantStimulus; + core_network = experiment.Network.IntersectionStimulus; + core_network_th = experiment.Network.IntersectionSignificantStimulus; + + try + cell_indices = experiment.Plot.CurrentIndices; + catch + neurons = experiment.Raster.Neurons; + cell_indices = 1:neurons; + current_sorting = 'no sorting'; + + % Write experiment + experiment.Plot.CurrentIndices = cell_indices; + experiment.Plot.CurrentSorting = current_sorting; + Write_Experiment(handles,experiment); + end + + % Get coordinates of network +% xy = Get_Force_XY(network); + xy = Get_Circular_XY(n); + + % Plot network + save_plot=get(handles.chkSavePlot,'value'); + + edge_color = [0.5 0.5 0.5]; + node_color = [0.8 0.8 0.8]; + network_plot = network.*network_th; + Plot_Adjacencies_And_Network(network(cell_indices,cell_indices),... + network_plot(cell_indices,cell_indices),['All networks (union) - ' name ' - stimulus'],... + xy,node_color,edge_color,save_plot) + lims_x = get(gca,'xlim'); + lims_y = get(gca,'ylim'); + + % Plot core network + network_plot = core_network.*core_network_th; + Plot_Adjacencies_And_Network(core_network(cell_indices,cell_indices),... + network_plot(cell_indices,cell_indices),... + ['Core network (intersection) - ' name ' - stimulus'],... + xy,node_color,edge_color,save_plot) + xlim(lims_x) + ylim(lims_y) + + % Plot network of each state + colors = Read_Colors(groups); + for i=1:groups + state_network = experiment.Network.Stimulus{i}; + state_network_th = experiment.Network.StimulusSignificant{i}; + + network_plot = state_network.*state_network_th; + Plot_Adjacencies_And_Network(state_network(cell_indices,cell_indices),... + network_plot(cell_indices,cell_indices),... + ['State ' num2str(i) ' network - ' name],... + xy,colors(i,:),edge_color,save_plot) + + xlim(lims_x) + ylim(lims_y) + end + + % Color black + set(hObject,'ForeGroundColor',[0 0 0]) +end + +function btnPlotRandomNetworks_Callback(hObject, ~, handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + + % Read experiment + experiment = Read_Experiment(handles); + name = strrep(experiment.Raster.Name,'_','-'); + width = experiment.Peaks.FixedWidthWindow; + n = experiment.Raster.Neurons; + vector_method = experiment.Peaks.VectorMethod; + network_method = experiment.Peaks.NetworkMethod; + groups = experiment.Clustering.Groups; + raster_states = experiment.Clustering.RasterStates; + + if(isfield(experiment,'Plot')) + cell_indices = experiment.Plot.CellIndices; + else + cell_indices = 1:experiment.Raster.Neurons; + end + + shuffled_joined = zeros(n); + shuffled_core = ones(n); + for i = 1:groups + raster_state = raster_states{i}; + samples=size(raster_state,2); + raster_shuffled = shuffle(raster_state,'time_shift'); + + % Get vectors from state + n=ceil(samples/width); + vector_indices=zeros(samples,1); + for j=1:n + ini=(j-1)*width+1; + fin=j*width; + if(fin>samples) + fin=samples; + end + vector_indices(ini:fin)=j; + end + % Create a matrix with all vector peaks + shuffled_networks=Get_Peak_Vectors(raster_shuffled,vector_indices,... + vector_method,network_method); + state_shuffled = squareform(mean(shuffled_networks>0)); + state_shuffled_th = state_shuffled==1; + + % Get + shuffled_joined = shuffled_joined + state_shuffled_th; + shuffled_core = shuffled_core.*state_shuffled_th; + + % Write in experiment + experiment.Network.StateShuffled{i}=state_shuffled; + experiment.Network.StateSignificantShuffled{i}=state_shuffled_th; + end + shuffled_joined_th = logical(shuffled_joined); + + % Get coordinates of network + xy = Get_Force_XY(shuffled_joined); + xy = xy(cell_indices,:); + + % Plot networks + save_plot=get(handles.chkSavePlot,'value'); + edge_color = [0.5 0.5 0.5]; + node_color = [0.8 0.8 0.8]; + Plot_Adjacencies_And_Network(shuffled_joined(cell_indices,cell_indices),... + shuffled_joined_th(cell_indices,cell_indices),... + ['Shuffled all networks (union) - ' name],... + xy,node_color,edge_color,save_plot) + lims_x = get(gca,'xlim'); + lims_y = get(gca,'ylim'); + + % Plot core network + Plot_Adjacencies_And_Network(shuffled_core(cell_indices,cell_indices),... + shuffled_core(cell_indices,cell_indices),... + ['Shuffled core network (intersection) - ' name],... + xy,node_color,edge_color,save_plot) + xlim(lims_x) + ylim(lims_y) + + % Plot network of each state + colors = Read_Colors(groups); + for i=1:groups + state_shuffled = experiment.Network.StateShuffled{i}; + state_shuffled_th = experiment.Network.StateSignificantShuffled{i}; + + Plot_Adjacencies_And_Network(state_shuffled(cell_indices,cell_indices),... + state_shuffled_th(cell_indices,cell_indices),... + ['Shuffled state ' num2str(i) ' network - ' name],... + xy,colors(i,:),edge_color,save_plot) + + xlim(lims_x) + ylim(lims_y) + end + + % Write experiment + experiment.Network.NetworkShuffled=shuffled_joined; + experiment.Network.ShuffledSignificant=shuffled_joined_th; + experiment.Network.ShuffledCoreSignificant=shuffled_core; + Write_Experiment(handles,experiment); + + % Color black + set(hObject,'ForeGroundColor',[0 0 0]) +end + +%% --- Save --- + +%% Save raster by x seconds +function btnSaveBySeconds_Callback(hObject,~,handles) + % Color yellow + set(hObject,'ForeGroundColor',[0.5 0.5 0]); pause(0.1); pause on + + % Ask for division + prompt = {'Enter the seconds of windows to divide the raster:'}; + title = 'Enter parameters'; + dims = [1 50]; + default_input = {'60'}; + answer = inputdlg(prompt,title,dims,default_input); + if(~isempty(answer)) + window_sec=str2num(answer{1}); + if(isempty(window_sec)) + window_sec=60; + end + + % Read experiment + experiment=Read_Experiment(handles); + name=experiment.Raster.Name; + samples_per_second=experiment.Raster.SamplesPerSecond; + samples=experiment.Raster.Samples; + final_sec=samples/samples_per_second; + use_spikes=experiment.Peaks.UseSpikes; + + % Check if raster was ploted + h=findobj('name',['Raster (' name ')']); + if isempty(h) + % Plot + if(use_spikes) + btnPlotRaster_Callback([],[],handles); + else + btnPlotFrequencies_Callback([],[],handles); + end + end + + % Check if coactivity was ploted + h=findobj('name',['Coactivity (' name ')']); + i=findobj('tag',['CoactiveAxes' name]); + if (isempty(h) && isempty(i)) + btnPlotPeaks_Callback([],[],handles); + end + + % Save + Save_Raster_By_Windows(name,samples_per_second,window_sec,final_sec) + end + % Color black + set(hObject,'ForeGroundColor',[0 0 0]) +end diff --git a/PeaksVectors_JP_OLD.m b/PeaksVectors_JP_OLD.m new file mode 100644 index 0000000..572451a --- /dev/null +++ b/PeaksVectors_JP_OLD.m @@ -0,0 +1,33 @@ +% Peaks Vectors +% Join the vectors of the same peak. +% +% [XPeaks] = PeaksVectors_JP(X,Xid,mode) +% +% Inputs +% X = data as F x C matrix (F = #frames, C = #cells) +% Xid = Fx1 vector containing the peaks indexes +% sumMode = set 1 if you want to sum the activity or set 0 if you only want +% binary data +% +% Outputs +% XPeaks = data as matrix PxC (P = #peaks) +% +% ..:: by Jesús E. Pérez-Ortega ::.. Jun-2012 +% +% V 2.0 3rd input added: 'sumMode' (binary or sum) jun-2012 +% modified march-2018 + +function [XPeaks] = PeaksVectors_JP_OLD(X,Xid,binary) + +C=size(X,1); +peaks=max(Xid); + +XPeaks=zeros(peaks,C); +for i=1:peaks + peak_n=find(Xid==i); + XPeak=X(:,peak_n); + XPeaks(i,:)=sum(XPeak,2); + if binary + XPeaks(i,:)=XPeaks(i,:)&1; + end +end diff --git a/Plot_Activation_Order.m b/Plot_Activation_Order.m new file mode 100644 index 0000000..dc91963 --- /dev/null +++ b/Plot_Activation_Order.m @@ -0,0 +1,95 @@ +function [id_position,entropy,sequences_sorted] = Plot_Activation_Order(raster,network,window,... + bin,colors,name,save) +% Plot activation order of neurons based on a network +% +% By Pérez-Ortega Jesús, Feb 2019 + +switch(nargin) + case 2 + window = size(raster,2); + bin = 1; + colors = autumn(size(raster,1)); + name = ''; + save = false; + case 3 + bin = 1; + colors = autumn(size(raster,1)); + name = ''; + save = false; + case 4 + colors = autumn(size(raster,1)); + name = ''; + save = false; + case 5 + name = ''; + save = false; + case 6 + save = false; +end + +% Get sequences +n = size(raster,1); +sequences = randperm(n); + +n_seqs = size(raster,2)/window; +for i = 1:n_seqs + ini = (i - 1) * window + 1; + fin = i * window; + r = Reshape_Raster(raster(:,ini:fin),bin); + sequence = Get_Network_Activation_Sequence(r,network); + rest = n-length(sequence); + if (rest>0) + resting = setdiff(1:n,sequence); + resting = resting(randperm(rest)); + sequence = [sequence resting]; + end + sequences(i,1:length(sequence)) = sequence; +end + +% Identify activation position +[sequences_sorted,id_position,probabilities] = Sort_Activation_Position(sequences); + +% Get entropy +for i = 1:n + p = probabilities(:,i); + entropy(i) = Get_Entropy(p); +end +% Normalize entropy +%entropy = entropy/Get_Entropy(ones(1,n+1)); % with "neuron 0" +entropy = entropy/Get_Entropy(ones(1,n)); % without "neuron 0" + +Set_Figure(['Sequences - ' name],[0 0 500 800]); +ax=subplot(6,1,1:4); +imagesc(sequences_sorted); hold on +%colormap(ax,[1 1 1;colors])% with "neuron 0" +colormap(ax,colors) % without "neuron 0" +ylabel('sequence #') +title(name) + +ax=subplot(6,1,5); +imagesc(probabilities,[0 0.5]) +set(gca,'ytick',[1 n+1],'yticklabel',[0 n]) +colormap(ax,flipud(gray(100))) +ylabel('neuron #') + +subplot(6,1,6) +plot(entropy,'o-','color',colors(1,:)) +ylabel('entropy') +xlabel('activation order') +xlim([0.5 n+0.5]) +ylim([0 1]) +if(save) + Save_Figure([name ' - activation order']) +end + +% color neurons +% Set_Figure(['Color neurons - ' name],[0 0 20 800]); +% imagesc((1:n)') +% set(gca,'xtick',[],'ytick',[]) +% colormap(colors) +% for i = 1:n +% text(1,i,num2str(i),'horizontalalignment','center') +% end +% if(save) +% Save_Figure([name ' - color neurons']) +% end \ No newline at end of file diff --git a/Plot_Adjacencies_And_Network.m b/Plot_Adjacencies_And_Network.m new file mode 100644 index 0000000..d0fe538 --- /dev/null +++ b/Plot_Adjacencies_And_Network.m @@ -0,0 +1,55 @@ +% Plot weighted and binary adjacency and its network + +function Plot_Adjacencies_And_Network(adjacency,adjacency_threshold,name,xy,node_color,... + edge_color,save_plot) + + if(nargin==1) + adjacency_threshold = adjacency; + name = 'network'; + xy = []; + node_color = [0.2 0.2 0.2]; + edge_color = [0.5 0.5 0.5]; + save_plot = false; + elseif(nargin==2) + name = 'network'; + xy = []; + node_color = [0.2 0.2 0.2]; + edge_color = [0.5 0.5 0.5]; + save_plot = false; + elseif(nargin==3) + xy = []; + node_color = [0.2 0.2 0.2]; + edge_color = [0.5 0.5 0.5]; + save_plot = false; + elseif(nargin==4) + node_color = [0.2 0.2 0.2]; + edge_color = [0.5 0.5 0.5]; + save_plot = false; + elseif(nargin==5) + edge_color = [0.5 0.5 0.5]; + save_plot = false; + elseif(nargin==6) + save_plot = false; + end + + curved = 0; + size_node = 20; + numbers = false; + + Set_Figure(name,[0 0 1000 300]); + % Plot weighted adjacency matrix + subplot(1,3,1) + imagesc(adjacency); pbaspect([1 1 1]) + title(name) + % Plot significant adjacency matrix + ax=subplot(1,3,2); + imagesc(adjacency_threshold>0); colormap(ax,[1 1 1; 0 0 0]); + pbaspect([1 1 1]) + % Plot weighted network + subplot(1,3,3) + Plot_WU_Network(adjacency_threshold,xy,node_color,curved,edge_color,size_node,numbers) + % Save Plot + if(save_plot) + Save_Figure(name) + end +end \ No newline at end of file diff --git a/Plot_Changes_In_Neurons.m b/Plot_Changes_In_Neurons.m new file mode 100644 index 0000000..9c13c4e --- /dev/null +++ b/Plot_Changes_In_Neurons.m @@ -0,0 +1,19 @@ +% Plot change in each neuron activity +% +% Jesús Pérez-Ortega, nov 2018 + +function Plot_Changes_In_Neurons(change,change_title,name) + if(nargin==2) + name=''; + end + n_neurons=length(change); + plot(find(change>0),change(change>0),'.r','markersize',10); hold on + plot(find(change==0),change(change==0),'.k','markersize',10) + plot(find(change<0),change(change<0),'.b','markersize',10) + plot([0 n_neurons],[mean(change) mean(change)],'--k','markersize',10) + title(name) + xlabel('neuron'); ylabel(change_title) + xlim([1 n_neurons]); %ylim([-7000 7000]) + view([90 90]); + set(gca,'xdir','reverse') +end \ No newline at end of file diff --git a/Plot_Coactivity.m b/Plot_Coactivity.m new file mode 100644 index 0000000..8b26e1e --- /dev/null +++ b/Plot_Coactivity.m @@ -0,0 +1,56 @@ +% Plot Coactivity +% +% P?rez-Ortega Jes?s - March 2018 +% Modified Nov 2018 +function Plot_Coactivity(coactivity,name,threshold,samples_per_second) + if(nargin==1) + name = ''; + threshold = []; + samples_per_second = 1; + elseif(nargin==2) + threshold = []; + samples_per_second = 1; + elseif(nargin==3) + samples_per_second = 1; + end + + h=findobj('name',['Raster (' name ')']); + if isempty(h) + Set_Figure(['Coactivity (' name ')'], [0 0 1220 230]); + Set_Axes(['CoactiveAxes' name],[0 0 1 1]); + else + figure(h); + h=Hold_Axes(['CoactiveAxes' name]); + if(isempty(h)) + Set_Axes(['CoactiveAxes' name],[0 0 1 0.34]); + else + cla; + end + end + F=length(coactivity); + + plot((1:F)/samples_per_second,coactivity,'k');hold on + ylabel('coactivity'); ylim([min(coactivity(3:end-2)) max(coactivity(3:end-2))]) + + if(~isempty(threshold)) + plot([1 F]/samples_per_second,[threshold threshold],'--','color',[0.5 0.5 0.5],'lineWidth',2) + end + + if(F/samples_per_second<30) + set(gca,'box','off','xtick',0:F/samples_per_second,'tag',['CoactiveAxes' name]) + xlabel('Time (s)'); + elseif(F/samples_per_second/60<3) + set(gca,'box','off','xtick',0:10:F/samples_per_second,'tag',['CoactiveAxes' name]) + xlabel('Time (s)'); + else + set(gca,'box','off','xtick',0:60:F/samples_per_second,'xticklabel',0:F/samples_per_second/60,'tag',['CoactiveAxes' name]) + xlabel('Time (min)'); + end + + xlim([0 F/samples_per_second]) + +% ylims=get(gca,'ylim'); +% if(ylims(1)<0) +% ylim([0 ylims(2)]) +% end +end \ No newline at end of file diff --git a/Plot_Degree_Distribution.m b/Plot_Degree_Distribution.m new file mode 100644 index 0000000..8177b69 --- /dev/null +++ b/Plot_Degree_Distribution.m @@ -0,0 +1,76 @@ +% Plot degree distribution +% +% Jesús Pérez-Ortega +% modified nov 2018 + +function [slope,R2,intercept] = Plot_Degree_Distribution(adjacency,name,save) + + if(nargin==2) + save=false; + end + + cells=length(adjacency); + links=sum(adjacency,1); + + if ~links + warning('There are no links.') + return; + end + + % Get links histogram (only neurons with at least one link) + links_1=links(links>0); + N=length(links_1); + nbins = round(log2(N)+1); % Sturges rule + %nbins = 2*iqr(links_1)*N^(-1/3); % Freedman-Diaconis rule for heavy tailed + + [y,x]=hist(links_1,nbins); + y=y/sum(y); + + % Plot power law fit + try + [slope, intercept, R2] = Fit_Power_Law(x',y'); + catch + slope=NaN; + intercept=NaN; + R2=NaN; + return + end + xfit=min(x):(max(x)-min(x))/100:max(x); + yfit=intercept*xfit.^(slope); + + % Plot Links + Set_Figure(['Links - ' name],[0 0 600 200]); + Set_Axes(['Links ' name],[0 0 1 1]) + plot(find(links>0),links_1,'ob'); hold on + plot(find(links==0),links(links==0),'or') + set(gca,'xlim',[0 cells],'ylim',[0 max(links)+1]) + title(name); xlabel('# neuron'); ylabel('Count of links') + + % Save + if(save) + Save_Figure([name ' - Links' ]); + end + + % Plot linear links histogram + Set_Figure(['Distribution - ' name],[0 0 600 200]); + Set_Axes(['Links ' name],[0 0 0.5 1]) + plot(x,y,'ob','markersize',10); hold on + plot(xfit,yfit,'k') + xlabel('k'); ylabel('P(k)'); title(['Degree distribution - ' name]) + + % Write the coefficients + text(max(x)/2,max(y)*0.8,['\lambda=' num2str(-slope,'%0.1f')],'FontSize',14,'FontWeight','bold') + text(max(x)/2,max(y)*0.6,['R^2=' num2str(R2,'%0.3f')],'FontSize',14,'FontWeight','bold') + + % Plot loglog links histogram + Set_Axes(['HistLinksPL ' name],[0.5 0 0.5 1]) + loglog(x,y,'ob','markersize',10); hold on + loglog(xfit,yfit,'k') + xlabel('k'); ylabel('P(k)'); title(['Degree distribution - ' name]) + + % Save + if(save) + Save_Figure([name ' - Distribution']); + end + +end \ No newline at end of file diff --git a/Plot_Edge.m b/Plot_Edge.m new file mode 100644 index 0000000..e5b7860 --- /dev/null +++ b/Plot_Edge.m @@ -0,0 +1,84 @@ +%% Plot edge with or without arrow +% +% Jesús Pérez-Ortega - june 2018 +% Modified nov 2018 + +function [x,y] = Plot_Edge(XY_initial,XY_final,radius,length_arrow,line_width,color) + + if(nargin==5) + color=[0 0 0]; + elseif(nargin==4) + color=[0 0 0]; + line_width=2; + elseif(nargin==3) + color=[0 0 0]; + line_width=2; + length_arrow=3.5; + elseif(nargin==2) + color=[0 0 0]; + line_width=2; + length_arrow=3.5; + radius=0.15; + end + + if(XY_initial(1)>0 && XY_final(1)>0 && XY_initial(2)>0 && XY_final(2)<0) + XY_temporal=XY_initial; + XY_initial=XY_final; + XY_final=XY_temporal; + end + + % If loop + if (XY_initial==XY_final) + r=0.2; + [theta,rho] = cart2pol(XY_initial(1),XY_initial(2)); + rho=rho+r; + [x_center,y_center] = pol2cart(theta,rho); + + % Circle + ang=0:0.01:2*pi; + x=r*cos(ang)+x_center; + y=r*sin(ang)+y_center; + else + if(radius) + % Curve + l=norm(XY_initial-XY_final); % length of line + dx=XY_final(1)-XY_initial(1); + dy=XY_final(2)-XY_initial(2); + + alpha=atan2(dy,dx); % angle of rotation + cosa=cos(alpha); + sina=sin(alpha); + + points=linspace(pi/4,3*pi/4); + a=(0.5-cos(points)/2^0.5)*l; + b=((sin(points)-2^0.5/2)/(1-2^0.5/2))*radius; + + x=a*cosa-b*sina+XY_initial(1); + y=a*sina+b*cosa+XY_initial(2); + else + x=[XY_initial(1) XY_final(1)]; + y=[XY_initial(2) XY_final(2)]; + end + end + plot(x,y,'-','color',color,'linewidth',line_width); hold on + + if(length_arrow) + % Arrow end + xai=x([end end-20]); + yai=y([end end-20]); + xa1=length_arrow*[0 0.1 0.08 0.1 0]'; + ya1=length_arrow*[0 0.03 0 -0.03 0]'; + dx=diff(xai); + dy=diff(yai); + alpha=atan2(dy,dx); % angle of rotation + cosa=cos(alpha); + sina=sin(alpha); + xa=xa1*cosa-ya1*sina+xai(1); + ya=xa1*sina+ya1*cosa+yai(1); + + if (XY_initial~=XY_final) + fill(xa,ya,color) + plot(xa,ya,'-','color',color) + end + end +end diff --git a/Plot_Ensembles.m b/Plot_Ensembles.m new file mode 100644 index 0000000..b658c40 --- /dev/null +++ b/Plot_Ensembles.m @@ -0,0 +1,65 @@ +function xy = Plot_Ensembles(network,xy,xy_colors,structure,name,save) +% Plot ensembles +% +% xy = Plot_Ensembles(network,xy,xy_colors,structure,name) +% +% Jesus Perez-Ortega April-19 + +switch(nargin) + case 4 + name = ''; + save = false; + case 3 + save = false; +end + +% Get data +[n_networks,n_neurons] = size(structure); + +if isempty(xy) + % Identify inactive + n_inactive = length(find(sum(network)==0)); + + % Generate internal circular XY for inactive + xy_1 = Get_Circular_XY(n_neurons-n_inactive); + + % Generate external circular XY for inactive + xy_2 = Get_Circular_XY(n_inactive,1.1); + xy = [xy_1;xy_2]; +end + +% Plot structure +Set_Figure(['Structure - ' name],[0 0 1200 300]); +pcolor([structure; zeros(1,n_neurons)]) +colormap([1 1 1;Read_Colors(n_networks)]) +yticks(1.5:n_networks+0.5) +yticklabels(1:n_networks) +xlabel('neuron #') +ylabel('ensemble #') +if save + Save_Figure(['Structure - ' name]) +end + +% Plot ensembles +Set_Figure(['Network - ' name],[0 0 600 600]); +Plot_Network(network,'undirected',xy,xy_colors) +title(['Network - ' name]) +if save + Save_Figure(['Network - ' name]) +end + +Set_Figure(['Single networks' name],[0 0 900 900]); +rows = ceil(n_networks/2); +for i = 1:n_networks + id = find(structure(i,:)); + subplot(rows,2,i) + Plot_Network(network(id,id),'undirected',xy(id,:),xy_colors(id,:)) + title(['# ' num2str(i)]) + xlim([-1 1]) + ylim([-1 1]) +end +if save + Save_Figure(['Single networks - ' name]) +end + + diff --git a/Plot_Hierarchy.m b/Plot_Hierarchy.m new file mode 100644 index 0000000..ee96400 --- /dev/null +++ b/Plot_Hierarchy.m @@ -0,0 +1,60 @@ +% Plot clustering coefficient in function of degree C(k) vs k +% +% Jesús Pérez-Ortega +% Modified Jan 2019 + +function [slope,R2,intercept,links,clocal] = Plot_Hierarchy(adjacency,name,save) + switch(nargin) + case 1 + name = ''; + save = false; + case 2 + save= false; + end + + % Number of cells + N=length(adjacency); + + if ~N + links=0; + clocal=0; + return; + end + + links=sum(adjacency)'; + clocal=clustering_coef_bu(adjacency); + l=links; + c=clocal; + + % delete zeros clustering coefficients + x=links(clocal>0); + y=clocal(clocal>0); + + Set_Figure([name ' - Hierarchy'],[0 0 600 200]); + + [slope, intercept, R2] = Fit_Power_Law(x,y); + xfit=min(x):(max(x)-min(x))/100:max(x); + yfit=intercept*xfit.^slope; + + % linear plot + Set_Axes(['Axes - C(k) lineal' name],[0 0 0.5 1]); hold on + plot(xfit,yfit,'k','linewidth',2,'markersize',10); hold on + plot(l,c,'or','linewidth',2,'markersize',10) + xlabel('k'); ylabel('C(k)') + title('Local Clustering Coeficient C(k)') + + % Write the coefficients + text(max(x)/2,max(y)*0.8,['\alpha=' num2str(-slope,'%0.1f')],'FontSize',14,'FontWeight','bold') + text(max(x)/2,max(y)*0.6,['R^2=' num2str(R2,'%0.3f')],'FontSize',14,'FontWeight','bold') + + % loglog plot + Set_Axes(['Axes - C(k) loglog' name],[0.5 0 0.5 1]) + loglog(xfit,yfit,'k','linewidth',2,'markersize',10); hold on + loglog(l,c,'or','linewidth',2,'markersize',10) + xlabel('k'), ylabel('C(k)') + title('Local Clustering Coeficient C(k)') + + if(save) + Save_Figure([name ' - Hierarchy']) + end +end \ No newline at end of file diff --git a/Plot_Network.m b/Plot_Network.m new file mode 100644 index 0000000..0442e96 --- /dev/null +++ b/Plot_Network.m @@ -0,0 +1,111 @@ +function Plot_Network(adjacency,type,xy,xy_colors,link_color,node_size,node_number) +% Plot weighted and undirected network with loops +% +% Plot_Network(adjacency,xy,xy_colors,link_color,node_size,node_number) +% +% Pérez-Ortega Jesús - april 2018 +% Modified Jan 2019 +% Modified April 2019 + +n=length(adjacency); +switch (nargin) + case 6 + node_number=false; + case 5 + node_number=false; + node_size=10; + case 4 + node_number=false; + node_size=10; + link_color=[0.5 0.5 0.5]; + case 3 + node_number=false; + node_size=10; + link_color=[0.5 0.5 0.5]; + xy_colors = Read_Colors(n); + case 2 + node_number=false; + node_size=10; + link_color=[0.5 0.5 0.5]; + xy_colors = Read_Colors(n); + xy = []; + case 1 + node_number=false; + node_size=10; + link_color=[0.5 0.5 0.5]; + xy_colors = Read_Colors(n); + xy = []; + type = 'undirected'; +end + +if(size(xy_colors,1)==1) + xy_colors=repmat(xy_colors,n,1); +end + +lims = []; +if(nargin==1||isempty(xy)) + nodes=length(adjacency); + xy=Get_Circular_XY(nodes); + lims = [-1.5 1.5]; +% xy = Get_Force_XY(adjacency); +end + +C=length(adjacency); + +% Plot edges +curved = 0; +if(sum(adjacency(:))) + line_widths = zeros(C); + line_widths(adjacency>0) = Scale(adjacency(adjacency>0),0.5,5); + hold on + for a=1:C + for b=a:C + if (adjacency(a,b)) + if strcmp(type,'directed') + length_arrow=3+0.5*line_widths(a,b); + else + length_arrow = 0; + end + link_color = mean([xy_colors(a,:);xy_colors(b,:)]); + Plot_Edge(xy(a,:),xy(b,:),curved,length_arrow,line_widths(a,b),... + link_color); + end + end + end +end + +links=sum(adjacency); +if(length(node_size)==1) + if(sum(links) && node_size) + sizes_in=Scale(sum(adjacency),node_size,node_size+40); + else + sizes_in=ones(1,C)*30; + end +else + sizes_in=Scale(node_size,10,50); +end + +% Plot nodes +for i=1:C + if(links(i)) + plot(xy(i,1),xy(i,2),'.k','MarkerSize',sizes_in(i)+5) + plot(xy(i,1),xy(i,2),'.','color',xy_colors(i,:),'MarkerSize',sizes_in(i)) + if(node_number) + %text(xy(i,1)*1.1,xy(i,2)*1.1,num2str(i),'HorizontalAlignment','Center') + text(xy(i,1),xy(i,2),num2str(i),'HorizontalAlignment','left') + end + else + plot(xy(i,1),xy(i,2),'.','color',mean([0.5 0.5 0.5; xy_colors(i,:)]),... + 'MarkerSize',sizes_in(i)*.2) + end +end +set(gca,'xtick',[],'ytick',[],'xcolor',[1 1 1],'ycolor',[1 1 1]) +if(lims) + xlim(lims) + ylim(lims) +else + xlim([min(xy(:,1)) max(xy(:,1))]) + ylim([min(xy(:,2)) max(xy(:,2))]) +end + +pbaspect([1 1 1]) diff --git a/Plot_Peaks_On_Coactivity.m b/Plot_Peaks_On_Coactivity.m new file mode 100644 index 0000000..e3cb86b --- /dev/null +++ b/Plot_Peaks_On_Coactivity.m @@ -0,0 +1,34 @@ +% Plot Peaks +% +% Pérez-Ortega Jesús - March 2018 + +function Plot_Peaks_On_Coactivity(coactivity,peak_indices,group_indices,noise_group,samples_per_second,plot_peak_number,group_colors,width) + n=max(group_indices); + if(nargin==5) + plot_peak_number=false; + width=1; + group_colors=Read_Colors(n); + elseif(nargin==6) + width=1; + group_colors=Read_Colors(n); + elseif(nargin==7) + width=1; + end + + F=length(peak_indices); + peaks=max(peak_indices); + for i=1:peaks + peak=find(peak_indices==i); + group=group_indices(i); + if(~ismember(group,noise_group)) + x=peak(1)-width:peak(end)+width; + x(x<=0)=[]; + x(x>=F)=[]; + plot(x/samples_per_second,coactivity(x),'-','color',group_colors(group,:),'linewidth',1.5) + if(plot_peak_number) + max_peak=max(coactivity(x)); + text(peak(1)/samples_per_second,max_peak,num2str(i)) + end + end + end +end \ No newline at end of file diff --git a/Plot_Raster.m b/Plot_Raster.m new file mode 100644 index 0000000..c4497e2 --- /dev/null +++ b/Plot_Raster.m @@ -0,0 +1,52 @@ +% Plot Raster +% +% Pérez-Ortega Jesús 2018 +% modified Nov 2018 +% modified Mar 2019 +function Plot_Raster(data,name,spikes,reshape) + switch(nargin) + case 1 + name=''; + spikes = true; + reshape = true; + case 2 + spikes = true; + reshape = true; + case 3 + reshape = true; + end + + [C,F]=size(data); + Set_Figure(['Raster (' name ')'],[0 0 1220 460]); + Set_Axes(['RasterAxes' name],[0 0.34 1 0.66]); + axis([0 F 0.5 C+0.5]); box; hold on + if(spikes) + if(reshape) + if(F<6001) + win = 1; + elseif(F<20001) + win = 10; + elseif(F<100001) + win = 20; + else + win = 50; + end + else + win = 1; + end + + if(sum(data(:,1))>floor(sum(data(:,1)))) + imagesc(data); colormap(flipud(gray)) + else + data = Reshape_Raster(data,win); xlim([0 F/win]) + imagesc(data,[0,1]); colormap([1 1 1; 0 0 0]) + end + else +% data=data>0; +% imagesc(data,[0,1]); colormap([1 1 1; 0.8 0 0]) +% imagesc(data,[-0.2,0.2]); Set_Colormap_Blue_White_Red(); + imagesc(data,[-4 4]); Set_Colormap_Blue_White_Red(); + end + set(gca,'XTicklabel','','XTick',[0 C]) + title(strrep(name,'_','-')) +end \ No newline at end of file diff --git a/Plot_Rich_Club.m b/Plot_Rich_Club.m new file mode 100644 index 0000000..ae61149 --- /dev/null +++ b/Plot_Rich_Club.m @@ -0,0 +1,52 @@ +% Plot rich-club ordering +% +% By Jesús Pérez-Ortega, Jan 2019 + +%function [R_index,R_index_rnd] = Plot_Rich_Club(adjacency,name,save) +function [R_index] = Plot_Rich_Club(adjacency,name,save) + + switch(nargin) + case 1 + name = ''; + save = false; + case 2 + save = false; + end + + % Get basic properties + N = length(adjacency); + links = sum(adjacency); + K = sum(links)/2; + K_mean = 2*K/N; + level = round(sqrt(N*K_mean)); + + % Get rich-club index + Rclub = rich_club_bu(adjacency,level); +% Rclub_rnd = rich_club_bu(makerandCIJ_und(N,K),level); + + % 1,000 iterations + for i = 1:1000 + Rclub_rnds(i,:) = rich_club_bu(makerandCIJ_und(N,K),level); + end + Rclub_rnd_avg = nanmean(Rclub_rnds); + + % Get rich-club ordering + R_index = Rclub./Rclub_rnd_avg; + %R_index_rnd = Rclub_rnd./Rclub_rnd_avg; + + Set_Figure([name ' - Rich club'],[0 0 300 200]); + semilogx(R_index,'o','color',[0 0.5 0],'markersize',10); hold on + %semilogx(R_index_rnd,'or','markersize',10) + plot([1 length(R_index)],[1 1],'--k') + %legend({'Real','ER'}) + xlabel('k'), ylabel('\rho_{ran}(k)') + title([ name ' - Rich-club ordering']) + + ylims = get(gca,'ylim'); + ylims(1) = 0; + ylim(ylims) + + if(save) + Save_Figure([ name ' - Rich-club ordering']) + end +end \ No newline at end of file diff --git a/Plot_Sequence_Graph_By_Width.m b/Plot_Sequence_Graph_By_Width.m new file mode 100644 index 0000000..5bd7895 --- /dev/null +++ b/Plot_Sequence_Graph_By_Width.m @@ -0,0 +1,50 @@ +% Plot sequence graph by width +% +% Pérez-Ortega Jesús - April 2018 + +function adjacency_vectors = Plot_Sequence_Graph_By_Width(sequence,width,name,save,group_colors) + n=max(sequence); + if (nargin==4) + group_colors = Read_Colors(n); + elseif(nargin==3) + save=false; + group_colors = Read_Colors(n); + end + + n_seq=length(sequence); + n_div=floor(n_seq/width); + n_nodes=max(sequence); + + sub_fig=1; + for i=1:n_div + if(~mod(i-1,20)) + Set_Figure(['Sequence graph - ' name ' (' num2str(ceil(i/20)) ')' ],[0 0 1400 800]); + sub_fig=1; + end + subplot(4,5,sub_fig) + if(i*width>n_seq) + seq_i=sequence(((i-1)*width+1):end); + else + seq_i=sequence(((i-1)*width+1):i*width); + end + + adjacency = Get_Adjacency_From_Sequence(seq_i); + Plot_WD_Network(adjacency,[],group_colors); + title(['Peak ' num2str(i)]) + sub_fig=sub_fig+1; + + % Convert directed adjacency matrix to vector + adjacency_complete=zeros(n_nodes); + adjacency_complete(1:length(adjacency),1:length(adjacency))=adjacency; + adjacency_vector=adjacency_complete(:); + adjacency_vectors(i,:)=adjacency_vector; + end + + if(save) + n_win=ceil(n_div/20); + for i=1:n_win + Hold_Figure(['Sequence graph - ' name ' (' num2str(i) ')' ]); + Save_Figure(['Sequence graph - ' name ' (' num2str(i) ').png' ]); + end + end +end diff --git a/Plot_Sequences.m b/Plot_Sequences.m new file mode 100644 index 0000000..b82f17a --- /dev/null +++ b/Plot_Sequences.m @@ -0,0 +1,71 @@ +% Plot all sequences in vertical way +% +% Pérez-Ortega Jesús - May 2018 + +function Plot_Sequences(sequences,division_ms,save,name,colors) + n_colors=max(sequences(:)); + default_colors= Read_Colors(n_colors); + if(nargin==2) + save=false; + name='Sequences'; + colors=default_colors; + elseif(nargin==3) + name='Sequences'; + colors=default_colors; + elseif(nargin==4) + colors=default_colors; + end + + [n,n_in_seq]=size(sequences); + + if(division_ms*n_in_seq>=3000) + times=(division_ms:division_ms:division_ms*n_in_seq)/1000; + elseif(division_ms==0) + times=1:n_in_seq; + else + times=division_ms:division_ms:division_ms*n_in_seq; + end + + Set_Figure(name,[0 0 900 200]); + map=winter(n); + for i=1:n + plot(times,sequences(i,:)+1*i,'color',map(i,:));hold on + end + + sequence_mode=mode(sequences); + errors=0; + for i=1:n_in_seq + errors=errors+length(find(sequences(:,i)~=sequence_mode(i))); + end + + if(division_ms*n_in_seq>=3000) + xlabel('time (s)') + elseif(division_ms==0) + xlabel('peak # (t)') + else + xlabel('time (ms)') + end + ylabel('sequence #') + + errors_percetage=errors/(n*n_in_seq)*100; + + title([name ' - ' num2str(errors) ' errors from the mode - '... + num2str(errors_percetage) '%']) + + % to save + if(save) + Save_Figure(name); + end + + % Sequences in image + Set_Figure([name ' - image'],[0 0 1000 800]); + imagesc(sequences) + colormap(colors) + title([name ' - ' num2str(errors) ' errors from the mode - '... + num2str(errors_percetage) '%']) + + % to save + if(save) + Save_Figure([name ' - image']); + end +end \ No newline at end of file diff --git a/Plot_Similarity.m b/Plot_Similarity.m new file mode 100644 index 0000000..d141ca9 --- /dev/null +++ b/Plot_Similarity.m @@ -0,0 +1,31 @@ +% Plot Similarity +% +% Pérez-Ortega Jesús - March 2018 +function Plot_Similarity(sim,name,tree,sim_method,linkage_method) + if(nargin==3) + sim_method=''; + linkage_method=''; + end + + Set_Figure(['Clustering (' name ')'],[0 0 900 300]); + + % Plot similarity in time function + Set_Axes(['SimAxes' name],[0 0 0.33 1]); + imagesc(sim) + set(gca,'YDir','normal') + title([sim_method ' similarity']) + xlabel('# peak (t)') + + % Plot dendrogram + Set_Axes(['TreeAxes' name],[0.66 0 0.33 1]); + dendrogram(tree,0,'orientation','rigth','ColorThreshold','default'); + Tid = str2num(get(gca,'YTicklabel')); + set(gca,'xtick',[],'ytick',[]) + title([linkage_method ' linkage']) + + % Plot similarity sort by dendrogram + Set_Axes(['SimSortAxes' name],[0.33 0 0.33 1]); + imagesc(sim(Tid,Tid)) + set(gca,'YDir','normal','xtick',[],'ytick',[]) + title('Sort by similarity') +end \ No newline at end of file diff --git a/Plot_Similarity_With_Dendrogram.m b/Plot_Similarity_With_Dendrogram.m new file mode 100644 index 0000000..467b629 --- /dev/null +++ b/Plot_Similarity_With_Dendrogram.m @@ -0,0 +1,39 @@ +% Plot Similarity +% +% Pérez-Ortega Jesús - Nov 2018 +function Plot_Similarity_With_Dendrogram(data,name,sim_method,linkage_method) + if(nargin==1) + name=''; + sim_method='euclidean'; + linkage_method='ward'; + end + + [similarity, distance]= Get_Peaks_Similarity(data,sim_method); + if(~isempty(similarity)) + tree=linkage(squareform(distance,'tovector'),linkage_method); + else + error('Similarity matrix could not be created.') + end + + Set_Figure(['Clustering (' name ')'],[0 0 900 300]); + + % Plot similarity in time function + Set_Axes(['SimAxes' name],[0 0 0.33 1]); + imagesc(similarity) + set(gca,'YDir','normal') + title([sim_method ' similarity']) + xlabel('# peak (t)') + + % Plot dendrogram + Set_Axes(['TreeAxes' name],[0.66 0 0.33 1]); + dendrogram(tree,0,'orientation','rigth','ColorThreshold','default'); + Tid=str2num(get(gca,'YTicklabel')); + %set(gca,'xtick',[],'ytick',[]) + title([linkage_method ' linkage']) + + % Plot similarity sort by dendrogram + Set_Axes(['SimSortAxes' name],[0.33 0 0.33 1]); + imagesc(similarity(Tid,Tid)) + set(gca,'YDir','normal','xtick',[],'ytick',[]) + title('Sort by similarity') +end \ No newline at end of file diff --git a/Plot_Spikes.m b/Plot_Spikes.m new file mode 100644 index 0000000..8a96713 --- /dev/null +++ b/Plot_Spikes.m @@ -0,0 +1,35 @@ +% Plot spikes +% +% Jesús Pérez-Ortega, Nov 2018 + +function Plot_Spikes(spikes,value,color) + switch(nargin) + case 1 + color = [0 0 0]; + value = 1; + case 2 + color = [0 0 0]; + end + + samples = length(spikes); + + id = find(spikes); + n_spikes = length(id); + + time=repmat(id',1,2); + + spikes_1=ones(n_spikes,1)*value-1; + spikes_2=ones(n_spikes,1)*value; + spikes=[spikes_1 spikes_2]; + + plot(time',spikes','color',color) + if(samples<1000) + label = [num2str(samples) ' ms']; + elseif(samples<=60000) + label = [num2str(samples/1000) ' s']; + else + label = [num2str(samples/60000) ' min']; + end + + set(gca,'ytick',[],'xtick',[0 samples],'xticklabel',{'0',label}) +end \ No newline at end of file diff --git a/Plot_Spikes_From_ISI.m b/Plot_Spikes_From_ISI.m new file mode 100644 index 0000000..fa5b14d --- /dev/null +++ b/Plot_Spikes_From_ISI.m @@ -0,0 +1,24 @@ +% Plot spikes from ISI +% +% Jesús Pérez-Ortega, Nov 2018 + +function Plot_Spikes_From_ISI(isi) + + n_spikes = length(isi); + + % Get spikes from ISI + if(size(isi,2)==1) + time=repmat(cumsum(isi),1,2); + else + time=repmat(cumsum(isi)',1,2); + end + + spikes_1=zeros(n_spikes,1); + spikes_2=ones(n_spikes,1); + spikes=[spikes_1 spikes_2]; + + plot(time',spikes','k') +% title(['Spikes = ' num2str(n_spikes)]) +% xlabel('time [s]') + set(gca,'ytick',[]) +end \ No newline at end of file diff --git a/Plot_States_By_Width.m b/Plot_States_By_Width.m new file mode 100644 index 0000000..a4ebdd3 --- /dev/null +++ b/Plot_States_By_Width.m @@ -0,0 +1,53 @@ +% Plot Peak states +% +% Pérez-Ortega Jesús - March 2018 + +function seqs = Plot_States_By_Width(sequence,width,name,save,group_colors) + n=max(sequence); + if (nargin==4) + group_colors = Read_Colors(n); + elseif(nargin==3) + save=false; + group_colors = Read_Colors(n); + end + + states=max(sequence); + n_seq=length(sequence); + n_div=floor(n_seq/width); + + sub_fig=1; + for i=1:n_div + if(~mod(i-1,20)) + Set_Figure(['Sequences by width - ' name ' (' num2str(ceil(i/20)) ')' ],[0 0 1400 800]); + sub_fig=1; + end + subplot(4,5,sub_fig) + if(i*width>n_seq) + seq_i=sequence(((i-1)*width+1):end); + seqs(i,1:length(seq_i))=seq_i; + else + seq_i=sequence(((i-1)*width+1):i*width); + seqs(i,:)=seq_i; + end + + plot(seq_i,'-k');hold on + for j=1:states + idx=find(seq_i==j); + ns=length(idx); + plot(idx,repmat(j,ns,1),'o','color',group_colors(j,:),'markersize',10,... + 'linewidth',2) + end + ylim([1 states]) + title(['Sequence ' num2str(i)]) + + sub_fig=sub_fig+1; + end + + if(save) + n_win=ceil(n_div/20); + for i=1:n_win + Hold_Figure(['Sequences by width - ' name ' (' num2str(i) ')']); + Save_Figure(['Sequences by width - ' name ' (' num2str(i) ')']); + end + end +end diff --git a/Plot_States_Sequence.m b/Plot_States_Sequence.m new file mode 100644 index 0000000..047b769 --- /dev/null +++ b/Plot_States_Sequence.m @@ -0,0 +1,54 @@ +% Plot sequence +% +% P?rez-Ortega Jes?s - March 2018 +% modified Aug-2018 +function count_states = Plot_States_Sequence(sequence,noise_group,name,group_colors) + n=max(sequence); + default_colors= Read_Colors(n); + switch(nargin) + case 1 + noise_group=0; + name='untitled'; + group_colors=default_colors; + case 2 + name='untitled'; + group_colors=default_colors; + case 3 + group_colors=default_colors; + end + groups=max(sequence); + + Set_Figure(['States Sequence - ' name],[0 0 1220 200]); + Set_Axes('States sequence',[-.05 0 .85 1]); hold on + for i=noise_group + sequence=sequence(sequence~=i); + end + for i=1:max(sequence) + idx=find(sequence==i); + ns=length(idx); + count_states(i)=ns; + plot(idx,repmat(i,ns,1),'o','color',group_colors(i,:),'markersize',10,... + 'linewidth',2) + end + plot(sequence,'-k') + ylim([0 groups+1]); ylabel('State') + l=length(sequence); + xlim([0 l+1]); xlabel('# peak') + title(['States sequence - ' strrep(name,'_','-')]) + + Set_Axes('States count',[0.71 0.1 .15 .85]); hold on + for i=1:max(sequence) + bar(i,count_states(i),'FaceColor',group_colors(i,:)); hold on + end + set(gca,'xtick',[]); ylabel('count'); + view(90,-90) + + n=length(sequence)-1; + transitions=sum(diff(sequence)~=0)/n; + same=1-transitions; + Set_Axes('transitions',[0.86 0 .14 1]); hold on + bar([1 2],[same transitions]); hold on + set(gca,'xtick',1:2,'xticklabel',{'same','between'}); + xlim([.5 2.5]); ylim([0 1]) + title('transitions'); ylabel('fraction'); +end \ No newline at end of file diff --git a/Plot_WD_Network.m b/Plot_WD_Network.m new file mode 100644 index 0000000..f26da66 --- /dev/null +++ b/Plot_WD_Network.m @@ -0,0 +1,70 @@ +% Plot weighted and directed network with loops +% +% Pérez-Ortega Jesús - april 2018 +% modified june 2018 + +function Plot_WD_Network(adjacency,xy,xy_colors,node_size,node_number) + n=length(adjacency); + if (nargin==4) + node_number=false; + elseif (nargin==3) + node_number=false; + node_size=35; + elseif (nargin==2) + node_number=false; + node_size=35; + xy_colors = Read_Colors(n); + elseif (nargin==1) + node_number=false; + node_size=35; + xy_colors = Read_Colors(n); + end + + if(nargin==1||isempty(xy)) + nodes=length(adjacency); + xy=Get_Circular_XY(nodes); + end + + % Plot edges with arrow + C=length(adjacency); + size_factor=5/sum(adjacency(:)); + radius=.15; % Arrow property + hold on + for a=1:C + for b=a:C + if (adjacency(a,b)) + line_width=adjacency(a,b)*size_factor+.5; + length_arrow=3+0.5*line_width; + Plot_Edge(xy(a,:),xy(b,:),radius,length_arrow,line_width); + end + end + end + for a=1:C-1 + for b=a+1:C + if (adjacency(b,a)) + line_width=adjacency(b,a)*size_factor+.5; + length_arrow=3+0.5*line_width; + Plot_Edge(xy(b,:),xy(a,:),radius,length_arrow,line_width); + end + end + end + if(sum(adjacency(:))) + sizes_in=sum(adjacency)./sum(adjacency(:)).*35+node_size; + else + sizes_in=ones(1,C)*.35+node_size; + end + + % Plot nodes + for i=1:C + plot(xy(i,1),xy(i,2),'.k','MarkerSize',sizes_in(i)) + plot(xy(i,1),xy(i,2),'.','color',xy_colors(i,:),'MarkerSize',sizes_in(i)*.7) + if(node_number) + text(xy(i,1)*1.1,xy(i,2)*1.1,num2str(i),'HorizontalAlignment','Center') + end + end + set(gca,'xtick',[],'ytick',[]) + xlim([-1.5 1.5]) + ylim([-1.5 1.5]) + box on + pbaspect([1 1 1]) +end \ No newline at end of file diff --git a/Plot_WU_Network.m b/Plot_WU_Network.m new file mode 100644 index 0000000..438308b --- /dev/null +++ b/Plot_WU_Network.m @@ -0,0 +1,105 @@ +function Plot_WU_Network(adjacency,xy,xy_colors,curved,link_color,node_size,node_number) +% Plot weighted and undirected network with loops +% +% Plot_WU_Network(adjacency,curved,xy,xy_colors,link_color,node_size,node_number) +% +% Pérez-Ortega Jesús - april 2018 +% Modified Jan 2019 +% Modified April 2019 + +n=length(adjacency); +switch (nargin) + case 6 + node_number=false; + case 5 + node_number=false; + node_size=10; + case 4 + node_number=false; + node_size=10; + link_color=[0.5 0.5 0.5]; + case 3 + node_number=false; + node_size=10; + link_color=[0.5 0.5 0.5]; + curved=0; + case 2 + node_number=false; + node_size=10; + link_color=[0.5 0.5 0.5]; + curved=0; + xy_colors = Read_Colors(n); + case 1 + node_number=false; + node_size=10; + link_color=[0.5 0.5 0.5]; + curved=0; + xy_colors = Read_Colors(n); + xy = []; +end + +if(size(xy_colors,1)==1) + xy_colors=repmat(xy_colors,n,1); +end + +lims = []; +if(nargin==1||isempty(xy)) + nodes=length(adjacency); + xy=Get_Circular_XY(nodes); + lims = [-1.5 1.5]; +% xy = Get_Force_XY(adjacency); +end + +C=length(adjacency); + +% Plot edges +if(sum(adjacency(:))) + line_widths = zeros(C); + line_widths(adjacency>0) = Scale(adjacency(adjacency>0),0.5,5); + hold on + for a=1:C + for b=a:C + if (adjacency(a,b)) + length_arrow=0; + Plot_Edge(xy(a,:),xy(b,:),curved,length_arrow,line_widths(a,b),... + link_color); + end + end + end +end + +links=sum(adjacency); +if(length(node_size)==1) + if(sum(links) && node_size) + sizes_in=Scale(sum(adjacency),node_size,node_size+40); + else + sizes_in=ones(1,C)*30; + end +else + sizes_in=Scale(node_size,10,50); +end + +% Plot nodes +for i=1:C + if(links(i)) + plot(xy(i,1),xy(i,2),'.k','MarkerSize',sizes_in(i)+5) + plot(xy(i,1),xy(i,2),'.','color',xy_colors(i,:),'MarkerSize',sizes_in(i)) + if(node_number) + %text(xy(i,1)*1.1,xy(i,2)*1.1,num2str(i),'HorizontalAlignment','Center') + text(xy(i,1),xy(i,2),num2str(i),'HorizontalAlignment','left') + end + else + plot(xy(i,1),xy(i,2),'.','color',mean([0.5 0.5 0.5; xy_colors(i,:)]),... + 'MarkerSize',sizes_in(i)*.2) + end +end +set(gca,'xtick',[],'ytick',[],'xcolor',[1 1 1],'ycolor',[1 1 1]) +if(lims) + xlim(lims) + ylim(lims) +else + xlim([min(xy(:,1)) max(xy(:,1))]) + ylim([min(xy(:,2)) max(xy(:,2))]) +end + +pbaspect([1 1 1]) diff --git a/RV_coefficient.m b/RV_coefficient.m new file mode 100644 index 0000000..6928d62 --- /dev/null +++ b/RV_coefficient.m @@ -0,0 +1,12 @@ +% RV coefficient +% +% Compare the correlation between matrices of differents column size. +% tr(XX'YY')/sqrt(tr((XX')^2)*tr((YY')^2)) +% +% Pérez-Ortega Jesús march-2018 + +function RV = RV_coefficient(X,Y) + XX=X*X'; + YY=Y*Y'; + RV=trace(XX*YY)/sqrt(trace(abs(XX^2))*trace(abs(YY^2))); +end \ No newline at end of file diff --git a/Read_Colors.m b/Read_Colors.m new file mode 100644 index 0000000..8c1fd71 --- /dev/null +++ b/Read_Colors.m @@ -0,0 +1,45 @@ +% Return colors +% +% Pérez-Ortega Jesús - May 2018 + +function colors = Read_Colors(num_colors) + 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 + +end \ No newline at end of file diff --git a/Remove_Oscillations.m b/Remove_Oscillations.m new file mode 100644 index 0000000..d1475fe --- /dev/null +++ b/Remove_Oscillations.m @@ -0,0 +1,15 @@ +% Remove smoothed signal from original +% +% Remove the smoothed loess signal from the original +% +% By Jes?s P?rez-Ortega april-2018 + +function [smoothed,removed]= Remove_Oscillations(coactivity,percentage) + if(nargin==1) + %percentage=0.2; + percentage=0.01; % 15 min at 1 kH --> 9 s + end + + removed=smooth(coactivity,percentage,'loess'); % quadratic fit + smoothed=coactivity-removed; +end \ No newline at end of file diff --git a/Reshape_Raster.m b/Reshape_Raster.m new file mode 100644 index 0000000..9674454 --- /dev/null +++ b/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/Save_Figure.m b/Save_Figure.m new file mode 100644 index 0000000..c9f19e5 --- /dev/null +++ b/Save_Figure.m @@ -0,0 +1,43 @@ +% Save current figure with default style +% +% Pérez-Ortega Jesús - June 2018 +% modified Feb 2019 +function Save_Figure(name,format) + + if (nargin==1) + s.Format = 'png'; + else + s.Format = format; + end + s.Version='1'; + s.Preview='none'; + s.Width='auto'; + s.Height='auto'; + s.Units='inches'; + s.Color='rgb'; + s.Background='w'; + s.FixedFontSize='14'; + s.ScaledFontSize='auto'; + s.FontMode='fixed'; + s.FontSizeMin='8'; + s.FixedLineWidth='1'; + s.ScaledLineWidth='auto'; + s.LineMode='none'; + s.LineWidthMin='0.5'; + s.FontName='auto'; + s.FontWeight='auto'; + s.FontAngle='auto'; + s.FontEncoding='latin1'; + s.PSLevel='3'; + s.Renderer='auto'; + s.Resolution='auto'; + s.LineStyleMap='none'; + s.ApplyStyle='0'; + s.Bounds='loose'; + s.LockAxes='on'; + s.LockAxesTicks='off'; + s.ShowUI='on'; + s.SeparateText='off'; + + hgexport(gcf,[name '.' s.Format ],s); +end \ No newline at end of file diff --git a/Save_Raster_By_Windows.m b/Save_Raster_By_Windows.m new file mode 100644 index 0000000..968399d --- /dev/null +++ b/Save_Raster_By_Windows.m @@ -0,0 +1,20 @@ +% Save figures second by second from raster +% +% Pérez-Ortega Jesús - May 2018 + +function Save_Raster_By_Windows(name,samples_per_second,window_sec,final_sec) + + Hold_Axes(['CoactiveAxes' name]); + set(gca,'xtick',0:window_sec:final_sec) + set(gca,'xticklabel',0:window_sec:final_sec) + xlabel('time (s)') + for i=1:window_sec:final_sec + Hold_Axes(['RasterAxes' name]); + xlim([i-1 i+window_sec-1]*samples_per_second) + Hold_Axes(['CoactiveAxes' name]); + xlim([i-1 i+window_sec-1]) + + % Configure and save image + Save_Figure([name '_' num2str(i-1,'%.2f') '-' num2str(i+window_sec-1,'%.2f') 'sec']); + end +end diff --git a/Save_Raster_Min_By_Min.m b/Save_Raster_Min_By_Min.m new file mode 100644 index 0000000..09380e8 --- /dev/null +++ b/Save_Raster_Min_By_Min.m @@ -0,0 +1,16 @@ +% Save figures minute by minute from raster +% +% Pérez-Ortega Jesús - March 2018 + +function Save_Raster_Min_By_Min(data_name,samples_per_minute,initial_min,final_min) + initial_min=initial_min+1; + for i=initial_min:final_min + Hold_Axes(['RasterAxes' data_name]); + xlim([(i-1)*samples_per_minute+1 i*samples_per_minute]) + Hold_Axes(['CoactiveAxes' data_name]); + xlim([(i-1)*60 i*60]) + + % Configure and save image + Save_Figure([data_name '_' num2str(i-1,'%.2i') '-' num2str(i,'%.2i') 'min']); + end +end diff --git a/Scale.m b/Scale.m new file mode 100644 index 0000000..3cfab2b --- /dev/null +++ b/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/Set_Axes.m b/Set_Axes.m new file mode 100644 index 0000000..e23082b --- /dev/null +++ b/Set_Axes.m @@ -0,0 +1,5 @@ +%% Set Axes +function Set_Axes(AxesName,Position) + axes('outerposition',Position) + set(gca,'Tag',AxesName) +end \ No newline at end of file diff --git a/Set_Colormap_Blue_White_Red.m b/Set_Colormap_Blue_White_Red.m new file mode 100644 index 0000000..81b8da5 --- /dev/null +++ b/Set_Colormap_Blue_White_Red.m @@ -0,0 +1,11 @@ +% Set colormap blue-white-red +% +% By Jesús Pérez-Ortega jan-2018 + +function Set_Colormap_Blue_White_Red() + bluemap=colormap(gray(32))+repmat([0 0 1],32,1); + bluemap(bluemap>1)=1; + redmap=colormap(flipud(gray(32)))+repmat([1 0 0],32,1); + redmap(redmap>1)=1; + colormap(gca,[bluemap;redmap]) +end \ No newline at end of file diff --git a/Set_Figure.m b/Set_Figure.m new file mode 100644 index 0000000..5c82403 --- /dev/null +++ b/Set_Figure.m @@ -0,0 +1,23 @@ +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 + +h = findobj('name',title_name); +if isempty(h) + h = figure('name',title_name,'numbertitle','off','position',position); +else + figure(h); +end +clf + +if nargout + fig = h; +end diff --git a/Set_Label_Time.m b/Set_Label_Time.m new file mode 100644 index 0000000..3eb76d5 --- /dev/null +++ b/Set_Label_Time.m @@ -0,0 +1,29 @@ +function Set_Label_Time(samples,sample_frequency) +% Set convenient label of time +% +% Set_Label_Time(samples,sample_frequency) +% +% Jesus Perez-Ortega April-19 + +% less than 30 seconds +if(samples/sample_frequency<10) + step = 1; + xlabel('Time (s)'); +elseif(samples/sample_frequency<30) + step = 3; + xlabel('Time (s)'); +elseif(samples/sample_frequency/60<3) + step = 10; + xlabel('Time (s)'); +elseif(samples/sample_frequency/60<60) + step = 60; + xlabel('Time (min)'); +else + step = 60*60; + xlabel('Time (h)'); +end + +set(gca,'box','off','xtick',0:step*sample_frequency:samples,... + 'xticklabel',0:samples/sample_frequency/step) + + diff --git a/Shuffle_Test.m b/Shuffle_Test.m new file mode 100644 index 0000000..5b1c467 --- /dev/null +++ b/Shuffle_Test.m @@ -0,0 +1,113 @@ +% Test for significant number of coactive neurons +% +% P?rez Ortega Jes?s E. - March 2018 +% Modified April 2018 + +function Th = Shuffle_Test(raster,smooth_window,n,shuffle_method,alpha,... + integer_interval,remove_oscillations,z_score,only_mean,zscore_window,plot_shuffled_figures) + + if(nargin==7) + plot_shuffled_figures=false; + elseif(nargin==6) + plot_shuffled_figures=false; + remove_oscillations=false; + elseif(nargin==5) + plot_shuffled_figures=false; + remove_oscillations=false; + integer_interval=true; + elseif(nargin==4) + plot_shuffled_figures=false; + remove_oscillations=false; + integer_interval=true; + alpha=0.05; + elseif(nargin==3) + plot_shuffled_figures=false; + remove_oscillations=false; + integer_interval=true; + alpha=0.05; + shuffle_method='time_shift'; + end + + [c,f]=size(raster); + smooth_coactivity = Get_And_Filter_Coactivity(raster,smooth_window); + + % Remove Oscillations + if(remove_oscillations) + smooth_coactivity = Remove_Oscillations(smooth_coactivity); + end + + % Z-score + if(z_score) + smooth_coactivity=Z_Score_Coactivity(smooth_coactivity,zscore_window,only_mean); + end + + if(integer_interval) + interval=1:c; + PS=zeros(n,c); + else + interval=0:.1:max(smooth_coactivity); + PS=zeros(n,length(interval)); + end + + % Plot only if necessary + if(plot_shuffled_figures) + 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(smooth_coactivity) + title('Observed') + y_limits=get(gca,'ylim'); + end + + % SHUFFLED data + for i=1:n + R_shuffled = shuffle(raster,shuffle_method); % Get raster shuffled + smooth_co = Get_And_Filter_Coactivity(R_shuffled,smooth_window); % Get coactivity + % Remove Oscillations + if(remove_oscillations) + smooth_co = Remove_Oscillations(smooth_co); + end + + % Z-score + if(z_score) + smooth_co=Z_Score_Coactivity(smooth_co,0,only_mean); + end + + HS = histc(smooth_co,interval); % Histogram of coactive cells + %PS(i,:) = cumsum(HS)/f; % Synchrony probability + PS(i,:) = cumsum(HS)/max(cumsum(HS)); % Synchrony probability + + % Plot only if necessary + if(plot_shuffled_figures && i<5) + Hold_Figure('Raster | real VS shuffled'); + subplot(5,1,i+1) + imagesc(R_shuffled) + title(['Shuffled ' num2str(i)]) + Hold_Figure('Coactivity | real VS shuffled') + subplot(5,1,i+1) + plot(smooth_co) + title(['Shuffled ' num2str(i)]) + ylim(y_limits) + end + end + PSmean=mean(PS); % Mean of Probability + + % Set threshold from Probability distribution < alpha + th_idx=find(PSmean>=(1-alpha),1,'first'); + Th=interval(th_idx); + + % Plot distribution shuffled + Set_Figure('Shuffled test',[0 0 600 300]); + plot(interval,PSmean,'o-b'); hold on + title({['Cumulative distribution of raster shuffled (' num2str(n) ' iterations; method '... + strrep(shuffle_method,'_','-') ')'];['th = ' num2str(Th) ' at alpha=' num2str(alpha)]}) + ylabel('P(x)') + xlabel('x') + if(~isnan(Th)) + xlim([0 Th*3]) + end +end \ No newline at end of file diff --git a/Shuffle_Test_For_Every_Link.m b/Shuffle_Test_For_Every_Link.m new file mode 100644 index 0000000..2765be1 --- /dev/null +++ b/Shuffle_Test_For_Every_Link.m @@ -0,0 +1,46 @@ +% Weight threshold +% +% Montecarlo for determining weight threshold in weighted matrix. +% +% Threshold is determined when the edges reach the 1-alpha percentage in +% its randomized versions (N times). +% +% Adjacency matrix is determined by coactivity, weight of connection means +% the times which are two nodes coactive. +% +% w_ths = Shuffle_Test_For_Every_Link(raster,shuffle_method,connectivity_method,iterations,alpha) +% +% Jesus E. Perez-Ortega - Sep 2018 +% modified june 2018 +% modified nov 2018 + +function w_ths = Shuffle_Test_For_Every_Link(raster,shuffle_method,connectivity_method,iterations,alpha) + + % Data size + n = size(raster,1); + adjacency_shuffle=zeros(iterations,n,n); + + for i=1:iterations + % Shuffled version + raster_shuffled = shuffle(raster,shuffle_method); + + % Adjacent from shuffled version + adjacency_shuffle(i,:,:)=Get_Adjacency_From_Raster(raster_shuffled,connectivity_method); + end + + % Cumulative distribution of individual link + w_ths = zeros(n); + for j=1:n-1 + for k=j+1:n + try + y=hist(adjacency_shuffle(:,j,k),0:max(adjacency_shuffle(:,j,k))); + cum=cumsum(y)/max(cumsum(y)); + th=find(cum>(1-alpha),1,'first')-1; + w_ths(j,k)=th; + w_ths(k,j)=w_ths(j,k); + catch + disp(th) + end + end + end +end \ No newline at end of file diff --git a/Shuffle_Test_For_Weighted_Matrix.m b/Shuffle_Test_For_Weighted_Matrix.m new file mode 100644 index 0000000..5807894 --- /dev/null +++ b/Shuffle_Test_For_Weighted_Matrix.m @@ -0,0 +1,67 @@ +% Weight threshold +% +% Montecarlo for determining weight threshold in weighted matrix. +% +% Threshold is determined when the edges reach the 1-alpha percentage in +% its randomized versions (N times). +% +% Adjacency matrix is determined by coactivity, weight of connection means +% the times which are two nodes coactive. +% +% W_Th = WeightTh(Data, N, alpha) +% +% ..:: by Jes?s E. P?rez-Ortega ::.. Jun-2013 +% modified june 2018 + +function w_threshold = Shuffle_Test_For_Weighted_Matrix(raster,shuffle_method,connectivity_method,iterations,alpha) + + % Data size + samples=size(raster,2); + + for i=1:iterations + % Shuffled version + raster_shuffled = shuffle(raster,shuffle_method); % Get raster shuffled + + % Adjacent from shuffled version + adjacency_shuffle=Get_Adjacency_From_Raster(raster_shuffled,connectivity_method); + + % Number of edges at each threshold + edges_th_shuffled=zeros(1,100); + if(strcmp(connectivity_method,'jaccard')) + interval=0:0.001:1; + else + interval=0:samples; + end + j=1; + for th=interval + edges_th_shuffled(j)=sum(sum(adjacency_shuffle>th)); + if(~edges_th_shuffled(j)) + break; + end + j=j+1; + end + + % Coactivation probability + p(i,1:length(edges_th_shuffled)) = cumsum(edges_th_shuffled)/max(cumsum(edges_th_shuffled)); + end + p(p==0)=1; + p_mean=mean(p); + interval=interval(1:length(p_mean)); + % Weight threshold + idx=find(p_mean>(1-alpha),1,'first'); + w_threshold=interval(idx); + + % Plot distribution shuffled + Set_Figure('Shuffled weigthed matrix test',[0 0 600 300]); + %plot(interval/samples,p_mean,'o-b'); + plot(interval,p_mean,'o-b'); + title({['Cumulative distribution of raster shuffled (' num2str(iterations) ' iterations; method '... + strrep(shuffle_method,'_','-') ')'];['th = ' num2str(w_threshold) ' at alpha=' num2str(alpha)]}) + ylabel('P(x)') + xlabel('x') + ylim([0 1.05]) + if(~isnan(w_threshold) && w_threshold>0) + xlim([0 (w_threshold*2)]) + end + %w_threshold=w_threshold/100; +end \ No newline at end of file diff --git a/Sort_Activation_Position.m b/Sort_Activation_Position.m new file mode 100644 index 0000000..4eab137 --- /dev/null +++ b/Sort_Activation_Position.m @@ -0,0 +1,63 @@ +% Sort activation position +% +% By Jesús Pérez-Ortega, Jan 2019 + +function [sequences_sorted,position,probabilities] = Sort_Activation_Position(sequences) + [trials,n] = size(sequences); + + % without "neuron 0" + % get the probability of activation + for i = 1:n + p = sum(sequences==i)/trials; + probabilities(i,:) = p; + end + + % Set the activation position + % C + [~,position] = Sort_Raster(probabilities); + + % sort the sequence + sequences_sorted = zeros(trials,n); + for i = 1:n + sequences_sorted(sequences==position(i))=i; + end + + % sort preference + probabilities = probabilities(position,:); + +% % set activation position (other methods) +% % A +% remaining = 1:n; +% for i = 1:n +% [~,id] = max(preferences(remaining,i)); +% position(i) = remaining(id); +% remaining = setdiff(remaining,remaining(id)); +% end +% +% % B +% for i = 1:n +% [~,position(i)] = max(preferences(i,:)); +% end +% [~,position] = sort(position); +end + +% % with "neuron 0" probability +% % get the probability of activation +% for i = 0:n +% p = sum(sequences==i)/trials; +% probabilities(i+1,:) = p; +% end +% +% % Set the activation position +% [~,position] = Sort_Raster(probabilities(2:end,:)); +% position = [1 position+1]; +% +% % sort the sequence +% sequences_sorted = zeros(trials,n); +% for i = 0:n +% sequences_sorted(sequences==position(i+1)-1)=i; +% end +% +% % sort preference +% probabilities = probabilities(position,:); +% position = position(2:end)-1; \ No newline at end of file diff --git a/Sort_Raster.m b/Sort_Raster.m new file mode 100644 index 0000000..d13418b --- /dev/null +++ b/Sort_Raster.m @@ -0,0 +1,21 @@ +function [sorted,sorted_id] = Sort_Raster(raster,direction) +% Sort raster +% +% [sorted,sorted_id] = Sort_Raster(raster,direction) +% +% Jesús Pérez-Ortega - Dic 2018 +% Modified Jan 2019 + + if(nargin==1) + direction = 'descend'; + end + + sorted = raster; + sorted_id = 1:size(raster,1); + n_seq = size(sorted,2); + for i = n_seq:-1:1 + [~,id] = sort(sorted(:,i),direction); + sorted = sorted(id,:); + sorted_id = sorted_id(id); + end +end \ No newline at end of file diff --git a/Z_Score_Coactivity.m b/Z_Score_Coactivity.m new file mode 100644 index 0000000..bfb1a40 --- /dev/null +++ b/Z_Score_Coactivity.m @@ -0,0 +1,51 @@ +% Get Z-score of coactivity with an specific time window +% +% +% By Jes?s P?rez-Ortega jan-2018 +% modified april-2018 +% modified august-2018 + +function z_score_coactivity = Z_Score_Coactivity(coactivity,zscore_window,only_mean) + if(nargin<2) + zscore_window=length(coactivity); + only_mean=false; + end + + if(zscore_window==0) + zscore_window=length(coactivity); + end + + z_score_coactivity=zeros(size(coactivity)); + F=length(coactivity); + if(F>zscore_window) + n_final=round(F/zscore_window); + for i=1:n_final + inicio=zscore_window*(i-1)+1; + fin=inicio+zscore_window-1; + if(fin>F) + fin=F; + end + if(only_mean) + z_score_coactivity(inicio:fin)=coactivity(inicio:fin)-mean(coactivity(inicio:fin)); + else + z_score_coactivity(inicio:fin)=(coactivity(inicio:fin)-mean(coactivity(inicio:fin)))/std(coactivity(inicio:fin)); + end + end + n_Smooth=length(z_score_coactivity); + if(fin=2 %degree must be at least 2 + S=G(V,V); + C(u)=sum(S(:))/(k^2-k); + end +end \ No newline at end of file diff --git a/density_und.m b/density_und.m new file mode 100644 index 0000000..87debf4 --- /dev/null +++ b/density_und.m @@ -0,0 +1,28 @@ +function [kden,N,K] = density_und(CIJ) +% DENSITY_UND Density +% +% kden = density_und(CIJ); +% [kden,N,K] = density_und(CIJ); +% +% Density is the fraction of present connections to possible connections. +% +% Input: CIJ, undirected (weighted/binary) connection matrix +% +% Output: kden, density +% N, number of vertices +% K, number of edges +% +% Notes: Assumes CIJ is undirected and has no self-connections. +% Weight information is discarded. +% +% +% Olaf Sporns, Indiana University, 2002/2007/2008 + + +% Modification history: +% 2009-10: K fixed to sum over one half of CIJ [Tony Herdman, SFU] + +N = size(CIJ,1); +K = nnz(triu(CIJ)); +kden = K/((N^2-N)/2); + diff --git a/dist_corr.m b/dist_corr.m new file mode 100644 index 0000000..9774871 --- /dev/null +++ b/dist_corr.m @@ -0,0 +1,6 @@ +function d = dist_corr(x,y) + A=x-mean(x); + B=y-mean(y); + + d = 1 - sum(A.*B)/sqrt(sum(A.*A).*sum(B.*B)); +end \ No newline at end of file diff --git a/distance_bin.m b/distance_bin.m new file mode 100644 index 0000000..5210b78 --- /dev/null +++ b/distance_bin.m @@ -0,0 +1,45 @@ +function D=distance_bin(A) +%DISTANCE_BIN Distance matrix +% +% D = distance_bin(A); +% +% The distance matrix contains lengths of shortest paths between all +% pairs of nodes. An entry (u,v) represents the length of shortest path +% from node u to node v. The average shortest path length is the +% characteristic path length of the network. +% +% Input: A, binary directed/undirected connection matrix +% +% Output: D, distance matrix +% +% Notes: +% Lengths between disconnected nodes are set to Inf. +% Lengths on the main diagonal are set to 0. +% +% Algorithm: Algebraic shortest paths. +% +% +% Mika Rubinov, U Cambridge +% Jonathan Clayden, UCL +% 2007-2013 + +% Modification history: +% 2007: Original (MR) +% 2013: Bug fix, enforce zero distance for self-connections (JC) + +A=double(A~=0); %binarize and convert to double format + +l=1; %path length +Lpath=A; %matrix of paths l +D=A; %distance matrix + +Idx=true; +while any(Idx(:)) + l=l+1; + Lpath=Lpath*A; + Idx=(Lpath~=0)&(D==0); + D(Idx)=l; +end + +D(~D)=inf; %assign inf to disconnected nodes +D(1:length(A)+1:end)=0; %clear diagonal \ No newline at end of file diff --git a/distcorr.m b/distcorr.m new file mode 100644 index 0000000..ff89017 --- /dev/null +++ b/distcorr.m @@ -0,0 +1,35 @@ +function dcor = distcorr(x,y) + % This function calculates the distance correlation between x and y. + % Reference: http://en.wikipedia.org/wiki/Distance_correlation + % Date: 18 Jan, 2013 + % Author: Shen Liu (shen.liu@hotmail.com.au) + % Check if the sizes of the inputs match + if size(x,1) ~= size(y,1) + error('Inputs must have the same number of rows') + end + % Delete rows containing unobserved values + N = any([isnan(x) isnan(y)],2); + x(N,:) = []; + y(N,:) = []; + % Calculate doubly centered distance matrices for x and y + a = pdist2(x,x); + mcol = mean(a); + mrow = mean(a,2); + ajbar = ones(size(mrow))*mcol; + akbar = mrow*ones(size(mcol)); + abar = mean(mean(a))*ones(size(a)); + A = a - ajbar - akbar + abar; + b = pdist2(y,y); + mcol = mean(b); + mrow = mean(b,2); + bjbar = ones(size(mrow))*mcol; + bkbar = mrow*ones(size(mcol)); + bbar = mean(mean(b))*ones(size(b)); + B = b - bjbar - bkbar + bbar; + % Calculate squared sample distance covariance and variances + dcov = sum(sum(A.*B))/(size(mrow,1)^2); + dvarx = sum(sum(A.*A))/(size(mrow,1)^2); + dvary = sum(sum(B.*B))/(size(mrow,1)^2); + % Calculate the distance correlation + dcor = sqrt(dcov/sqrt(dvarx*dvary)); +end \ No newline at end of file diff --git a/makerandCIJ_und.m b/makerandCIJ_und.m new file mode 100644 index 0000000..e03dbf1 --- /dev/null +++ b/makerandCIJ_und.m @@ -0,0 +1,25 @@ +function [CIJ] = makerandCIJ_und(N,K) +%MAKERANDCIJ_UND Synthetic directed random network +% +% CIJ = makerandCIJ_und(N,K); +% +% This function generates an undirected random network +% +% Inputs: N, number of vertices +% K, number of edges +% +% Output: CIJ, undirected random connection matrix +% +% Note: no connections are placed on the main diagonal. +% +% +% Olaf Sporns, Indiana University, 2007/2008 + +ind = triu(~eye(N)); +i = find(ind); +rp = randperm(length(i)); +irp = i(rp); + +CIJ = zeros(N); +CIJ(irp(1:K)) = 1; +CIJ = CIJ+CIJ'; % symmetrize diff --git a/rich_club_bu.m b/rich_club_bu.m new file mode 100644 index 0000000..e7e7651 --- /dev/null +++ b/rich_club_bu.m @@ -0,0 +1,49 @@ +function [R,Nk,Ek] = rich_club_bu(CIJ,varargin) +%RICH_CLUB_BU Rich club coefficients (binary undirected graph) +% +% R = rich_club_bu(CIJ) +% [R,Nk,Ek] = rich_club_bu(CIJ,klevel) +% +% The rich club coefficient, R, at level k is the fraction of edges that +% connect nodes of degree k or higher out of the maximum number of edges +% that such nodes might share. +% +% Input: CIJ, connection matrix, binary and undirected +% klevel, optional input argument. klevel sets the +% maximum level at which the rich club +% coefficient will be calculated. If klevel is +% not included the the maximum level will be +% set to the maximum degree of CIJ. +% +% Output: R, vector of rich-club coefficients for levels +% 1 to klevel. +% Nk, number of nodes with degree>k +% Ek, number of edges remaining in subgraph with +% degree>k +% +% Reference: Colizza et al. (2006) Nat. Phys. 2:110. +% +% Martijn van den Heuvel, University Medical Center Utrecht, 2011 + +Degree = sum(CIJ); %compute degree of each node + +if nargin == 1 + klevel = max(Degree); +elseif nargin == 2 + klevel = varargin{1}; +elseif nargin > 2 + error('number of inputs incorrect. Should be [CIJ], or [CIJ, klevel]') +end + +R = zeros(1,klevel); +Nk = zeros(1,klevel); +Ek = zeros(1,klevel); +for k = 1:klevel + SmallNodes=find(Degree<=k); %get 'small nodes' with degree <=k + subCIJ=CIJ; %extract subnetwork of nodes >k by removing nodes <=k of CIJ + subCIJ(SmallNodes,:)=[]; %remove rows + subCIJ(:,SmallNodes)=[]; %remove columns + Nk(k)=size(subCIJ,2); %number of nodes with degree >k + Ek(k)=sum(subCIJ(:)); %total number of connections in subgraph + R(k)=Ek(k)/(Nk(k)*(Nk(k)-1)); %unweighted rich-club coefficient +end \ No newline at end of file diff --git a/shuffle.m b/shuffle.m new file mode 100644 index 0000000..ff12568 --- /dev/null +++ b/shuffle.m @@ -0,0 +1,116 @@ +function [shuffled,rand_id] = shuffle(x,method) + +%shuffle Shuffles raster data using various different methods +% Shuffles spike data (0 or 1) using three differnt methods +% assumes rows are individual cells and columns are time frames +% +% 'frames' - shuffles the frames in time, maintains activity pattern of +% each frame +% +% 'time' - shuffles the activity of each individual cell in time +% each cell maintains its total level of activity +% +% 'time_shift' - shifts the time trace of each cell by a random amount +% each cell maintains its pattern of activity +% +% Methods fom synfire chains paper +% +% 'isi' - Inter-Spike Interval shuffling within cells +% each cell maintains in level of activity +% +% 'cell' - shuffles the activity at a given time between cells +% each frame maintains the number of active cells +% +% 'exchange' - exchange pairs of spikes across cells +% slow, but each cell and frame maintains level of activity +% +% jzaremba 01/2012 +% +% modified Perez-Ortega Jesus - Aug 2018 +% modified - Feb 2019 + +if nargin < 2 + method = 'frames'; + warning('Method frames was selected.') +end + +if ~any(strcmp(method, {'frames','time','time_shift','isi','cell','exchange'})) + method = 'frames'; + warning('Method frames was selected.') +end + +rand_id =[]; +shuffled=x; + +switch method + case 'frames' + n = size(x,2); + randp = randperm(n); + shuffled = sortrows([randp;x]')'; + shuffled = shuffled(2:end,:); + + case 'time' + n = size(x,2); + for i = 1:size(x,1) + randp = randperm(n); + temp = sortrows([randp; x(i,:)]')'; + shuffled(i,:) = temp(2,:); + end + + case 'time_shift' + n = size(x,2); + for i = 1:size(x,1) + randp = randi(n); + shuffled(i,:) = [ x(i,n-randp+1:n) x(i,1:n-randp) ]; + rand_id(i) = randp; + end + + case 'isi' + n = size(x,2); + shuffled = zeros(size(x,1),n); + + for i = 1:size(x,1) + % Pull out indices of spikes, get ISIs, buffer at start and end + isi = diff(find([1 x(i,:) 1])); + isi = isi(randperm(length(isi))); % Randomize ISIs + + temp = cumsum(isi); + temp = temp(1:end-1); % Removes the end spikes that were added + % Put the spikes back + shuffled(i,temp) = true; + end + + + case 'cell' + [n,len] = size(x); + for i = 1:len + randp = randperm(n); + temp = sortrows([randp' x(:,i)]); + shuffled(:,i) = temp(:,2); + end + + case 'exchange' + n = sum(x(:)); + for i = 1:2*n + randp = randi(n,1,2); + [r,c] = find(shuffled); + + % Make sure that the swap will actually do something + while randp(1)==randp(2) || r(randp(1))==r(randp(2)) || c(randp(1))==c(randp(2)) ||... + shuffled(r(randp(2)),c(randp(1)))==true ||... + shuffled(r(randp(1)),c(randp(2)))==true + randp = randi(n,1,2); + end + + % Swap + shuffled(r(randp(2)),c(randp(1))) = true; + shuffled(r(randp(1)),c(randp(1))) = false; + + shuffled(r(randp(1)),c(randp(2))) = true; + shuffled(r(randp(2)),c(randp(2))) = false; + + end + +end + + \ No newline at end of file