Re-initialization

Author: hornauerp

.gitattributes

@@ -0,0 +1 @@
+*.npy filter=lfs diff=lfs merge=lfs -text

.gitignore

@@ -0,0 +1,3 @@
+*.asv
+*.mexa64
+*.nfs*

Classes/MEArecording.m

@@ -0,0 +1,128 @@
+classdef Unit < handle
+    properties
+        MEArecording
+        ReferenceWaveform
+        ReferenceElectrode
+        TemplateID %Added
+        SpikeTimes
+        WaveformFeatures
+        ActivityFeatures
+        RegularityFeatures
+        GraphFeatures
+        Catch22
+    end
+
+    properties(SetObservable = true)
+        ClusterID
+    end
+
+    properties(Dependent)
+        unitID
+        ACG
+        FullACG
+    end
+    
+    methods (Static)
+                
+        function feature_names = returnFeatureNames(feature_group)
+            switch feature_group
+                case "act"
+                    feature_names = ["FiringRate","MeanInterSpikeInterval","VarianceInterSpikeInterval","CVInterSpikeInterval","PartialAutocorrelation"];
+                case "reg"
+                    feature_names = ["RegularityFrequency","RegularityMagnitude","RegularityFit"];
+                case "c22"
+                    feature_names = string(GetAllFeatureNames());
+            end
+        end
+        
+    end
+
+    
+    methods
+        
+        function unit = Unit(mearec, waveform, reference_electrode, spike_times, waveform_features)
+            if nargin > 0
+                unit.MEArecording = mearec;
+                unit.ReferenceWaveform = waveform;
+                unit.ReferenceElectrode = reference_electrode;
+                unit.SpikeTimes = spike_times;
+                unit.WaveformFeatures = struct2table(waveform_features);
+                unit.inferActivityFeatures();
+            end
+        end
+
+        function unit_id = get.unitID(unit)
+            unit_id = find(unit.MEArecording.Units == unit);
+        end
+
+        function acg = get.ACG(unit)
+            acg = unit.MEArecording.Connectivity.CCG.CCGs(:,unit.unitID,unit.unitID);
+        end
+
+        function acg = get.FullACG(unit)
+            acg = unit.MEArecording.Connectivity.FullCCG.CCGs(:,unit.unitID,unit.unitID);
+        end
+        
+        function [eCCG,iCCG] = getCCGconnections(unit)
+            % unit.MEArecording.Connectivity.C
+            % e_idx = 
+        end
+        
+        function inferActivityFeatures(unit)
+            if length(unit.SpikeTimes) > 2
+                isi = diff(unit.SpikeTimes);
+                act_feat.FiringRate = length(unit.SpikeTimes)/unit.MEArecording.RecordingInfo.Duration;
+                act_feat.MeanInterSpikeInterval = mean(isi);
+                act_feat.VarianceInterSpikeInterval = var(isi);
+                act_feat.CVInterSpikeInterval = std(isi)/mean(isi);
+                pacf = parcorr(isi,1);
+                act_feat.PartialAutocorrelation = pacf(2);
+                unit.ActivityFeatures = struct2table(act_feat);
+                if unit.MEArecording.Parameters.Analyses.Regularity
+                    unit.RegularityFeatures = unit.getRegularity();
+                end
+                if unit.MEArecording.Parameters.Analyses.Catch22
+                    catch_22_table = unit.MEArecording.run_catch_22(unit.SpikeTimes);
+                    catch_22_table.Properties.VariableNames = "SC_" + string(catch_22_table.Properties.VariableNames);
+                    unit.Catch22 = catch_22_table;
+                    
+                end
+            end
+        end
+        
+        function regularity_table = getRegularity(unit)
+           binned_activity = histcounts(unit.SpikeTimes,'BinLimits',[0 unit.MEArecording.RecordingInfo.Duration],'BinWidth',unit.MEArecording.Parameters.Regularity.Binning);
+           norm_activity = binned_activity/max(binned_activity);
+           NFFT = length(norm_activity);
+           F = (0 : 1/NFFT : 1/2-1/NFFT)*(1/unit.MEArecording.Parameters.Regularity.Binning);
+           TEMP = fft(norm_activity,NFFT);
+           TEMP(1) = 0;
+           freq_domain = abs(TEMP(1:NFFT/2));
+           [mag,freq_idx] = max(freq_domain);
+           reg_feat.RegularityFrequency = F(freq_idx);
+           reg_feat.RegularityMagnitude = mag;
+           
+           norm_freq_domain = freq_domain/max(freq_domain);
+           l = [];
+           p = [];
+           for i = 1:unit.MEArecording.Parameters.Regularity.N_peaks
+               [pks,locs,~,~] = findpeaks(norm_freq_domain,'NPeaks',1,'SortStr','descend');
+               norm_freq_domain = norm_freq_domain(locs:end);
+               l = [l; locs]; %#ok<AGROW>
+               p = [p; pks]; %#ok<AGROW>
+               if isempty(pks) || length(norm_freq_domain) < 10
+                   break
+               end
+           end
+           log_p = log10(p)-min(log10(p)); % Make data positive to be able to fit
+           try
+               f = fit(cumsum(l),log_p,'exp1');
+               reg_feat.RegularityFit = f.b;
+           catch
+               reg_feat.RegularityFit = NaN;
+           end
+           regularity_table = struct2table(reg_feat);
+       end
+       
+    end
+end

Classes/RecordingGroup.m

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:76a47cffdac296e3d4fd0a287a71fd13a975e37a5332b8abc7d4d562e8243436
+size 4556112

Classes/Unit.m

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:408e6bf9cb58cfc34c7d1b9f837ec3b69aef367ed7bb814ec24e660220294684
+size 9112144

Data/Figure3/age_pred_feature_names.mat

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ec07ee8ccc8aba16eb04a767ee2c3244ef5b3c2f2a6fd3d3eabadf2af7bd0a16
+size 536968880

Data/Figure3/base_age_prediction.mat

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:cb98c3ebc7ec0c5e4bbcb946cc9d2f34f8619ef0197fc9ab8c3976192281a430
+size 7517860

Data/Figure3/comb_clustering.mat

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d0c7f5c4c1fa30090692a3cb7aee7ee7f7347eabd9388a803b113766ab3fed17
+size 15035640

Data/Figure3/feature_group_heatmap.mat

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c9bcda01ab7c92d089a1a6049011dbfccc32c640279614b654abc1d6f0543c80
+size 717631280

Data/Figure3/network_clustering.mat

@@ -0,0 +1,164 @@
+%CCG - Compute multiple cross- and auto-correlograms
+%
+%  USAGE
+%
+%    ccg_params.duration = 0.05,<options>)
+%
+%    times          times of all events
+%                   NOTE: spiketimes in SECONDS.
+%    groups         group IDs for each event in time list (should be
+%                   integers 1:nGroups)
+%
+%   alternate buzcode usage: 
+%   times           {Ncells} array of [Nspikes]  spiketimes for each cell.
+%   groups          []
+%
+%                   spikes = bz_GetSpikes
+%                   [ccg,t] = CCG(spikes.times,[])  
+%
+%    <options>      optional list of property-value pairs (see table below)
+%
+%    =========================================================================
+%     Properties    Values
+%    -------------------------------------------------------------------------
+%     'binSize'     bin size in s (default = 0.01)
+%     'duration'    duration in s of each xcorrelogram (default = 2)
+%     'norm'        normalization of the CCG, 'counts' or 'rate' (DL added 8/1/17)
+%                   'counts' gives raw event/spike count,
+%                   'rate' returns CCG in units of spks/second (default: counts)
+%    =========================================================================
+%
+%
+%  OUTPUT
+%   ccg     [t x ngroups x ngroups] matrix where ccg(t,i,j) is the
+%           number (or rate) of events of group j at time lag t with  
+%           respect to reference events from group i
+%   t       time lag vector (units: seonds)
+%
+%  SEE
+%
+%    See also ShortTimeCCG.
+
+% Copyright (C) 2012 by Michaël Zugaro
+%
+% This program is free software; you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation; either version 3 of the License, or
+% (at your option) any later version.
+
+function [ccg,t] = CCG(times,groups,varargin)
+
+% Default values
+duration = 2;
+binSize = 0.01;
+Fs = 1/20000;
+normtype = 'counts';
+
+% Option for spike times to be in {Ncells} array of spiketimes DL2017
+if iscell(times) && isempty(groups)
+    numcells = length(times);
+    for cc = 1:numcells
+        groups{cc}=cc.*ones(size(times{cc}));
+    end
+    times = cat(1,times{:}); groups = cat(1,groups{:});
+end
+
+%Sort
+[times,sortidx] = sort(times);
+groups = groups(sortidx);
+
+% Check parameters
+if nargin < 2,
+  error('Incorrect number of parameters (type ''help <a href="matlab:help CCG">CCG</a>'' for details).');
+end
+%if ~isdvector(times),
+%	error('Parameter ''times'' is not a real-valued vector (type ''help <a href="matlab:help CCG">CCG</a>'' for details).');
+%end
+if ~isdscalar(groups) && ~isdvector(groups),
+	error('Parameter ''groups'' is not a real-valued scalar or vector (type ''help <a href="matlab:help CCG">CCG</a>'' for details).');
+end
+if ~isdscalar(groups) && length(times) ~= length(groups),
+	error('Parameters ''times'' and ''groups'' have different lengths (type ''help <a href="matlab:help CCG">CCG</a>'' for details).');
+end
+groups = groups(:);
+times = times(:);
+
+% Parse parameter list
+for i = 1:2:length(varargin),
+	if ~ischar(varargin{i}),
+		error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help CCG">CCG</a>'' for details).']);
+	end
+	switch(lower(varargin{i})),
+		case 'binsize',
+			binSize = varargin{i+1};
+			%if ~isdscalar(binSize,'>0'),
+		%		error('Incorrect value for property ''binSize'' (type ''help <a href="matlab:help CCG">CCG</a>'' for details).');
+	%		end
+		case 'duration',
+			duration = varargin{i+1};
+			if ~isdscalar(duration,'>0'),
+				error('Incorrect value for property ''duration'' (type ''help <a href="matlab:help CCG">CCG</a>'' for details).');
+            end
+            
+            case 'Fs',
+			Fs = varargin{i+1};
+			if ~isdscalar(Fs,'>0'),
+				error('Incorrect value for property ''Fs'' (type ''help <a href="matlab:help CCG">CCG</a>'' for details).');
+            end
+        case 'norm'
+            normtype = varargin{i+1};
+     
+  end
+end
+
+
+
+% Number of groups, number of bins, etc.
+if length(groups) == 1,
+	groups = ones(length(times),1);
+	nGroups = 1;
+else
+	nGroups = max(unique(groups));
+end
+
+
+halfBins = round(duration/binSize/2);
+nBins = 2*halfBins+1;
+t = (-halfBins:halfBins)'*binSize;
+times = round(times/Fs);
+binSize_Fs = round(binSize/Fs);
+if length(times) <= 1,
+	% ---- MODIFIED BY EWS, 1/2/2014 ----
+    % *** Use unsigned integer format to save memory ***
+    ccg = uint16(zeros(nBins,nGroups,nGroups));
+    % -----------------------------------
+    return
+end
+
+% Compute CCGs
+nEvents = length(times);
+% 
+
+counts = double(CCGHeart(times,uint32(groups),binSize_Fs,uint32(halfBins)));
+% -----------------------------------
+% 
+% Reshape the results
+n = max(groups);
+counts = reshape(counts,[nBins n n]);
+
+
+if n < nGroups,
+	counts(nBins,nGroups,nGroups) = 0;
+end
+
+%Rate normalization: counts/numREFspikes/dt to put in units of spikes/s. DL
+switch normtype
+    case 'rate'
+        for gg = 1:nGroups
+            numREFspikes = sum(groups==gg);%number of reference events for group
+            counts(:,gg,:) = counts(:,gg,:)./numREFspikes./binSize;
+        end
+end
+        
+        
+ccg = flipud(counts);