From e97bc4a5e7ea10140beb08d1eda954cede4072d4 Mon Sep 17 00:00:00 2001 From: Bradley Theilman Date: Fri, 25 Sep 2015 13:40:16 -0700 Subject: [PATCH 01/12] Added matlab routines to deal with toefiles coming from klustakwik. These include routines to generate the toefiles and make spike rasters --- toefile_utils/kwik2mat.m | 154 ++++++++++++++++++++++++++ toefile_utils/make_raster.m | 35 ++++++ toefile_utils/make_raster_all_cells.m | 46 ++++++++ toefile_utils/make_raster_stim.m | 37 +++++++ 4 files changed, 272 insertions(+) create mode 100755 toefile_utils/kwik2mat.m create mode 100755 toefile_utils/make_raster.m create mode 100755 toefile_utils/make_raster_all_cells.m create mode 100755 toefile_utils/make_raster_stim.m diff --git a/toefile_utils/kwik2mat.m b/toefile_utils/kwik2mat.m new file mode 100755 index 0000000..dce57a6 --- /dev/null +++ b/toefile_utils/kwik2mat.m @@ -0,0 +1,154 @@ +% kwik to mat. Get good units out of sorted kwik file with stimulus +% information and store in one well formatted mat file + +clear +close all + +kwikfile = './st1215/man_sort/pen2/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 + if cluster_group == 2 + 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 +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; +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) + + disp(unit_num) + unit_entry = struct; + unit_entry.id = good_cluster_ids(unit_num); + stims = cell(nstims, 1); + + % Loop through each stimulus + for stim_num = 1:nstims + + disp(stim_num) + 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 + disp(trialnum) + 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; + + if any(spiketimes_secs_thistrial_relstimonset <0) + disp('PRE-STIM ACTIVITY DETECTED'); + end + + 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 + +sav_date = datestr(now, 30); +outfile_name = strcat('./st1215/good_data/st1215_cat_P01_S01_2ndPen_', sav_date, '.mat'); +save(outfile_name, 'toedata'); + + + diff --git a/toefile_utils/make_raster.m b/toefile_utils/make_raster.m new file mode 100755 index 0000000..ea38659 --- /dev/null +++ b/toefile_utils/make_raster.m @@ -0,0 +1,35 @@ +clear +close all + +% Make a spike raster for a given cell + + +load('st1215_cat_P01_S01_2ndPen_fixalignment_20150924T141526.mat'); +unit_index = 2; +fs = 31250.0 + +% Get the data for the chose cell +unit_data = toedata{unit_index, 1}; +nstims = length(unit_data.stims); + +figure(); +for stimnum = 1:nstims + subplot(7, 7, stimnum) + + stim_data = unit_data.stims{stimnum, 1}; + stim_end_secs = double(stim_data.stim_end_times - stim_data.stim_start_times)/fs; + ntrials = stim_data.ntrials; + for trialnum = 1:ntrials + ys = [0 + trialnum, 1+trialnum]; + if ~isempty(stim_data.toes{trialnum, 1}) + for spikenum = 1:length( stim_data.toes{trialnum, 1}) + line([stim_data.toes{trialnum, 1}(spikenum), stim_data.toes{trialnum, 1}(spikenum)], ys); + end + end + xlim([-2, stim_end_secs(trialnum)+2]) + ylim([0, ntrials+1]); + line([0, 0], [0, ntrials+1], 'Color', 'red'); + line([stim_end_secs(trialnum), stim_end_secs(trialnum)], [0, ntrials+1], 'Color', 'red'); + end +end + diff --git a/toefile_utils/make_raster_all_cells.m b/toefile_utils/make_raster_all_cells.m new file mode 100755 index 0000000..d37b6d9 --- /dev/null +++ b/toefile_utils/make_raster_all_cells.m @@ -0,0 +1,46 @@ +clear +close all + +% Make a spike raster for a given cell + +datafile = 'st1215_cat_P01_S01_2ndPen_20150915T165826'; + +load(strcat(datafile, '.mat')); +fs = 31250.0; + + +for unit_index = 1:31 + + disp(unit_index) + % Get the data for the chose cell + unit_data = toedata{unit_index, 1}; + nstims = length(unit_data.stims); + + figure(); + for stimnum = 1:nstims + subplot(7, 7, stimnum) + + stim_data = unit_data.stims{stimnum, 1}; + stim_end_secs = double(stim_data.stim_end_times - stim_data.stim_start_times)/fs; + ntrials = stim_data.ntrials; + for trialnum = 1:ntrials + ys = [0 + trialnum, 1+trialnum]; + if ~isempty(stim_data.toes{trialnum, 1}) + for spikenum = 1:length(stim_data.toes{trialnum, 1}) + line([stim_data.toes{trialnum, 1}(spikenum), stim_data.toes{trialnum, 1}(spikenum)], ys); + end + end + xlim([-2, stim_end_secs(trialnum)+2]) + ylim([0, ntrials+1]); + line([0, 0], [0, ntrials+1], 'Color', 'red'); + line([stim_end_secs(trialnum), stim_end_secs(trialnum)], [0, ntrials+1], 'Color', 'red'); + end + end + fig = gcf; + fig.PaperPositionMode = 'auto'; + figfilename = strcat(datafile, '_cell', num2str(unit_index)); + print(figfilename, '-dpng', '-r0'); + close all + + +end diff --git a/toefile_utils/make_raster_stim.m b/toefile_utils/make_raster_stim.m new file mode 100755 index 0000000..aa56219 --- /dev/null +++ b/toefile_utils/make_raster_stim.m @@ -0,0 +1,37 @@ +clear +close all + +% Make a spike raster for a given cell + + +load('st1215_cat_P01_S01_2ndPen_fixalignment_20150924T141526.mat'); +unit_index = 27; +fs = 31250.0 + +% Get the data for the chose cell +stimnum = 49 + +nstims = 49; + +figure(); +for unit_index= 1:31 + subplot(8, 4, unit_index) + + unit_data = toedata{unit_index, 1}; + stim_data = unit_data.stims{stimnum, 1}; + stim_end_secs = double(stim_data.stim_end_times - stim_data.stim_start_times)/fs; + ntrials = stim_data.ntrials; + for trialnum = 1:ntrials + ys = [0 + trialnum, 1+trialnum]; + if ~isempty(stim_data.toes{trialnum, 1}) + for spikenum = 1:length( stim_data.toes{trialnum, 1}) + line([stim_data.toes{trialnum, 1}(spikenum), stim_data.toes{trialnum, 1}(spikenum)], ys); + end + end + xlim([-2, stim_end_secs(trialnum)+2]) + ylim([0, ntrials+1]); + line([0, 0], [0, ntrials+1], 'Color', 'red'); + line([stim_end_secs(trialnum), stim_end_secs(trialnum)], [0, ntrials+1], 'Color', 'red'); + end +end + From b61848bff41ef61405feaddb7882b0ddfbea55e5 Mon Sep 17 00:00:00 2001 From: Bradley Theilman Date: Fri, 25 Sep 2015 13:55:54 -0700 Subject: [PATCH 02/12] Cleaned things up a little, added documentation for the raster routines --- toefile_utils/kwik2mat.m | 10 ++-------- toefile_utils/make_raster.m | 17 ++++++++++++----- toefile_utils/make_raster_all_cells.m | 7 ++++++- toefile_utils/make_raster_stim.m | 18 ++++++++++-------- 4 files changed, 30 insertions(+), 22 deletions(-) diff --git a/toefile_utils/kwik2mat.m b/toefile_utils/kwik2mat.m index dce57a6..057be02 100755 --- a/toefile_utils/kwik2mat.m +++ b/toefile_utils/kwik2mat.m @@ -1,5 +1,6 @@ % 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 @@ -81,7 +82,7 @@ %% Go through each cell and build final data matrix -fs = 31250.0; +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; @@ -93,7 +94,6 @@ for unit_num = 1:length(good_cluster_ids) - disp(unit_num) unit_entry = struct; unit_entry.id = good_cluster_ids(unit_num); stims = cell(nstims, 1); @@ -101,7 +101,6 @@ % Loop through each stimulus for stim_num = 1:nstims - disp(stim_num) stim_unit_entry = struct; this_stim_data = stim_data{stim_num, 1}; @@ -120,7 +119,6 @@ %start toes = cell(this_stim_data.ntrials, 1); for trialnum = 1:this_stim_data.ntrials - disp(trialnum) trial_start = trial_start_times(trialnum); trial_end = trial_end_times(trialnum); @@ -129,10 +127,6 @@ spiketimes_samps_thistrial_relstimonset = spiketimes_samps_thistrial - this_stim_data.start_times(trialnum); spiketimes_secs_thistrial_relstimonset = double(spiketimes_samps_thistrial_relstimonset) / fs; - if any(spiketimes_secs_thistrial_relstimonset <0) - disp('PRE-STIM ACTIVITY DETECTED'); - end - toes{trialnum, 1} = spiketimes_secs_thistrial_relstimonset; end stim_unit_entry.toes = toes; diff --git a/toefile_utils/make_raster.m b/toefile_utils/make_raster.m index ea38659..ab52c28 100755 --- a/toefile_utils/make_raster.m +++ b/toefile_utils/make_raster.m @@ -1,21 +1,28 @@ +%% make_raster.m +% Given a toefile, plot a spike raster plot for all trials of all stimuli +% for a single unit +% Brad Theilman September 2015 + clear close all -% Make a spike raster for a given cell - +%% Make a spike raster for a given cell load('st1215_cat_P01_S01_2ndPen_fixalignment_20150924T141526.mat'); + +%TODO: Get sampling frequency from data file. +fs = 31250.0; + +% unit_index is which cell you wish to plot unit_index = 2; -fs = 31250.0 -% Get the data for the chose cell +% Get the data for the chosen cell unit_data = toedata{unit_index, 1}; nstims = length(unit_data.stims); figure(); for stimnum = 1:nstims subplot(7, 7, stimnum) - stim_data = unit_data.stims{stimnum, 1}; stim_end_secs = double(stim_data.stim_end_times - stim_data.stim_start_times)/fs; ntrials = stim_data.ntrials; diff --git a/toefile_utils/make_raster_all_cells.m b/toefile_utils/make_raster_all_cells.m index d37b6d9..b078fc5 100755 --- a/toefile_utils/make_raster_all_cells.m +++ b/toefile_utils/make_raster_all_cells.m @@ -1,3 +1,8 @@ +%% make_raster_all_cells.m +% Given a toefile, plot a spike raster plot for all trials of all stimuli +% Does this for every unit in the file, saving the figures to png +% Brad Theilman September 2015 + clear close all @@ -37,7 +42,7 @@ end end fig = gcf; - fig.PaperPositionMode = 'auto'; + fig.PaperPositionMode = 'auto'; % TODO: Save in a large resolution figfilename = strcat(datafile, '_cell', num2str(unit_index)); print(figfilename, '-dpng', '-r0'); close all diff --git a/toefile_utils/make_raster_stim.m b/toefile_utils/make_raster_stim.m index aa56219..5f570c0 100755 --- a/toefile_utils/make_raster_stim.m +++ b/toefile_utils/make_raster_stim.m @@ -1,17 +1,19 @@ -clear -close all +%% make_raster.m +% Given a toefile, plot a spike raster plot for all trials of one stimulus +% for all units +% Brad Theilman September 2015 -% Make a spike raster for a given cell +clear +close all load('st1215_cat_P01_S01_2ndPen_fixalignment_20150924T141526.mat'); -unit_index = 27; -fs = 31250.0 -% Get the data for the chose cell -stimnum = 49 +fs = 31250.0; % TODO: get this from the data file + +stimnum = 49; -nstims = 49; +nstims = 49; % TODO: get this from the data file figure(); for unit_index= 1:31 From 6c5108f49e0ccad25efc087258deab2a6d850105 Mon Sep 17 00:00:00 2001 From: Bradley Theilman Date: Sun, 27 Sep 2015 16:18:17 -0700 Subject: [PATCH 03/12] Refactored the extraction of kwik data to conform to new format from klusta-pipeline Also fixed raster plotting to remove empty rows --- toefile_utils/kwik2mat.m | 32 +++- toefile_utils/kwik2mat.m~ | 164 ++++++++++++++++++ toefile_utils/make_raster.m | 8 +- toefile_utils/make_raster.m~ | 42 +++++ .../make_raster_singlecell_singlestim.m | 42 +++++ 5 files changed, 276 insertions(+), 12 deletions(-) create mode 100644 toefile_utils/kwik2mat.m~ create mode 100644 toefile_utils/make_raster.m~ create mode 100644 toefile_utils/make_raster_singlecell_singlestim.m diff --git a/toefile_utils/kwik2mat.m b/toefile_utils/kwik2mat.m index 057be02..211c389 100755 --- a/toefile_utils/kwik2mat.m +++ b/toefile_utils/kwik2mat.m @@ -5,7 +5,7 @@ clear close all -kwikfile = './st1215/man_sort/pen2/st1215_cat_P01_S01_2ndPen.kwik'; +kwikfile = './st1215/kwikfiles/pen3/st1215_cat_P01_S01_3rdPen.kwik'; % get information kwikinfo = h5info(kwikfile); @@ -47,17 +47,33 @@ end %% Now get stimulus information -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'); + +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; @@ -141,7 +157,7 @@ %% Format output file name sav_date = datestr(now, 30); -outfile_name = strcat('./st1215/good_data/st1215_cat_P01_S01_2ndPen_', sav_date, '.mat'); +outfile_name = strcat('./st1215/good_data/st1215_cat_P01_S01_3rdPen_', sav_date, '.mat'); save(outfile_name, 'toedata'); diff --git a/toefile_utils/kwik2mat.m~ b/toefile_utils/kwik2mat.m~ new file mode 100644 index 0000000..9a62c7b --- /dev/null +++ b/toefile_utils/kwik2mat.m~ @@ -0,0 +1,164 @@ +% 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/man_sort/pen3/st1215_cat_P01_S01_3rdPen.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 + if cluster_group == 2 + 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_filenames = cell2mat(stim_text); +stim_start_filenames = stim_start_filenames(2:2:end, :); +stim_start_filenames = cellstr(stim_start_filenames); + +% 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 + +sav_date = datestr(now, 30); +outfile_name = strcat('./st1215/good_data/st1215_cat_P01_S01_2ndPen_', sav_date, '.mat'); +save(outfile_name, 'toedata'); + + + diff --git a/toefile_utils/make_raster.m b/toefile_utils/make_raster.m index ab52c28..85e4bca 100755 --- a/toefile_utils/make_raster.m +++ b/toefile_utils/make_raster.m @@ -8,13 +8,13 @@ %% Make a spike raster for a given cell -load('st1215_cat_P01_S01_2ndPen_fixalignment_20150924T141526.mat'); +load('st1215_cat_P01_S01_3rdPen_20150927T155818.mat'); %TODO: Get sampling frequency from data file. fs = 31250.0; % unit_index is which cell you wish to plot -unit_index = 2; +unit_index = 9; % Get the data for the chosen cell unit_data = toedata{unit_index, 1}; @@ -27,14 +27,14 @@ stim_end_secs = double(stim_data.stim_end_times - stim_data.stim_start_times)/fs; ntrials = stim_data.ntrials; for trialnum = 1:ntrials - ys = [0 + trialnum, 1+trialnum]; + ys = [trialnum-1, trialnum]; if ~isempty(stim_data.toes{trialnum, 1}) for spikenum = 1:length( stim_data.toes{trialnum, 1}) line([stim_data.toes{trialnum, 1}(spikenum), stim_data.toes{trialnum, 1}(spikenum)], ys); end end xlim([-2, stim_end_secs(trialnum)+2]) - ylim([0, ntrials+1]); + ylim([0, ntrials]); line([0, 0], [0, ntrials+1], 'Color', 'red'); line([stim_end_secs(trialnum), stim_end_secs(trialnum)], [0, ntrials+1], 'Color', 'red'); end diff --git a/toefile_utils/make_raster.m~ b/toefile_utils/make_raster.m~ new file mode 100644 index 0000000..25e3888 --- /dev/null +++ b/toefile_utils/make_raster.m~ @@ -0,0 +1,42 @@ +%% make_raster.m +% Given a toefile, plot a spike raster plot for all trials of all stimuli +% for a single unit +% Brad Theilman September 2015 + +clear +close all + +%% Make a spike raster for a given cell + +load('st1215_cat_P01_S01_3rdPen_20150927T155818.mat'); + +%TODO: Get sampling frequency from data file. +fs = 31250.0; + +% unit_index is which cell you wish to plot +unit_index = 4; + +% Get the data for the chosen cell +unit_data = toedata{unit_index, 1}; +nstims = length(unit_data.stims); + +figure(); +for stimnum = 1:nstims + subplot(7, 7, stimnum) + stim_data = unit_data.stims{stimnum, 1}; + stim_end_secs = double(stim_data.stim_end_times - stim_data.stim_start_times)/fs; + ntrials = stim_data.ntrials; + for trialnum = 1:ntrials + ys = [trialnum-1, trialnum]; + if ~isempty(stim_data.toes{trialnum, 1}) + for spikenum = 1:length( stim_data.toes{trialnum, 1}) + line([stim_data.toes{trialnum, 1}(spikenum), stim_data.toes{trialnum, 1}(spikenum)], ys); + end + end + xlim([-2, stim_end_secs(trialnum)+2]) + ylim([0, ntrials+1]); + line([0, 0], [0, ntrials+1], 'Color', 'red'); + line([stim_end_secs(trialnum), stim_end_secs(trialnum)], [0, ntrials+1], 'Color', 'red'); + end +end + diff --git a/toefile_utils/make_raster_singlecell_singlestim.m b/toefile_utils/make_raster_singlecell_singlestim.m new file mode 100644 index 0000000..cc3e6ac --- /dev/null +++ b/toefile_utils/make_raster_singlecell_singlestim.m @@ -0,0 +1,42 @@ +%% make_raster.m +% Given a toefile, plot a spike raster plot for all trials of all stimuli +% for a single unit +% Brad Theilman September 2015 + +clear +close all + +%% Make a spike raster for a given cell + +load('st1215_cat_P01_S01_3rdPen_20150927T155818.mat'); + +%TODO: Get sampling frequency from data file. +fs = 31250.0; + +% unit_index is which cell you wish to plot +unit_index = 9; + +% Get the data for the chosen cell +unit_data = toedata{unit_index, 1}; +nstims = length(unit_data.stims); + +figure(); +stimnum = 10; + +stim_data = unit_data.stims{stimnum, 1}; +stim_end_secs = double(stim_data.stim_end_times - stim_data.stim_start_times)/fs; +ntrials = stim_data.ntrials; +for trialnum = 1:ntrials + ys = [trialnum-1, trialnum]; + if ~isempty(stim_data.toes{trialnum, 1}) + for spikenum = 1:length( stim_data.toes{trialnum, 1}) + line([stim_data.toes{trialnum, 1}(spikenum), stim_data.toes{trialnum, 1}(spikenum)], ys); + end + end + xlim([-2, stim_end_secs(trialnum)+2]) + ylim([0, ntrials]); + line([0, 0], [0, ntrials+1], 'Color', 'red'); + line([stim_end_secs(trialnum), stim_end_secs(trialnum)], [0, ntrials+1], 'Color', 'red'); +end + + From 7385bd533f19f9cee63554447698db8002a8a9da Mon Sep 17 00:00:00 2001 From: Bradley Theilman Date: Sun, 27 Sep 2015 17:04:07 -0700 Subject: [PATCH 04/12] Fixed figure printing to print at a decent size/resolution --- toefile_utils/make_raster_all_cells.m | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/toefile_utils/make_raster_all_cells.m b/toefile_utils/make_raster_all_cells.m index b078fc5..40d3c69 100755 --- a/toefile_utils/make_raster_all_cells.m +++ b/toefile_utils/make_raster_all_cells.m @@ -8,7 +8,7 @@ % Make a spike raster for a given cell -datafile = 'st1215_cat_P01_S01_2ndPen_20150915T165826'; +datafile = './pen2/st1215_cat_P01_S01_2ndPen_fixalignment_20150924T141526'; load(strcat(datafile, '.mat')); fs = 31250.0; @@ -21,7 +21,8 @@ unit_data = toedata{unit_index, 1}; nstims = length(unit_data.stims); - figure(); + fig = gcf; + fig.Visible='off' for stimnum = 1:nstims subplot(7, 7, stimnum) @@ -29,22 +30,24 @@ stim_end_secs = double(stim_data.stim_end_times - stim_data.stim_start_times)/fs; ntrials = stim_data.ntrials; for trialnum = 1:ntrials - ys = [0 + trialnum, 1+trialnum]; + ys = [trialnum-1, trialnum]; if ~isempty(stim_data.toes{trialnum, 1}) for spikenum = 1:length(stim_data.toes{trialnum, 1}) line([stim_data.toes{trialnum, 1}(spikenum), stim_data.toes{trialnum, 1}(spikenum)], ys); end end xlim([-2, stim_end_secs(trialnum)+2]) - ylim([0, ntrials+1]); + ylim([0, ntrials]); line([0, 0], [0, ntrials+1], 'Color', 'red'); line([stim_end_secs(trialnum), stim_end_secs(trialnum)], [0, ntrials+1], 'Color', 'red'); end end - fig = gcf; - fig.PaperPositionMode = 'auto'; % TODO: Save in a large resolution + + fig.PaperUnits = 'inches'; + fig.PaperPosition = [0 0 11 8.5]; + fig.PaperPositionMode = 'manual'; % TODO: Save in a large resolution figfilename = strcat(datafile, '_cell', num2str(unit_index)); - print(figfilename, '-dpng', '-r0'); + print(figfilename, '-dpng', '-r600'); close all From 0115474303686c8bcd104fbc8507c9f9a559fb06 Mon Sep 17 00:00:00 2001 From: Bradley Theilman Date: Mon, 28 Sep 2015 22:40:20 -0700 Subject: [PATCH 05/12] Can now change colors of plot components Just cause. --- toefile_utils/make_raster.m | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/toefile_utils/make_raster.m b/toefile_utils/make_raster.m index 85e4bca..bc372fb 100755 --- a/toefile_utils/make_raster.m +++ b/toefile_utils/make_raster.m @@ -20,9 +20,13 @@ unit_data = toedata{unit_index, 1}; nstims = length(unit_data.stims); +% Plotting parameters +spike_color = 'blue'; +stim_color = 'red'; + figure(); for stimnum = 1:nstims - subplot(7, 7, stimnum) + subplot(7, 6, stimnum) stim_data = unit_data.stims{stimnum, 1}; stim_end_secs = double(stim_data.stim_end_times - stim_data.stim_start_times)/fs; ntrials = stim_data.ntrials; @@ -30,13 +34,13 @@ ys = [trialnum-1, trialnum]; if ~isempty(stim_data.toes{trialnum, 1}) for spikenum = 1:length( stim_data.toes{trialnum, 1}) - line([stim_data.toes{trialnum, 1}(spikenum), stim_data.toes{trialnum, 1}(spikenum)], ys); + line([stim_data.toes{trialnum, 1}(spikenum), stim_data.toes{trialnum, 1}(spikenum)], ys, 'Color', spike_color); end end xlim([-2, stim_end_secs(trialnum)+2]) ylim([0, ntrials]); - line([0, 0], [0, ntrials+1], 'Color', 'red'); - line([stim_end_secs(trialnum), stim_end_secs(trialnum)], [0, ntrials+1], 'Color', 'red'); + line([0, 0], [0, ntrials+1], 'Color', stim_color); + line([stim_end_secs(trialnum), stim_end_secs(trialnum)], [0, ntrials+1], 'Color', stim_color); end end From f3a314b2f0224d7741da5b2098c9a844cbfa37bb Mon Sep 17 00:00:00 2001 From: Bradley Theilman Date: Wed, 30 Sep 2015 12:36:55 -0700 Subject: [PATCH 06/12] Make neurogram added. Definitely could be improved.. --- toefile_utils/kwik2mat.m~ | 164 -------------------------------- toefile_utils/make_neurogram.m | 60 ++++++++++++ toefile_utils/make_neurogram.m~ | 56 +++++++++++ toefile_utils/make_raster.m~ | 42 -------- 4 files changed, 116 insertions(+), 206 deletions(-) delete mode 100644 toefile_utils/kwik2mat.m~ create mode 100644 toefile_utils/make_neurogram.m create mode 100644 toefile_utils/make_neurogram.m~ delete mode 100644 toefile_utils/make_raster.m~ diff --git a/toefile_utils/kwik2mat.m~ b/toefile_utils/kwik2mat.m~ deleted file mode 100644 index 9a62c7b..0000000 --- a/toefile_utils/kwik2mat.m~ +++ /dev/null @@ -1,164 +0,0 @@ -% 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/man_sort/pen3/st1215_cat_P01_S01_3rdPen.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 - if cluster_group == 2 - 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_filenames = cell2mat(stim_text); -stim_start_filenames = stim_start_filenames(2:2:end, :); -stim_start_filenames = cellstr(stim_start_filenames); - -% 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 - -sav_date = datestr(now, 30); -outfile_name = strcat('./st1215/good_data/st1215_cat_P01_S01_2ndPen_', sav_date, '.mat'); -save(outfile_name, 'toedata'); - - - diff --git a/toefile_utils/make_neurogram.m b/toefile_utils/make_neurogram.m new file mode 100644 index 0000000..88e384f --- /dev/null +++ b/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 + + + \ No newline at end of file diff --git a/toefile_utils/make_neurogram.m~ b/toefile_utils/make_neurogram.m~ new file mode 100644 index 0000000..a67744e --- /dev/null +++ b/toefile_utils/make_neurogram.m~ @@ -0,0 +1,56 @@ +% 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 = 25; %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)); + 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); + win_starts(winnum) = win_start_stim / fs; + winnum = winnum+1; + end + + end + + subplot(7, 6, stim_index) + imagesc(win_starts, 1:n_units, spike_density); + line([0, 0], [0, n_units], +end + + + \ No newline at end of file diff --git a/toefile_utils/make_raster.m~ b/toefile_utils/make_raster.m~ deleted file mode 100644 index 25e3888..0000000 --- a/toefile_utils/make_raster.m~ +++ /dev/null @@ -1,42 +0,0 @@ -%% make_raster.m -% Given a toefile, plot a spike raster plot for all trials of all stimuli -% for a single unit -% Brad Theilman September 2015 - -clear -close all - -%% Make a spike raster for a given cell - -load('st1215_cat_P01_S01_3rdPen_20150927T155818.mat'); - -%TODO: Get sampling frequency from data file. -fs = 31250.0; - -% unit_index is which cell you wish to plot -unit_index = 4; - -% Get the data for the chosen cell -unit_data = toedata{unit_index, 1}; -nstims = length(unit_data.stims); - -figure(); -for stimnum = 1:nstims - subplot(7, 7, stimnum) - stim_data = unit_data.stims{stimnum, 1}; - stim_end_secs = double(stim_data.stim_end_times - stim_data.stim_start_times)/fs; - ntrials = stim_data.ntrials; - for trialnum = 1:ntrials - ys = [trialnum-1, trialnum]; - if ~isempty(stim_data.toes{trialnum, 1}) - for spikenum = 1:length( stim_data.toes{trialnum, 1}) - line([stim_data.toes{trialnum, 1}(spikenum), stim_data.toes{trialnum, 1}(spikenum)], ys); - end - end - xlim([-2, stim_end_secs(trialnum)+2]) - ylim([0, ntrials+1]); - line([0, 0], [0, ntrials+1], 'Color', 'red'); - line([stim_end_secs(trialnum), stim_end_secs(trialnum)], [0, ntrials+1], 'Color', 'red'); - end -end - From bcfe23097f3b1037d9772ef5be0de5a771973712 Mon Sep 17 00:00:00 2001 From: Bradley Theilman Date: Thu, 3 Dec 2015 10:21:06 -0800 Subject: [PATCH 07/12] Up to speed --- .DS_Store | Bin 0 -> 6148 bytes toefile_utils/kwik2mat.m | 7 ++-- toefile_utils/make_neurogram.m~ | 56 -------------------------- toefile_utils/make_population_psth.m | 18 +++++++++ toefile_utils/make_raster_all_cells.m | 2 +- 5 files changed, 23 insertions(+), 60 deletions(-) create mode 100644 .DS_Store delete mode 100644 toefile_utils/make_neurogram.m~ create mode 100644 toefile_utils/make_population_psth.m diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..348ceb207a9ab9c0e52367aa261c035f4e11265b GIT binary patch literal 6148 zcmeHK%}N6?5Kd~ODCeq@f#VywOfGIlCW+UY<&IDPJ&)?6+|+CsJGxM5=mMtzVh5uetKrsD$HS7 zESz~Q+4E9A?IcydeWF^T83x9#`)O-vA4(aGY(F>>VXt8?@5wmzMI5%rQcsUX+?2hl z?8MEaVXxHNVHcw7ejFZ&L9tkB(y&+i{c>e{!zpi-2OCbmQYnI6+8GQiYi+%> zdvF~0Z*K4IA0D5c^&)WiUX?5v9Ks72-|*4(+HoY~ONbhdFw7(*28aP-;O{VCb}(!H z?>sU+DKS6{{15}UKRBR>uE9*B+B%@Z-$(4vBcgzfcL_wtpldMG2t6QNrvmC!Zmt+y zr-R=yajwBkqfTdB%?#t1nYp>4a5X#l9Tm>FtC4DAfEf75z?5#eEyegFTML@i=~ z82DEV@O;gyRpFJ)**fuZc-A_g=b$JUmuVcOfT50Jh{dCL6;uiM9W(%4gPBI~fY6VC MqJb)6;71vF2Yv%#@c;k- literal 0 HcmV?d00001 diff --git a/toefile_utils/kwik2mat.m b/toefile_utils/kwik2mat.m index 211c389..93fcba3 100755 --- a/toefile_utils/kwik2mat.m +++ b/toefile_utils/kwik2mat.m @@ -5,7 +5,7 @@ clear close all -kwikfile = './st1215/kwikfiles/pen3/st1215_cat_P01_S01_3rdPen.kwik'; +kwikfile = './st1215/kwdfiles/pen2_latestpipeline/st1215_cat_P01_S01_2ndPen.kwik'; % get information kwikinfo = h5info(kwikfile); @@ -27,7 +27,8 @@ cluster_group = h5readatt(kwikfile, cluster_attr_path, 'cluster_group'); disp(cluster_group); %cluster_group 2 is Good - if cluster_group == 2 + % 1 is MUA + if cluster_group == 1 good_cluster_inds(i) = 1; end @@ -157,7 +158,7 @@ %% Format output file name sav_date = datestr(now, 30); -outfile_name = strcat('./st1215/good_data/st1215_cat_P01_S01_3rdPen_', sav_date, '.mat'); +outfile_name = strcat('./st1215/good_data/pen2/more_merged/st1215_cat_P01_S01_2ndPen_moremerged_MUA', sav_date, '.mat'); save(outfile_name, 'toedata'); diff --git a/toefile_utils/make_neurogram.m~ b/toefile_utils/make_neurogram.m~ deleted file mode 100644 index a67744e..0000000 --- a/toefile_utils/make_neurogram.m~ +++ /dev/null @@ -1,56 +0,0 @@ -% 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 = 25; %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)); - 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); - win_starts(winnum) = win_start_stim / fs; - winnum = winnum+1; - end - - end - - subplot(7, 6, stim_index) - imagesc(win_starts, 1:n_units, spike_density); - line([0, 0], [0, n_units], -end - - - \ No newline at end of file diff --git a/toefile_utils/make_population_psth.m b/toefile_utils/make_population_psth.m new file mode 100644 index 0000000..a7634fd --- /dev/null +++ b/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 + \ No newline at end of file diff --git a/toefile_utils/make_raster_all_cells.m b/toefile_utils/make_raster_all_cells.m index 40d3c69..5b9897d 100755 --- a/toefile_utils/make_raster_all_cells.m +++ b/toefile_utils/make_raster_all_cells.m @@ -8,7 +8,7 @@ % Make a spike raster for a given cell -datafile = './pen2/st1215_cat_P01_S01_2ndPen_fixalignment_20150924T141526'; +datafile = './pen2/more_merged/st1215_cat_P01_S01_2ndPen_moremerged_MUA20151102T092630'; load(strcat(datafile, '.mat')); fs = 31250.0; From 24cb20ca861b1714107608e7fda8e7f17f633c66 Mon Sep 17 00:00:00 2001 From: Bradley Theilman Date: Thu, 3 Dec 2015 10:25:42 -0800 Subject: [PATCH 08/12] Added more fields to save metadata --- toefile_utils/kwik2mat.m | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/toefile_utils/kwik2mat.m b/toefile_utils/kwik2mat.m index 93fcba3..403b435 100755 --- a/toefile_utils/kwik2mat.m +++ b/toefile_utils/kwik2mat.m @@ -157,9 +157,21 @@ %% 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, 'toedata'); +save(outfile_name, 'data_to_save'); From d751333a68930ab22139dd4b202b4f46a702df5d Mon Sep 17 00:00:00 2001 From: Bradley Theilman Date: Thu, 3 Dec 2015 11:08:26 -0800 Subject: [PATCH 09/12] Generalizing the data structure --- toefile_utils/generate_population_vectors.m | 80 +++++++++++++++++++++ 1 file changed, 80 insertions(+) create mode 100644 toefile_utils/generate_population_vectors.m diff --git a/toefile_utils/generate_population_vectors.m b/toefile_utils/generate_population_vectors.m new file mode 100644 index 0000000..dd4a613 --- /dev/null +++ b/toefile_utils/generate_population_vectors.m @@ -0,0 +1,80 @@ +%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 = spiketimes{trial, 1}; + 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, 'popvec_data_to_save'); +clear +close all From 6398eb65ed8f8df1795ddc3ff4f7b3ccbf019412 Mon Sep 17 00:00:00 2001 From: Bradley Theilman Date: Mon, 7 Dec 2015 15:31:14 -0800 Subject: [PATCH 10/12] Calculating cell groups --- toefile_utils/generate_population_vectors.m | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/toefile_utils/generate_population_vectors.m b/toefile_utils/generate_population_vectors.m index dd4a613..2db4068 100644 --- a/toefile_utils/generate_population_vectors.m +++ b/toefile_utils/generate_population_vectors.m @@ -76,5 +76,9 @@ popvec_data_to_save.winstarts = windows_starts; popvecfilename = strcat(datafilename, '_popvec'); save(popvecfilename, 'popvec_data_to_save'); -clear -close all + +% From population vectors, generate cell groups + +%first compute average firing rate for cell: +avg_fr = sum(sum(sum(population_vectors, 4), 3), 2)/(nwin*nstim*n_bin_start); + From 395ed36d0b674f6a4e35c5953da92286e694523d Mon Sep 17 00:00:00 2001 From: Bradley Theilman Date: Mon, 7 Dec 2015 17:13:35 -0800 Subject: [PATCH 11/12] Cell groups --- toefile_utils/generate_cell_groups.m | 54 +++++++++++++++++++++ toefile_utils/generate_population_vectors.m | 18 +++---- 2 files changed, 61 insertions(+), 11 deletions(-) create mode 100644 toefile_utils/generate_cell_groups.m diff --git a/toefile_utils/generate_cell_groups.m b/toefile_utils/generate_cell_groups.m new file mode 100644 index 0000000..634a033 --- /dev/null +++ b/toefile_utils/generate_cell_groups.m @@ -0,0 +1,54 @@ +% From population vectors, generate cell groups +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) \ No newline at end of file diff --git a/toefile_utils/generate_population_vectors.m b/toefile_utils/generate_population_vectors.m index 2db4068..6bc9a41 100644 --- a/toefile_utils/generate_population_vectors.m +++ b/toefile_utils/generate_population_vectors.m @@ -7,7 +7,7 @@ datafilename = 'st1215_cat_P01_S01_2ndPen_moremerged_20151102T090824'; load(strcat(datafilename, '.mat')); -toedata = data_to_save.toedata; +%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. @@ -52,7 +52,7 @@ for bin_start = 1:n_bin_start spikecounts_thiscellstimbinst = zeros(1, 500); for trial = 1:ntrials - spikes_this_trial = spiketimes{trial, 1}; + 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]); @@ -70,15 +70,11 @@ 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; +% 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, 'popvec_data_to_save'); +save(popvecfilename, 'population_vectors', 'bin_width', 'windows_starts'); -% From population vectors, generate cell groups - -%first compute average firing rate for cell: -avg_fr = sum(sum(sum(population_vectors, 4), 3), 2)/(nwin*nstim*n_bin_start); From 5f512120b8d183861c7004dd03f4a5a0e3672f52 Mon Sep 17 00:00:00 2001 From: Bradley Theilman Date: Thu, 10 Dec 2015 11:34:01 -0800 Subject: [PATCH 12/12] Cleaning up --- toefile_utils/generate_cell_groups.m | 9 ++------- toefile_utils/generate_population_vectors.m | 5 ----- 2 files changed, 2 insertions(+), 12 deletions(-) diff --git a/toefile_utils/generate_cell_groups.m b/toefile_utils/generate_cell_groups.m index 634a033..b774966 100644 --- a/toefile_utils/generate_cell_groups.m +++ b/toefile_utils/generate_cell_groups.m @@ -1,4 +1,5 @@ % From population vectors, generate cell groups +% Must activate Plex library before running clear close all @@ -6,12 +7,12 @@ %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 @@ -19,23 +20,18 @@ % 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)) @@ -43,7 +39,6 @@ stream.addElement(cells_in_this_group); stream.ensureAllFaces(); end - end end end diff --git a/toefile_utils/generate_population_vectors.m b/toefile_utils/generate_population_vectors.m index 6bc9a41..2564eb8 100644 --- a/toefile_utils/generate_population_vectors.m +++ b/toefile_utils/generate_population_vectors.m @@ -25,8 +25,6 @@ % 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); @@ -38,7 +36,6 @@ end % build population vectors - ncells = length(toedata); nstim = 41; population_vectors = zeros(ncells, nwin, nstim, n_bin_start); @@ -57,11 +54,9 @@ 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;