Added matlab routines to deal with toefiles coming from klustakwik.

toefile_utils/generate_cell_groups.m

@@ -0,0 +1,49 @@
+% From population vectors, generate cell groups
+% Must activate Plex library before running
+clear
+close all
+
+load('st1215_cat_P01_S01_2ndPen_moremerged_20151102T090824_popvec.mat');
+
+%population_vectors(ncell, nwin, nstim, nbinstart);
+[ncell, nwin, nstim, n_bin_start] = size(population_vectors);
+
+%first compute average firing rate for cell:
+avg_fr = sum(sum(sum(population_vectors, 4), 3), 2)/(nwin*nstim*n_bin_start);
+
+%threshold population vectors with FR > 6*avgFR
+cell_thresh = repmat(6*avg_fr, [1, nwin, nstim, n_bin_start]);
+cell_groups = population_vectors > cell_thresh;
+
+% pull out activity during stimulus
+% windows_in_stim = (windows_starts >= 0 - bin_width & (windows_starts + bin_width) <= 1500 + bin_width)';
+% window_win = squeeze(windows_in_stim(:, 2));
+% cell_groups = cell_groups(:, window_win, :, :);
+
+% Now, build simplicial complex
+stream = api.Plex4.createExplicitSimplexStream();
+
+% Add vertices for each cell
+for i = 1:ncell
+    stream.addVertex(i);
+end
+
+%% Now go through each population vector and add simplicial complex
+for win = 1:nwin
+    for stim = 1:nstim
+        for binst = 1:n_bin_start
+            cell_group_this = squeeze(cell_groups(:, win, stim, binst));
+            cells_in_this_group = find(cell_group_this);
+            if(any(cells_in_this_group))
+                cells_in_this_group;
+                stream.addElement(cells_in_this_group);
+                stream.ensureAllFaces();
+            end
+        end
+    end
+end
+stream.finalizeStream();
+
+%% Compute homology
+persis = api.Plex4.getModularSimplicialAlgorithm(10, 7);
+persis.computeIntervals(stream)

toefile_utils/generate_population_vectors.m

@@ -0,0 +1,75 @@
+%cell groups
+% Generate population vectors as in Curto&Itskov 2008
+
+clear
+close all
+
+datafilename = 'st1215_cat_P01_S01_2ndPen_moremerged_20151102T090824';
+load(strcat(datafilename, '.mat'));
+
+%toedata = data_to_save.toedata;
+% parameters for analysis
+bin_width = 10; % time bin width for cell groups in ms
+n_bin_start = 5; % number of different bin start times.
+pre_stim_ms = 2000; % number of milliseconds pre-stim
+post_stim_ms = 2000; % number of milliseconds post-stim
+stim_ms = 1000; % stimuli duration in milliseconds
+
+fs = 31250;
+
+trial_duration = pre_stim_ms + post_stim_ms + stim_ms;
+nwin = floor(trial_duration / bin_width);
+
+starts_dt = floor(bin_width/n_bin_start); % equally spaced in each win.
+
+% get n_bin_start starting points for the windows within one bin time frame
+bin_starts = starts_dt.*[0:n_bin_start-1] - pre_stim_ms;
+
+windows_starts = zeros(n_bin_start, nwin);
+windows_ends = zeros(n_bin_start, nwin);
+
+windows_starts(:, 1) = bin_starts;
+windows_ends(:, 1) = bin_starts + bin_width;
+for i = 2:nwin
+    windows_starts(:, i) = windows_starts(:, i-1) + bin_width;
+    windows_ends(:, i) = windows_ends(:, i-1) + bin_width;
+end
+
+% build population vectors
+ncells = length(toedata);
+nstim = 41;
+population_vectors = zeros(ncells, nwin, nstim, n_bin_start);
+
+% For each cell, for each stim, for each trial, count the number of spikes
+% in each time bin to arrive at population vectors.
+for cell = 1:ncells
+    for stim = 1:nstim
+        spiketimes = toedata{cell,1}.stims{stim, 1}.toes;
+        ntrials = length(spiketimes);
+        for bin_start = 1:n_bin_start
+            spikecounts_thiscellstimbinst = zeros(1, 500);
+            for trial = 1:ntrials
+                spikes_this_trial = 1000*spiketimes{trial, 1}; %convert to ms
+                Nspikes = length(spikes_this_trial);
+                repspike = repmat(spikes_this_trial, [1, nwin]);
+                repbinst = repmat(windows_starts(bin_start, :), [Nspikes, 1]);
+                repbinend = repmat(windows_ends(bin_start, :), [Nspikes, 1]);
+                bool_spike = repspike >= repbinst & repspike < repbinend;
+                spikecounts_thistrial = sum(bool_spike, 1);
+                spikecounts_thiscellstimbinst = spikecounts_thiscellstimbinst + spikecounts_thistrial;
+            end
+            rate_this = spikecounts_thiscellstimbinst / bin_width;
+            population_vectors(cell, :, stim, bin_start) = rate_this;
+        end
+    end
+end
+
+% Save population vector data
+% popvec_data_to_save = rmfield(data_to_save, 'toedata');
+% popvec_data_to_save.popvecs = population_vectors;
+% popvec_data_to_save.binwidth = bin_width;
+% popvec_data_to_save.winstarts = windows_starts;
+popvecfilename = strcat(datafilename, '_popvec');
+save(popvecfilename, 'population_vectors', 'bin_width', 'windows_starts');
+
+

toefile_utils/kwik2mat.m

@@ -0,0 +1,177 @@
+% kwik to mat.  Get good units out of sorted kwik file with stimulus
+% information and store in one well formatted mat file
+% Brad Theilman August/September 2015
+
+clear
+close all
+
+kwikfile = './st1215/kwdfiles/pen2_latestpipeline/st1215_cat_P01_S01_2ndPen.kwik';
+
+% get information
+kwikinfo = h5info(kwikfile);
+
+% Get cluster and time information for each spike
+spike_clusters = h5read(kwikfile, '/channel_groups/0/spikes/clusters/main');
+spike_times = h5read(kwikfile, '/channel_groups/0/spikes/time_samples');
+
+% extract all the unique cluster ids:
+clusters = unique(spike_clusters);
+
+%% Build a boolean vector with 1 by a cluster id if it is a good cluster
+good_cluster_inds = zeros(1, length(clusters));
+
+for i = 1:length(clusters)
+
+    cluster_num = clusters(i);
+    cluster_attr_path = strcat('/channel_groups/0/clusters/main/', num2str(cluster_num));
+    cluster_group = h5readatt(kwikfile, cluster_attr_path, 'cluster_group');
+    disp(cluster_group);
+    %cluster_group 2 is Good
+    % 1 is MUA
+    if cluster_group == 1
+        good_cluster_inds(i) = 1;
+
+    end
+
+end
+% vector containing IDs of good clusters
+good_cluster_ids = clusters(logical(good_cluster_inds));
+
+%% Now, pull all spikes that belong to good clusters
+good_spikes = cell(length(good_cluster_ids), 2);
+for i = 1:length(good_cluster_ids)
+    cluster_num = good_cluster_ids(i);
+    %pull spikes from teh good cluster
+    this_cluster_spiketimes = spike_times(spike_clusters == cluster_num);
+    good_spikes{i, 1} = cluster_num;
+    good_spikes{i, 2} = this_cluster_spiketimes;
+end
+
+%% Now get stimulus information
+
+digmark_timesamples = h5read(kwikfile, '/event_types/DigMark/time_samples');
+digmark_codes = cell2mat(h5read(kwikfile, '/event_types/DigMark/codes'));
+stim_timesamples = h5read(kwikfile, '/event_types/Stimulus/time_samples');
+stim_codes = h5read(kwikfile, '/event_types/Stimulus/codes');
+stim_text = h5read(kwikfile, '/event_types/Stimulus/text');
+
+stim_start_times = double(digmark_timesamples(digmark_codes == '<'));
+stim_end_times = double(digmark_timesamples(digmark_codes == '>'));
+intertrial_start_times = double(digmark_timesamples(digmark_codes == '('));
+intertrial_end_times = double(digmark_timesamples(digmark_codes == ')'));
+
+stim_start_filename = cell2mat(stim_text);
+stim_start_filename = stim_start_filename(2:2:end, :);
+stim_start_filename = cellstr(stim_start_filename);
+
+% stim_start_times = h5read(kwikfile, '/event_types/60/time_samples');
+% stim_start_filename = h5read(kwikfile, '/event_types/60/stim_filename');
+% stim_end_times = h5read(kwikfile, '/event_types/62/time_samples');
+% intertrial_start_times = h5read(kwikfile, '/event_types/40/time_samples');
+% intertrial_end_times = h5read(kwikfile, '/event_types/41/time_samples');
+
+stim_files_unique = unique(stim_start_filename);
+nstims = length(stim_files_unique);
+
+
+%%
+% For each stimulus, get number of trials for that stimulus
+% start building final_mat
+final_data = struct;
+final_data.nstims = nstims;
+numtrials = zeros(1, nstims);
+stim_data = cell(nstims, 1);
+
+for i = 1:nstims
+
+    stim_entry = struct;
+    numtrials(i) = sum(strcmp(stim_files_unique(i), stim_start_filename));
+
+    stim_start_times_this_stim = stim_start_times(strcmp(stim_files_unique(i), stim_start_filename));
+    stim_end_times_this_stim = stim_end_times(strcmp(stim_files_unique(i), stim_start_filename));
+    stim_entry.name = stim_files_unique(i);
+    stim_entry.start_times = stim_start_times_this_stim;
+    stim_entry.end_times = stim_end_times_this_stim;
+    stim_entry.ntrials = sum(strcmp(stim_files_unique(i), stim_start_filename));
+    stim_data{i, 1} = stim_entry;
+
+end
+
+%% Go through each cell and build final data matrix
+
+fs = 31250.0; % TODO: get fs from data file itself 
+pre_stim_duration = 2; %in seconds
+post_stim_duration = 2; %seconds
+pre_stim_duration_samps = pre_stim_duration*fs;
+post_stim_duration_samps = post_stim_duration*fs;
+
+n_good_units = length(good_cluster_ids);
+toedata = cell(n_good_units, 1);
+
+
+for unit_num = 1:length(good_cluster_ids)
+
+    unit_entry = struct;
+    unit_entry.id = good_cluster_ids(unit_num);
+    stims = cell(nstims, 1);
+
+    % Loop through each stimulus
+    for stim_num = 1:nstims
+
+        stim_unit_entry = struct;
+
+        this_stim_data = stim_data{stim_num, 1};
+        stim_unit_entry.name = this_stim_data.name;
+        stim_unit_entry.ntrials = this_stim_data.ntrials;
+
+        %create trial times
+        trial_start_times = this_stim_data.start_times - pre_stim_duration_samps;
+        trial_end_times = this_stim_data.end_times + post_stim_duration_samps;
+        stim_unit_entry.stim_start_times = this_stim_data.start_times;
+        stim_unit_entry.stim_end_times = this_stim_data.end_times;
+        stim_unit_entry.trial_start_times = trial_start_times;
+        stim_unit_entry.trial_end_times = trial_end_times;
+
+        %Go through each trial to divy up spike times relative to stim
+        %start
+        toes = cell(this_stim_data.ntrials, 1);
+        for trialnum = 1:this_stim_data.ntrials
+            trial_start = trial_start_times(trialnum);
+            trial_end = trial_end_times(trialnum);
+
+            spiketimes_thisunit = double(good_spikes{unit_num, 2});
+            spiketimes_samps_thistrial = spiketimes_thisunit(spiketimes_thisunit >= trial_start & spiketimes_thisunit <= trial_end);
+            spiketimes_samps_thistrial_relstimonset = spiketimes_samps_thistrial - this_stim_data.start_times(trialnum);
+            spiketimes_secs_thistrial_relstimonset = double(spiketimes_samps_thistrial_relstimonset) / fs;
+
+            toes{trialnum, 1} = spiketimes_secs_thistrial_relstimonset;
+        end
+        stim_unit_entry.toes = toes;
+        stims{stim_num, 1} = stim_unit_entry;     
+    end
+
+    unit_entry.stims = stims;
+    toedata{unit_num, 1} = unit_entry;
+
+end
+
+%% Format output file name
+
+data_to_save = struct();
+data_to_save.birdID = 'st1215';
+data_to_save.penetrationID = 1;
+data_to_save.siteID = 1;
+data_to_save.siteID = 1;
+data_to_save.XYpos = [0, 0];
+data_to_save.Zpos = 0;
+data_to_save.target_structure = 'NCM';
+
+data_to_save.fs = 31250;
+data_to_save.toedata = toedata;
+
+sav_date = datestr(now, 30);
+outfile_name = strcat('./st1215/good_data/pen2/more_merged/st1215_cat_P01_S01_2ndPen_moremerged_MUA', sav_date, '.mat');
+save(outfile_name, 'data_to_save');
+
+
+

toefile_utils/make_neurogram.m

@@ -0,0 +1,60 @@
+% make_neurogram.m
+% Brad Theilman September 2015
+
+clear
+close all
+
+load('st1215_cat_P01_S01_3rdPen_20150927T155818.mat');
+
+%TODO:  Get sampling frequency from data file.
+fs = 31250.0;
+
+nstims = 41;
+n_units = length(toedata);
+
+psth_winsize = 150; %ms
+psth_winsize_samps = floor((psth_winsize/1000)*fs);
+
+figure();
+for stim_index = 1:nstims
+    % unit_index is which cell you wish to plot
+
+    spike_density = zeros(n_units, 1);
+    for unit_index = 1:n_units
+
+        % Get the data for the chosen cell
+        unit_data = toedata{unit_index, 1};
+
+        % Collapse across all trials.  
+        trialdata = unit_data.stims{stim_index, 1};
+        ntrials = length(trialdata.toes);
+        spiketimes_across_trials = [];
+        for trial = 1:ntrials
+            spiketimes_across_trials = [spiketimes_across_trials; trialdata.toes{trial, 1}];
+        end
+        trial_length = floor(mean(trialdata.trial_end_times - trialdata.trial_start_times));
+        stim_start = floor(mean(trialdata.stim_start_times - trialdata.trial_start_times));
+        stim_end_secs = floor(mean(trialdata.stim_end_times - trialdata.stim_start_times))/fs;
+        winnum = 1;
+        for win_start = 1:psth_winsize_samps:trial_length
+
+            win_end = min(win_start+psth_winsize_samps, trial_length) - stim_start;
+            win_start_stim = win_start - stim_start;
+            nspikes_in_win = sum(spiketimes_across_trials <= win_end/fs & spiketimes_across_trials >= win_start_stim/fs);
+            spike_density(unit_index, winnum) = nspikes_in_win / (ntrials * (psth_winsize/1000));
+            win_starts(winnum) = win_start_stim / fs;
+            winnum = winnum+1;
+        end
+
+    end
+    max_firing_rate = max(max(spike_density))
+    subplot(7, 6, stim_index)
+    imagesc(win_starts, 1:n_units, spike_density, [0 85] );
+    colormap('gray');
+
+    line([0, 0], [0, n_units+1], 'Color', 'red')
+    line([stim_end_secs, stim_end_secs], [0, n_units+1], 'Color', 'red');
+end
+
+
+

toefile_utils/make_population_psth.m

@@ -0,0 +1,18 @@
+% make_population_psth
+
+clear
+close all
+
+load('st1215_cat_P01_S01_3rdPen_20150927T155818.mat');
+
+%TODO:  Get sampling frequency from data file.
+fs = 31250.0;
+
+nstims = 41;
+n_units = length(toedata);
+
+psth_winsize = 150; %ms
+psth_winsize_samps = floor((psth_winsize/1000)*fs);
+
+for stim_index = 1:nstims
+