Update

Change_Sequence_Numbers.m

@@ -0,0 +1,24 @@
+function newSequence = Change_Sequence_Numbers(sequence,new)
+% Replace the numbers of the sequence by the new numbers given
+%
+%       newSequence = Change_Sequence_Numbers(sequence,new)
+%
+% By Jesus Perez-Ortega, Jun 2020
+
+nSeq = length(unique(sequence));
+nNew = length(unique(new));
+
+if nSeq == nNew
+    if iscolumn(new)
+        new = new';
+    end
+    newSequence = zeros(size(sequence));
+    j = 1;
+    for i = new
+        newSequence(sequence==i) = j;
+        j = j+1;
+    end
+else
+    error(['number of elements in "new" variable are not the same number os elements'...
+    'in "sequence" variable'])
+end

Cluster_Test.m

@@ -0,0 +1,102 @@
+function [recommended,indices] = Cluster_Test(treeOrData,similarity,metric,clusteringMethod,groups,fig)
+% Clustering indexes
+% Get indexes for evaluating clustering from hierarchical cluster tree
+%
+%       [recommended,indices] = Cluster_Test(treeOrData,similarity,metric,clusteringMethod,groups,fig)
+%
+%       default: metric = 'contrast'; clusteringMethod = 'hierarchical';
+%                groups = 2:30; fig = []
+%
+% Inputs:
+%      treeOrData = hierarchical cluster tree, or data for k-means
+%      similarity = 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')
+%      clusteringMethod,groups
+%      numFig = number of the figure to plot
+%
+% Outputs:
+% indices = clustering indices of 'metric' from 2 to 10 groups
+% recommended = recommended number of clusters
+%
+% ..:: by Jesús E. Pérez-Ortega ::.. Jun-2012
+% modified March-2018
+% modified Apr 2020
+
+switch(nargin)
+    case 2
+        metric = 'contrast';
+        clusteringMethod = 'hierarchical';
+        groups = 2:30;
+        fig = [];
+
+    case 3
+        clusteringMethod = 'hierarchical';
+        groups = 2:30;
+        fig = [];
+    case 4 
+        groups = 2:30;
+        fig = [];
+    case 5 
+        fig = [];
+end
+
+dist = 1-similarity;
+j = 1;
+for i = groups
+    switch(clusteringMethod)
+        case 'hierarchical'
+            T = cluster(treeOrData,'maxclust',i);
+        case 'kmeans'
+            T = kmeans(treeOrData,i);
+    end
+    g = max(T);
+
+    switch(metric)
+        case 'dunn'
+            indices(j) = DunnIdx_JP(g,dist,T);
+        case 'connectivity'
+            indices(j) = CCidx2_JP(g,similarity,T);
+        case 'davies'
+            indices(j) = DaviesIdx_JP(g,similarity,T);
+        case 'contrast'
+            indices(j) = ContrastIdx_JP(g,similarity,T);
+    end
+    j = j+1;
+end
+
+% Old way
+% [gIdx,g]=max(indices);
+% maximum = groups(g);
+
+% New way
+[~,id] = find(diff(indices)>0,1,'first');
+if isempty(id) || id==length(groups)-1
+    % The indices are decreasing, so select the first
+    recommended = groups(1);
+    id = 1;
+else
+    % Find the first peak of the indices
+    indicesCopy = indices;
+    indicesCopy(1:id) = 0;
+    [~,id] = find(diff(indicesCopy)<0,1,'first');
+    if isempty(id)
+        % If there is no peak find the first sudden increase
+        [~,id] = find(diff(diff(indicesCopy))<0,1,'first');
+        id = id+1;
+    end
+    recommended = groups(id);
+end
+
+% Plot 
+if ~isempty(fig)
+    plot(groups,indices)
+    hold on
+    plot(recommended,indices(id),'*r')
+    hold off
+    title([metric '''s index (' num2str(recommended) ' groups recommended)'])
+    xlabel('number of groups')
+    ylabel('index value')
+end

Compare_Ensembles.m

@@ -0,0 +1,50 @@
+function ids = Compare_Ensembles(structureA,structureB,simTh)
+% Compare ensembles by its neural identity using a threshold of imilarity
+%
+%       ids = Compare_Ensembles(structureA,structureB,simTh)
+%
+% By Jesus Perez-Ortega, May 2020
+
+% Compute similarity between ensemble neurons
+simNeurons = 1-pdist2(double(structureA>0),double(structureB>0),'Jaccard')';
+[maxSim,maxSimID] = max(simNeurons,[],1);
+
+% Identify ensembles from A above similarity threshold
+sameA = find(maxSim>=simTh);
+
+% Identify similar ensembles from B
+sameB = maxSimID(sameA);
+nSameB = length(sameB);
+nSameUniqueB = length(unique(sameB));
+
+% Check for more than one similar
+if nSameB>nSameUniqueB
+    disp('Some ensembles from B are similar to more than one ensemble from A.')
+    idRemove = [];
+    for i = 1:nSameUniqueB
+        id = find(sameB==sameB(i));
+        if length(id)>1
+            % Detect more similar
+            [~,maxID] = max(maxSim(id));
+
+            % Keep more similar and remove the others
+            idKeep = id(maxID);
+            idRemove = [idRemove setdiff(id,idKeep)];
+        end
+    end
+    sameA(idRemove) = [];
+    sameB(idRemove) = [];
+    disp([num2str(length(idRemove)) ' ensembles were removed.'])
+end
+
+% Get number of ensembles
+nA = size(structureA,1);
+nB = size(structureB,1);
+
+% Add values to ids structure
+ids.SameA = sameA;
+ids.SameB = sameB;
+ids.DiffA = setdiff(1:nA,sameA);
+ids.DiffB = setdiff(1:nB,sameB);
+ids.nA = nA;
+ids.nB = nB;

Ensemble_Correlation.m

@@ -0,0 +1,27 @@
+function [correlation,cell_indices] = Ensemble_Correlation(raster,indices,sequence,group)
+% Compute correlation between neurons and single ensemble
+%
+%       [correlation,cell_indices] = Ensemble_Correlation(raster,indices,sequence,group)
+%
+% By Jesus Perez-Ortega, Apr 2020
+
+neurons = size(raster,1);
+
+% find indices of group peaks
+idGroup = find(sequence==group)';
+idIndices = zeros(size(indices));
+for i = idGroup
+    idIndices(indices==i) = 1;
+end
+%signal = signal.*idIndices;  % old
+signal = idIndices;           % correlation with binary
+
+% compute correlation
+for i=1:neurons
+    D(i) = pdist([signal';raster(i,:)],'correlation');
+end
+correlation = 1-D;
+correlation(isnan(correlation)) = 0;
+correlation(correlation<0) = 0;
+% sort
+[~, cell_indices] = sort(correlation,'descend');

Find_Peaks_Or_Valleys.m

@@ -1,16 +1,19 @@
-function [indices,widths,amplitudes,ini_fin] = Find_Peaks_Or_Valleys(data,threshold,join,...
+function [indices,widths,amplitudes,ini_fin_times] = 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,...
+%   [indices,widths,amplitudes,ini_fin_times] = Find_Peaks_Or_Valleys(data,threshold,join,...
 %    detect_peaks,minimum_width,fixed_width,ignore_ini_fin)
 %
+%           default: threshold = 0; join = true; detect_peaks = true;
+%                    minimum_width = 0; fixed_width = 0; ignore_ini_fin = false;
+%
 % Inputs
 % data = data as vector Fx1 (F = #frames)
 % threshold = threshold
 % join = set mode to get peaks (0 = each vector above threshold is a peak;
 %        1 = joining of adjacent vectors above the threshold is a peak)
-% 
+%
 % Outputs
 % indices = Fx1 vector containing the peak indices
 %
@@ -21,26 +24,37 @@
 % 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
+% by Jesus E. Perez-Ortega, Feb-2012
+% last modification July-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;
+switch nargin
+    case 6
+        ignore_ini_fin = false;
+    case 5      
+        ignore_ini_fin = false;
+        fixed_width = 0;
+    case 4
+        ignore_ini_fin = false;
+        fixed_width = 0;
+        minimum_width = 0;
+    case 3
+        ignore_ini_fin = false;
+        fixed_width = 0;
+        minimum_width = 0;
+        detect_peaks = true;
+    case 2
+        ignore_ini_fin = false;
+        fixed_width = 0;
+        minimum_width = 0;
+        detect_peaks = true;
+        join = true;
+    case 1
+        ignore_ini_fin = false;
+        fixed_width = 0;
+        minimum_width = 0;
+        detect_peaks = true;
+        join = true;
+        threshold = 0;
 end
 
 % 0. Correct signal data
@@ -59,12 +73,12 @@
         disp('No peaks found!')
         widths = [];
         amplitudes = [];
-        ini_fin = [];
+        ini_fin_times = [];
     else
         disp('No valleys found!')
         widths = [];
         amplitudes = [];
-        ini_fin = [];
+        ini_fin_times = [];
     end
     return
 end
@@ -90,7 +104,7 @@
 
     % Delete if ends above threshold
     last = F;
-    idx = flipud(idx);
+    idx = fliplr(idx);
     for i = idx
         if last==i
             if detect_peaks
@@ -250,11 +264,11 @@
 
 % Get peaks or valleys width
 widths = zeros(count,1);
-ini_fin = zeros(count,2);
+ini_fin_times = zeros(count,2);
 for i = 1:count
     id = find(indices==i);
-    ini_fin(i,1) = id(1);
-    ini_fin(i,2) = id(end);
+    ini_fin_times(i,1) = id(1);
+    ini_fin_times(i,2) = id(end);
     widths(i) = length(id);
 end
 

Get_Adjacency_From_Raster.m

@@ -1,27 +1,34 @@
+function adjacency = Get_Adjacency_From_Raster(raster,connectivity_method)
 % Get adjacency matrix from raster peaks
 %
+%       adjacency = Get_Adjacency_From_Raster(raster,connectivity_method)
+%
+%       default: connectivity_method = 'coactivity'
+%
 % Negative numbers and NaNs are set to 0s
 %
-% Pérez-Ortega Jesús - march 2018
+% Perez-Ortega Jesus - 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    
+if nargin==1
+    connectivity_method = 'coactivity';
+end
+
+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