From 1ddf2d7f383706a47549c76d873cd05cf54fe4ea Mon Sep 17 00:00:00 2001 From: Konstantinos Date: Thu, 25 Apr 2019 22:29:41 -0400 Subject: [PATCH 1/5] Buzsaki functions for NWB file format --- NWB/bz_GetLFP_NWB.m | 201 +++++++++ NWB/bz_GetSpikes_NWB.m | 245 +++++++++++ NWB/bz_LoadBehavior_NWB.m | 436 +++++++++++++++++++ NWB/bz_LoadEvents_NWB.m | 121 +++++ NWB/get_SessionInfoFromNWB.m | 181 ++++++++ NWB/testing_BUZSAKI_NWB_functions_fidelity.m | 130 ++++++ 6 files changed, 1314 insertions(+) create mode 100644 NWB/bz_GetLFP_NWB.m create mode 100644 NWB/bz_GetSpikes_NWB.m create mode 100644 NWB/bz_LoadBehavior_NWB.m create mode 100644 NWB/bz_LoadEvents_NWB.m create mode 100644 NWB/get_SessionInfoFromNWB.m create mode 100644 NWB/testing_BUZSAKI_NWB_functions_fidelity.m diff --git a/NWB/bz_GetLFP_NWB.m b/NWB/bz_GetLFP_NWB.m new file mode 100644 index 00000000..28e75894 --- /dev/null +++ b/NWB/bz_GetLFP_NWB.m @@ -0,0 +1,201 @@ +function [lfp] = bz_GetLFP_NWB(varargin) +% bz_GetLFP - Get local field potentials. + +% Load local field potentials from NWB file + + +% USAGE +% +% [lfp] = bz_GetLFP(channels,) +% +% INPUTS +% +% channels(required) -must be first input, numeric +% list of channels to load (use keyword 'all' for all) +% channID is 0-indexing, a la neuroscope +% Name-value paired inputs: + +% 'intervals' -list of time intervals [0 10; 20 30] to read from +% the LFP file (default is [0 inf]) +% 'downsample' -factor to downsample the LFP (i.e. 'downsample',5 +% will load a 1250Hz .lfp file at 250Hz) +% +% OUTPUT +% +% lfp struct of lfp data. Can be a single struct or an array +% of structs for different intervals. lfp(1), lfp(2), +% etc for intervals(1,:), intervals(2,:), etc +% .data [Nt x Nd] matrix of the LFP data +% .timestamps [Nt x 1] vector of timestamps to match LFP data +% .interval [1 x 2] vector of start/stop times of LFP interval +% .channels [Nd X 1] vector of channel ID's +% .samplingRate LFP sampling rate [default = 1250] +% .duration duration, in seconds, of LFP interval + +% Copyright (C) 2004-2011 by Michaël Zugaro +% edited by David Tingley, 2017 + + +%% EXAMPLE CALLS +% lfp_NWB = bz_GetLFP_NWB(5); % Loads the entire signal from channel 5 (#6) from the current folder +% lfp_NWB = bz_GetLFP_NWB(5, 'nwb_file', nwb_file); % Loads the entire signal from channel 5 (#6) from the current folder from the specified NWB file +% lfp_NWB = bz_GetLFP_NWB([5 10], 'nwb_file', nwb_file); % Loads channels 5 (#6) and 10 (#11) +% lfp_NWB = bz_GetLFP_NWB(5, 'restrict',[0 120]); % Loads channel 5 (#6) between [0, 120] seconds +% lfp_NWB = bz_GetLFP_NWB(5, 'restrict',[0 120;240.2 265.23]); % Loads channel 5 (#6) between [0, 120] and [240.2 265.23] seconds +% lfp_NWB = bz_GetLFP_NWB(5, 'restrict',[0 120], 'downsample', 2); % Loads channel 5 (#6) between [0, 120] seconds and downsamples by a factor of 2 + + +% Added support for NWB: Konstantinos Nasiotis 2019 + + +%% Parse the inputs! + +channelsValidation = @(x) isnumeric(x) || strcmp(x,'all'); + +% parse args +p = inputParser; +addParameter(p,'nwb_file','',@isstr); +addRequired(p,'channels',channelsValidation) +addParameter(p,'basename','',@isstr) +addParameter(p,'intervals',[],@isnumeric) +addParameter(p,'restrict',[],@isnumeric) +addParameter(p,'basepath',pwd,@isstr); +addParameter(p,'downsample',1,@isnumeric); +addParameter(p,'saveMat',false,@islogical); +addParameter(p,'forceReload',false,@islogical); +addParameter(p,'noPrompts',false,@islogical); + +parse(p,varargin{:}) +nwb_file = p.Results.nwb_file; +channels = p.Results.channels; +downsamplefactor = p.Results.downsample; +basepath = p.Results.basepath; +noPrompts = p.Results.noPrompts; + +% doing this so you can use either 'intervals' or 'restrict' as parameters to do the same thing +intervals = p.Results.intervals; +restrict = p.Results.restrict; +if isempty(intervals) && isempty(restrict) % both empty + intervals = [0 Inf]; +elseif isempty(intervals) && ~isempty(restrict) % intervals empty, restrict isn't + intervals = restrict; +end + +[filepath,name,ext] = fileparts(nwb_file); + + +%% let's check that there is an appropriate NWB file +if isempty(nwb_file) + %disp('No nwb_file given, so we look for a *nwb file in the current directory') + d = dir('*nwb'); + if length(d) > 1 % we assume one .nwb file or this should break + error('there is more than one .nwb file in this directory?'); + elseif length(d) == 0 + d = dir('*nwb'); + if isempty(d) + error('could not find an nwb file..') + end + end + nwb_file = fullfile(d.folder, d.name); + + lfp.Filename = d.name; +else + lfp.Filename = [name ext]; +end + +%% things we can parse from sessionInfo or xml file + +sessionInfo = bz_getSessionInfo(basepath, 'noPrompts', noPrompts); + +try + samplingRate = sessionInfo.lfpSampleRate; +catch + samplingRate = sessionInfo.rates.lfp; % old ugliness we need to get rid of +end +samplingRateLFP = samplingRate./downsamplefactor; + +if mod(samplingRateLFP,1)~=0 + error('samplingRate/downsamplefactor must be an integer') +end +%% Channel load options +%Right now this assumes that all means channels 0:nunchannels-1 (neuroscope +%indexing), we could also add options for this to be select region or spike +%group from the xml... +if strcmp(channels,'all') + channels = sessionInfo.channels; +end + +%% get the data +disp('loading LFP file...') +nIntervals = size(intervals,1); + + +nwb2 = nwbRead(nwb_file); + + +try + % Check if the data is in LFP format + all_lfp_keys = keys(nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries); + + for iKey = 1:length(all_lfp_keys) + if ismember(all_lfp_keys{iKey}, {'lfp','bla bla bla'}) %%%%%%%% ADD MORE HERE, DON'T KNOW WHAT THE STANDARD FORMATS ARE + iLFPDataKey = iKey; + LFPDataPresent = 1; + break % Once you find the data don't look for other keys/trouble + else + LFPDataPresent = 0; + end + end +catch + LFPDataPresent = 0; +end + +if ~LFPDataPresent + error ('No LFP signals exist in this .nwb file') +end + + + + + +% returns lfp/bz format +for i = 1:nIntervals + lfp(i).duration = (intervals(i,2)-intervals(i,1)); + lfp(i).interval = [intervals(i,1) intervals(i,2)]; + + % Load data and put into struct + % we assume 0-indexing like neuroscope, but bz_LoadBinary uses 1-indexing to + % load.... + + if intervals(i,2) == Inf + use_this_bracket = sessionInfo.samples_NWB/sessionInfo.lfpSampleRate; + data_temp = zeros(length(channels),ceil((use_this_bracket - intervals(i,1))*samplingRateLFP)); + else + use_this_bracket = intervals(i,2); + data_temp = zeros(length(channels),floor((intervals(i,2) - intervals(i,1))*samplingRateLFP)); + end + + % This is not optimized yet + + for iChannel = 1:length(channels) + disp(['Now concatenating channel: ' num2str(channels(iChannel))]) + data_temp(iChannel,:) = nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries.get(all_lfp_keys{iLFPDataKey}).data.load([channels(iChannel)+1, intervals(i,1)*samplingRate+1],[1, downsamplefactor],[channels(iChannel)+1, use_this_bracket*samplingRate]); + end + + lfp(i).data = data_temp'; clear data_temp + lfp(i).timestamps = [lfp(i).interval(1):(1/samplingRateLFP):... + (lfp(i).interval(1)+(length(lfp(i).data)-1)/... + samplingRateLFP)]'; + lfp(i).channels = channels; + lfp(i).samplingRate = samplingRateLFP; + % check if duration is inf, and reset to actual duration... + if lfp(i).interval(2) == inf + lfp(i).interval(2) = length(lfp(i).timestamps)/lfp(i).samplingRate; + lfp(i).duration = (lfp(i).interval(i,2)-lfp(i).interval(i,1)); + end + + if isfield(sessionInfo,'region') && isfield(sessionInfo,'channels') + [~,~,regionidx] = intersect(lfp(i).channels,sessionInfo.channels,'stable'); + lfp(i).region = sessionInfo.region(regionidx); % match region order to channel order.. + end +end diff --git a/NWB/bz_GetSpikes_NWB.m b/NWB/bz_GetSpikes_NWB.m new file mode 100644 index 00000000..525674b3 --- /dev/null +++ b/NWB/bz_GetSpikes_NWB.m @@ -0,0 +1,245 @@ +function spikes = bz_GetSpikes_NWB(varargin) +% bz_getSpikes - Get spike timestamps. +% +% USAGE +% +% spikes = bz_getSpikes(varargin) +% +% INPUTS +% +% spikeGroups -vector subset of shank IDs to load (Default: all) +% region -string region ID to load neurons from specific region +% (requires sessionInfo file or units->structures in xml) +% UID -vector subset of UID's to load +% getWaveforms -logical (default=true) to load mean of raw waveform data +% forceReload -logical (default=false) to force loading from +% res/clu/spk files +% saveMat -logical (default=false) to save in buzcode format +% noPrompts -logical (default=false) to supress any user prompts + + +% spikes_NWB = bz_GetSpikes_NWB('UID',[2 4 7]); % Loads neurons [2, 4 7] from an NWB at the current directory +% spikes_NWB = bz_GetSpikes_NWB('nwb_file', nwb_file, 'UID',[2 4 7]); % Loads neurons [2, 4 7] from a specific NWB file +% spikes_NWB = bz_GetSpikes_NWB('nwb_file', nwb_file, 'spikeGroups', [2,4]); % Loads the neurons from a specific NWB file and spikeGroups [2,4] +% spikes_NWB = bz_GetSpikes_NWB('nwb_file', nwb_file, 'region', 'unknown'); % Loads the neurons from a specific NWB file and a specified region +% spikes_NWB = bz_GetSpikes_NWB('nwb_file', nwb_file, 'saveMat',true); % saves a cellInfo.mat file with all the Neurons + + +% Konstantinos Nasiotis 2019 + + + +%% TODO +disp(' 1. The template waveform should be filled by: nwb2.units.waveform_mean The example dataset didnt have any. I added a template') +disp(' 2. The maxWaveformCh should be added on the example dataset') +disp(' 3. CluID is probably filled by Kilosort - should be added on the example dataset') +disp(' 4. getWaveforms is not supported. Not sure nwb holds the spiking waveforms - Havent seen a field - UPDATE, THIS IS NOW SUPPORTED - CHECK nwb2.processing.get("ecephys").nwbdatainterface.get("SpikeEventSeries1")') + + + +%% Deal With Inputs +spikeGroupsValidation = @(x) assert(isnumeric(x) || strcmp(x,'all'),... + 'spikeGroups must be numeric or "all"'); + +p = inputParser; +addParameter(p,'nwb_file','',@isstr); +addParameter(p,'spikeGroups','all',spikeGroupsValidation); +addParameter(p,'region','',@isstr); % won't work without sessionInfodata +addParameter(p,'UID',[],@isvector); +addParameter(p,'basepath',pwd,@isstr); +addParameter(p,'getWaveforms',true,@islogical) +addParameter(p,'forceReload',false,@islogical); +addParameter(p,'saveMat',false,@islogical); +addParameter(p,'noPrompts',false,@islogical); + +parse(p,varargin{:}) + +%% Adding support only for spikeGroups, region and UID for now + +nwb_file = p.Results.nwb_file; +spikeGroups = p.Results.spikeGroups; +region = p.Results.region; +UID = p.Results.UID; +% basepath = p.Results.basepath; +% getWaveforms = p.Results.getWaveforms; +% forceReload = p.Results.forceReload; +saveMat = p.Results.saveMat; +% noPrompts = p.Results.noPrompts; + + + + +%% let's check that there is an appropriate NWB file +if isempty(nwb_file) + %disp('No nwb_file given, so we look for a *nwb file in the current directory') + d = dir('*nwb'); + if length(d) > 1 % we assume one .nwb file or this should break + error('there is more than one .nwb file in this directory?'); + elseif length(d) == 0 + d = dir('*nwb'); + if isempty(d) + error('could not find an nwb file..') + end + end + nwb_file = fullfile(d.folder, d.name); +end + + + +nwb2 = nwbRead(nwb_file); +[the_path, ~, ~] = fileparts(nwb_file); + + +% load([new_path_for_files filesep name '.sessionInfo.mat']) +sessionInfo = get_SessionInfoFromNWB(nwb2); + + + nNeurons = length(nwb2.units.id.data.load); + + +% % % % % % % % % % % % % % % % % % % % % % +% CHEK THIS nwb2.units.electrode_group.data(1).refresh(nwb2) +% % % % % % % % % % % % % % % % % % % % + + + % Assign neuron to region based on the region that its Shank belongs to + + % Get a single electrode that belongs in that shank + electrodes2Shank = nwb2.general_extracellular_ephys_electrodes.vectordata.get('group_name').data; + electrodes2Region = nwb2.general_extracellular_ephys_electrodes.vectordata.get('location').data; + + neurons2ShankNames = cell(1,nNeurons); + all_region = cell(1,nNeurons); + + for iNeuron = 1:length(nwb2.units.electrode_group.data) + path = nwb2.units.electrode_group.data(iNeuron).path; + inds = strfind(path, '/'); + neurons2ShankNames{iNeuron} = path(inds(end)+1:end); + + ElectrodesInThatShank = find(strcmp(electrodes2Shank,neurons2ShankNames{iNeuron})); + all_region{iNeuron} = electrodes2Region{ElectrodesInThatShank(1)}; + + end + + + uniqueShanks = unique(neurons2ShankNames,'legacy'); + shankID_of_selected_Neurons = zeros(1, nNeurons); + for iNeuron = 1:nNeurons + shankID_of_selected_Neurons(iNeuron) = find(strcmp(uniqueShanks, neurons2ShankNames{iNeuron})); + end + + + %% Make the check here of what will be loaded + + selected_neurons_UID = false(1,nNeurons); + selected_neurons_spikeGroups = false(1,nNeurons); + selected_neurons_region = false(1,nNeurons); + + % Check which neurons were selected + if ~isempty(UID) + selected_neurons_UID(UID) = true; + end + %Check which Shanks were selected + if ~isempty(spikeGroups) + for iShank = spikeGroups + selected_neurons_spikeGroups(find(shankID_of_selected_Neurons==iShank)) = true; + end + end + %Check which regions were selected + if ~isempty(region) + selected_neurons_region(find(contains(all_region,region))) = true; + end + + UID = find(selected_neurons_UID | selected_neurons_spikeGroups | selected_neurons_region); + + if isempty(UID) + warning(['Warning: The selection made didnt specify any subgroup of neurons. Loading all neurons!']) + UID = 1:nNeurons; + end + + + %% The first time that this function is loaded, make sure to save a .spikes.cellInfo.mat that contains information from all neurons + + if saveMat + UID_forSaveMAt = 1:nNeurons; + temp = get_the_spikes_from_selected_UIDs(UID_forSaveMAt, nwb2, sessionInfo, saveMat, all_region, neurons2ShankNames, the_path); clear temp + end + + spikes = get_the_spikes_from_selected_UIDs(UID, nwb2, sessionInfo, 0, all_region, neurons2ShankNames, the_path); + +end + + + + +function spikes = get_the_spikes_from_selected_UIDs(UID, nwb2, sessionInfo, saveMat, all_region, neurons2ShankNames, the_path) + + %% Get Spikes + spikes = struct; + spikes.samplingRate = sessionInfo.rates.wideband; + spikes.UID = UID; + + % The template waveform should be filled by: nwb2.units.waveform_mean + + template_Waveform = [21.2963012297339,21.2206998551635,21.5609060407305,... % This is from a template I used on Kilosort. + 22.8392565561944,28.2241362812803,30.1107342194246,... % I assign the same on every neuron. + 23.9595314702837,24.3822118826549,23.4681225355758,... % Check if I can get this from nwb + 18.0866792366068,11.9973321575690,-2.96486715514583,... + -110.074832790885,-186.628097395696,-240.085142069235,... + -183.947685024562,-143.562805299476,-104.992358564081,... + -3.29132763624548,22.3478476214865,44.7121087898714,... + 73.1141706455415,79.3615933259538,82.9182943568817,... + 75.0076414359195,70.1691534634109,65.2413184118645,... + 51.9698407486342,45.6296345630672,41.3822118826549,... + 37.1519713328267,31.5334146317958]; + + times = cell(1,length(UID)); + rawWaveform = cell(1,length(UID)); + spindices = []; + + entry = 0; + for iNeuron = UID + entry = entry+1; + if iNeuron == 1 + times_temp = nwb2.units.spike_times.data.load(1:sum(nwb2.units.spike_times_index.data.load(iNeuron))); + else + times_temp = nwb2.units.spike_times.data.load(sum(nwb2.units.spike_times_index.data.load(iNeuron-1))+1:sum(nwb2.units.spike_times_index.data.load(iNeuron))); + end + times{entry} = times_temp(times_temp~=0)'; + + rawWaveform{entry} = template_Waveform; % The template waveform should be filled by: nwb2.units.waveform_mean The example dataset didn't have any + spindices = [spindices ; times{entry} ones(length(times{entry}),1)*iNeuron]; + end + + % Spindices have to be sorted according to when each spike occured + [~,sortedIndices] = sort(spindices(:,1)); + spindices = spindices(sortedIndices,:); + + + uniqueShanks = unique(neurons2ShankNames,'legacy'); + shankID_of_selected_Neurons = zeros(1, length(UID)); + for iNeuron = 1:length(UID) + shankID_of_selected_Neurons(iNeuron) = find(strcmp(uniqueShanks, neurons2ShankNames{UID(iNeuron)})); + end + + + spikes.times = times; + spikes.shankID = shankID_of_selected_Neurons; + spikes.cluID = ones(1,length(UID))*(-1); % THESE ARE THE SPIKING TEMPLATES. THEY ARE FILLED FROM KILOSORT. I ADD A NEGATIVE VALUE TO SEE IF IT CAUSES AN ERROR SOMEWHERE + spikes.rawWaveform = rawWaveform; + spikes.maxWaveformCh = ones(1,length(UID))*(-1); % THESE ASSIGN THE MAXIMUM WAVEFORM TO A ACHANNEL. CHECK HOW TO ADD THIS. I ADD A NEGATIVE VALUE TO SEE IF IT CAUSES AN ERROR SOMEWHERE + spikes.sessionName = sessionInfo.FileName; + spikes.numcells = length(UID); + spikes.spindices = spindices; % This holds the timing of each spike, sorted, and the neuron it belongs to. + + spikes.region = all_region(UID); + + if saveMat + save(fullfile(the_path, [sessionInfo.FileName '.spikes.cellinfo.mat']),'spikes') + end +end + + + + + diff --git a/NWB/bz_LoadBehavior_NWB.m b/NWB/bz_LoadBehavior_NWB.m new file mode 100644 index 00000000..bc48115a --- /dev/null +++ b/NWB/bz_LoadBehavior_NWB.m @@ -0,0 +1,436 @@ +function behavior = bz_LoadBehavior_NWB(varargin) + +%% Load Behavior from NWB files following the Buzcode format +% behavior = bz_LoadBehavior_NWB; % Checks the current directory for a NWB file. Gives a list with the available behavior fields in the NWB file for the user to choose which one to load +% behavior = bz_LoadBehavior_NWB('nwb_file', nwb_file); % Checks the specified NWB file. Gives a list with the available behavior fields in the NWB file for the user to choose which one to load +% behavior = bz_LoadBehavior_NWB('nwb_file', nwb_file, 'EightMazePosition_norm_spatial_series'); % Loads the specific behavior from the specified NWB file +% behavior = bz_LoadBehavior_NWB('EightMazePosition_norm_spatial_series'); % Loads the specific behavior from a NWB file in the current directory + + +% The behaviors are loaded straight from the NWB file. + +% This code creates a file with the behavior selected the first time it's used. +% This is done in order to speed up the process of calling it later +% (.map takes a long time to compute if the trials are long). + + +% This code is based on the bz_LoadBehavior +% Konstantinos Nasiotis 2019 + + +%% TODO + +disp('1. Check if the BAD trials should be included in the computation of the .map field') + + + + +%% + p = inputParser; + addParameter(p,'nwb_file','',@isstr); + + if size(varargin,2)>2 + behaviorName = varargin{3}; + varargin = {varargin{1,1:2}}; + parse(p,varargin{:}) + nwb_file = p.Results.nwb_file; + elseif size(varargin,2)==1 + nwb_file = []; + behaviorName = varargin{1}; + elseif size(varargin,2)==2 + parse(p,varargin{:}) + nwb_file = p.Results.nwb_file; + behaviorName = []; + elseif size(varargin,2)==0 + nwb_file = []; + behaviorName = []; + end + + + %% let's check that there is an appropriate NWB file + if isempty(nwb_file) + %disp('No nwb_file given, so we look for a *nwb file in the current directory') + d = dir('*nwb'); + if length(d) > 1 % we assume one .nwb file or this should break + error('there is more than one .nwb file in this directory?'); + elseif length(d) == 0 + d = dir('*nwb'); + if isempty(d) + error('could not find an NWB file..') + end + end + nwb_file = fullfile(d.folder, d.name); + end + + + nwb2 = nwbRead(nwb_file); + + % Check if behavior fields exists in the dataset + try + allBehaviorKeys = keys(nwb2.processing.get('behavior').nwbdatainterface); + + behavior_exist_here = ~isempty(allBehaviorKeys); + if ~behavior_exist_here + disp('No behavior in this .nwb file') + return + + else + disp(' ') + disp('The following behavior types are present in this dataset') + disp('------------------------------------------------') + for iBehavior = 1:length(allBehaviorKeys) + disp(allBehaviorKeys{iBehavior}) + end + disp(' ') + end + catch + disp('No behavior in this .nwb file') + return + end + + %% Check if the behavior file has already been computed (it takes some time to fill the .mapping field, especially when the trials are long) + % If the behavior file is already computed, just load the file and return + + % FIRST CHECK - CHECK IF THE behaviorName REQUESTED EXISTS + + [the_path, name, ext] = fileparts(nwb_file); + + if exist(fullfile(the_path, [name '.' behaviorName '.behavior.mat']) , 'file') == 2 + load(fullfile(the_path, [name '.' behaviorName '.behavior.mat'])) + disp('Just loading the already computed behavior struct') + return + else + disp('A file with the specified Behavior hasnt already been created. Select the behavior from the list') + end + + %% If a the behavior file doesn't exist, create a new one + + % Check if a specific behavior was called to be loaded. If not, display a + % pop-up list with the available behaviors for selection + [iBehavior, ~] = listdlg('PromptString','Which behavior type would you like to load?',... + 'ListString',allBehaviorKeys,'SelectionMode','single'); + + % Fill the behavior + behavior = struct; % Initialize + + + % Select the key that is within the Behavior selection + allBehaviorSUBKeys = keys(nwb2.processing.get('behavior').nwbdatainterface.get(allBehaviorKeys(iBehavior)).spatialseries); + if length(allBehaviorSUBKeys)>1 + [iSubKeyBehavior, ~] = listdlg('PromptString','Multiple Behavior channel-selections within the selected behavior',... + 'ListString',allBehaviorSUBKeys,'SelectionMode','single'); + else + iSubKeyBehavior = 1; + end + + %% Check if the behavior file has already been computed (it takes some time to fill the .mapping field, especially when the trials are long) + % If the behavior file is already computed, just load the file and return + + % SECOND CHECK - HERE A SUBKEY HAS ALREADY BEEN SET FROM THE SELECTION + + if exist(fullfile(the_path, [name '.' allBehaviorSUBKeys{iSubKeyBehavior} '.behavior.mat']) , 'file') == 2 + load(fullfile(the_path, [name '.' allBehaviorSUBKeys{iSubKeyBehavior} '.behavior.mat'])) + disp('Just loading the already computed behavior struct') + return + end + + + %% Get the behavior info + disp('computing') + + selected_behavior = nwb2.processing.get('behavior').nwbdatainterface.get(allBehaviorKeys(iBehavior)).spatialseries.get(allBehaviorSUBKeys{iSubKeyBehavior}); % This is for easier reference through the code + + behavior.samplingRate = 1000 ./ nanmean(diff(selected_behavior.timestamps.load))./1000; + behavior.units = selected_behavior.data_unit; + + Info = allBehaviorKeys{iBehavior}; +% behavior.description = selected_behavior.description; + + % Check if timestamps and data have values. If not something weird + % is going on + if ~isempty(selected_behavior.timestamps) + behavior.timestamps = selected_behavior.timestamps.load; + else + behavior.timestamps = []; + warning(['Behavior: ' Info ' --- Timestamps are empty: weird']) + end + if ~isempty(selected_behavior.data) + + % Specifically for position signals, seperate to x and y + if ~isempty(strfind(allBehaviorKeys(iBehavior),'position')) + data = selected_behavior.data.load; + datasubstruct.x = data(1,:)'; + if size(data,1)>1 + datasubstruct.y = data(2,:)'; + end + if size(data,1)>2 + datasubstruct.z = data(3,:)'; + end + else + datasubstruct.data = selected_behavior.data.load; + + end + + + else + datasubstruct.data = []; + warning(['Behavior: ' Info ' --- Data is empty: weird']) + end + + + % position: .x, .y, and .z + % units: millimeters, centimeters, meters[default] + % orientation: .x, .y, .z, and .w + % rotationType: euler or quaternion[default] + % pupil: .x, .y, .diameter + + +% datasubstruct.reference_frame = nwb2.acquisition.get(all_raw_keys(allBehaviorKeys(iBehavior))).reference_frame; This seems to be present only on the position_sensor channels + + % Conside removing these + datasubstruct.starting_time_unit = selected_behavior.starting_time_unit; + datasubstruct.timestamps_interval = selected_behavior.timestamps_interval; + datasubstruct.comments = selected_behavior.comments; + datasubstruct.control = selected_behavior.control; + datasubstruct.control_description = selected_behavior.control_description; + datasubstruct.data_resolution = selected_behavior.data_resolution; + datasubstruct.starting_time = selected_behavior.starting_time; + datasubstruct.help = selected_behavior.help; + + + + %% Exclude special characters from the behavior name + BehaviorLabel = allBehaviorKeys{iBehavior}; + + clean_string = regexprep(BehaviorLabel,'[^a-zA-Z0-9_]',''); + if ~strcmp(clean_string, BehaviorLabel) + disp(['The variable name (' BehaviorLabel ') of the Behavior was changed to exclude special characters']) + BehaviorLabel = clean_string; + end + + %% If this is a position signal, rename to position to have compatibility with the Buzsaki functions + if ~isempty(strfind(allBehaviorKeys(iBehavior),'position')) + BehaviorLabel = 'position'; + end + + + behavior.(BehaviorLabel) = datasubstruct; + + % BehaviorInfo + behaviorinfo.description = selected_behavior.description; + behaviorinfo.acquisitionsystem = 'Fill Me'; + behaviorinfo.substructnames = {BehaviorLabel}; + behavior.behaviorinfo = behaviorinfo; + + + + %% Fill the events substructure + + nTrials = length(nwb2.intervals_trials.start_time.data.load); + uniqueConditions = unique(nwb2.intervals_trials.vectordata.get('condition').data); + + + all_conditions = nwb2.intervals_trials.vectordata.get('condition').data; + for iCondition = 1:length(all_conditions) + + events.trialConditions(iCondition) = find(strcmp(all_conditions{iCondition}, uniqueConditions));%1x221 double which unique condition each trialCondition belongs to - INDEX + + end + events.trialIntervals = [nwb2.intervals_trials.start_time.data.load nwb2.intervals_trials.stop_time.data.load]; % 221x2 double start-end of trial in seconds + events.conditionType = unique(nwb2.intervals_trials.vectordata.get('condition').data)'; % 1x10 cell (central, wheel ...) + +% nwb2.intervals_trials.id.data.load +% nwb2.intervals_trials.colnames +% nwb2.intervals_trials.vectordata +% nwb2.intervals_trials.vectordata.get('both_visit').data.load +% nwb2.intervals_trials.vectordata.get('condition').data +% nwb2.intervals_trials.vectordata.get('error_run').data.load +% nwb2.intervals_trials.vectordata.get('stim_run').data.load +% nwb2.intervals_trials.vectorindex + + trialTimestampBounds = [nwb2.intervals_trials.start_time.data.load nwb2.intervals_trials.stop_time.data.load]; + position_timestamps = selected_behavior.timestamps.load; + + for iTrial = 1:nTrials + + %% Load only the samples from each trial + % Find the indices of the timestamps that are selected + + [~, iPositionTimestamps] = histc(trialTimestampBounds(iTrial,:), position_timestamps); + + if iPositionTimestamps(1)~=0 && iPositionTimestamps(2)==0 + iPositionTimestamps(2) = length(position_timestamps); + end + + try + position = selected_behavior.data.load([1, iPositionTimestamps(1)], [Inf, iPositionTimestamps(2)]); + catch + error ('Error loading selected behavior file. Are there datapoints within the trials selection time-period?') + end + + % The boundaries struct will hold the position boundaries of each + % TYPE of CONDITION, for all the trials on each condition Type + if ~exist('boundaries') + if size(position,1)==1 + boundaries = repmat(struct('x_min',0,'x_max',0),length(uniqueConditions),1); + elseif size(position,1)==2 + boundaries = repmat(struct('x_min',0,'x_max',0,'y_min',0,'y_max',0),length(uniqueConditions),1); + elseif size(position,1)==3 + boundaries = repmat(struct('x_min',0,'x_max',0,'y_min',0,'y_max',0,'z_min',0,'z_max',0),length(uniqueConditions),1); + end + end + + +% clr = rand(1,3); +% plot(position(1,:),position(2,:),'.','color',clr); +% drawnow + + % Get the index of the condition that the trials belong to - this will + % be used for defining the spatial boundaries of the condition + iConditionOfTrial = find(strcmp(uniqueConditions, nwb2.intervals_trials.vectordata.get('condition').data(iTrial)))'; + + + trials{1,iTrial}.x = position(1,:)';% 608 x 1 + if size(position,1)>1 + trials{1,iTrial}.y = position(2,:)';% 608 x 1 + end + if size(position,1)>2 + trials{1,iTrial}.z = position(3,:)';% 608 x 1 + end + + trials{1,iTrial}.timestamps = position_timestamps(iPositionTimestamps(1):iPositionTimestamps(2));% 608 x 1 + trials{1,iTrial}.direction = 'FILL ME'; % 'clockwise' 'counterclockwise' + trials{1,iTrial}.type = nwb2.intervals_trials.vectordata.get('condition').data{iTrial}; % 'central alternation' + % trials{1,iTrial}.orientation = % 1x1 struct + % trials{1,iTrial}.errorPerMarker = 608 x 1 + + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Map the trials to a 1x201 vector + bins = 200; + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % normalize positions to template + c=1; + + uniqueConditions = unique((events.trialConditions)); + + for iCondition = 1:length(uniqueConditions) + + map{iCondition}=[]; + t_conc=[]; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % MAYBE SELECT ONLY THE TRIALS THAT WEREN'T BAD HERE +% iTrialsInCondition = find(events.trialConditions == uniqueConditions(iCondition) & ~nwb2.intervals_trials.vectordata.get('error_run').data.load'); + iTrialsInCondition = find(events.trialConditions == uniqueConditions(iCondition)); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + + for iTrial = 1:length(iTrialsInCondition) + selected_trial = trials{iTrialsInCondition(iTrial)}; % I set this for faster typing + + if size(position,1)==1 + t_conc = [selected_trial.timestamps, selected_trial.x, 20*(selected_trial.timestamps - selected_trial.timestamps(1))]; % Check what this 20 is + elseif size(position,1)==2 + t_conc = [selected_trial.timestamps, selected_trial.x selected_trial.y, 20*(selected_trial.timestamps - selected_trial.timestamps(1))]; % Check what this 20 is + elseif size(position,1)==3 + t_conc = [selected_trial.timestamps, selected_trial.x selected_trial.y selected_trial.z, 20*(selected_trial.timestamps - selected_trial.timestamps(1))]; % Check what this 20 is + end + +% % % % % % templateVector = 1:nBins; +% % % % % % target = size(position,2); +% % % % % % n_ent = length(templateVector); +% % % % % % trials{1,iTrial}.mapping = round(interp1( 1:n_ent, templateVector, linspace(1, n_ent, target) ))'; + + disp('the computation below doesnt make the usage of this function practical for on the fly retrieval of the Behavior') + if length(t_conc)>=bins + while length(t_conc)>bins+1 + + di = pdist(t_conc); + s = squareform(di); + s(find(eye(size(s))))=nan; + [a b] = min(s(:)); + [coords blah] = find(s==a); + t_conc(coords(1),:) = (t_conc(coords(1),:)+t_conc(coords(2),:))./2; + t_conc(coords(2),:) = []; + end + t_conc_all(iTrial,:,:) = t_conc; + end + end + if length(iTrialsInCondition)>0 + map{iCondition} = squeeze(median(t_conc_all(:,:,:),1)); + end + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Assign to the right behavior field + if size(position,1)==1 % 1D + events.map{iCondition}.x = map{iCondition}(:,2); + elseif size(position,1)==2 % 2D + events.map{iCondition}.x = map{iCondition}(:,2); + events.map{iCondition}.y = map{iCondition}(:,3); + elseif size(position,1)==3 % 3D + events.map{iCondition}.x = map{iCondition}(:,2); + events.map{iCondition}.y = map{iCondition}(:,3); + events.map{iCondition}.z = map{iCondition}(:,4); + end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + clear t_conc_all + + disp('finding mapping...') + for iTrial =1:length(iTrialsInCondition) % all trial types (rotations) + selected_trial = trials{iTrialsInCondition(iTrial)}; % I set this for faster typing + for iPoint = 1:length(selected_trial.timestamps) + + if size(position,1)==1 % 1D + [a b] = min(nansum(abs([selected_trial.timestamps(iPoint)-map{iCondition}(:,1),... % Timestamp + selected_trial.x(iPoint)-map{iCondition}(:,2),... % X_COORDINATE + (selected_trial.timestamps(iPoint)-selected_trial.timestamps(1))*50-map{iCondition}(:,1),... % penalty for time differences + (40*(iPoint./length(selected_trial.timestamps)*length(map{iCondition}) - (1:length(map{iCondition})))')])')); % penalty for order differences + + elseif size(position,1)==2 % 2D + [a b] = min(nansum(abs([selected_trial.timestamps(iPoint)-map{iCondition}(:,1),... % Timestamp + selected_trial.x(iPoint)-map{iCondition}(:,2),... % X_COORDINATE + selected_trial.y(iPoint)-map{iCondition}(:,3),... % Y_COORDINATE + (selected_trial.timestamps(iPoint)-selected_trial.timestamps(1))*50-map{iCondition}(:,1),... % penalty for time differences + (40*(iPoint./length(selected_trial.timestamps)*length(map{iCondition}) - (1:length(map{iCondition})))')])')); % penalty for order differences + + elseif size(position,1)==3 % 3D + [a b] = min(nansum(abs([selected_trial.timestamps(iPoint)-map{iCondition}(:,1),... % Timestamp + selected_trial.x(iPoint)-map{iCondition}(:,2),... % X_COORDINATE + selected_trial.y(iPoint)-map{iCondition}(:,3),... % Y_COORDINATE + selected_trial.z(iPoint)-map{iCondition}(:,4),... % Z_COORDINATE + (selected_trial.timestamps(iPoint)-selected_trial.timestamps(1))*50-map{iCondition}(:,1),... % penalty for time differences + (40*(iPoint./length(selected_trial.timestamps)*length(map{iCondition}) - (1:length(map{iCondition})))')])')); % penalty for order differences + end + + mapping{iCondition}{iTrial}(iPoint,:) = [map{iCondition}(b,1:end) b selected_trial.timestamps(iPoint)]; + + end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Assign to the right behavior field + trials{iTrialsInCondition(iTrial)}.mapping = mapping{iCondition}{iTrial}(:,end-1); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + end + end + + events.trials = trials; + behavior.events = events; + + %% Check that the behavior structure meets buzcode standards + [isBehavior] = bz_isBehavior(behavior); + switch isBehavior + case false + warning('Your behavior structure does not meet buzcode standards. Sad.') + end + + %% Save the behavior files if they are not already saved + if exist(fullfile(the_path, [name '.' allBehaviorSUBKeys{iSubKeyBehavior} '.behavior.mat']) ,'file') ~=2 + save(fullfile(the_path, [name '.' allBehaviorSUBKeys{iSubKeyBehavior} '.behavior.mat']), 'behavior') + end + +end \ No newline at end of file diff --git a/NWB/bz_LoadEvents_NWB.m b/NWB/bz_LoadEvents_NWB.m new file mode 100644 index 00000000..76afc640 --- /dev/null +++ b/NWB/bz_LoadEvents_NWB.m @@ -0,0 +1,121 @@ +function events = bz_LoadEvents_NWB(varargin) + +%% Load events from NWB files following the Buzcode format + +% Example calls: +% events = bz_LoadEvents_NWB; % Checks the current directory for a NWB file. Gives a list with the available events in the NWB file for the user to choose which one to load +% events = bz_LoadEvents_NWB('nwb_file', nwb_file); % Checks the specified NWB file. Gives a list with the available events in the NWB file for the user to choose which one to load +% events = bz_LoadEvents_NWB('nwb_file', nwb_file, 'PulseStim_5V_77777ms_LD12'); % Loads the specified event type from the NWB file +% events = bz_LoadEvents_NWB('PulseStim_5V_77777ms_LD12'); % Checks the current directory for a NWB file. Loads the specified file + + +% The events are loaded straight from the NWB file. No extra files needed. + +% This code is based on the bz_LoadEvents +% Konstantinos Nasiotis 2019 + + +%% + p = inputParser; + addParameter(p,'nwb_file','',@isstr); + + if size(varargin,2)>2 + eventsName = varargin{3}; + varargin = {varargin{1,1:2}}; + parse(p,varargin{:}) + nwb_file = p.Results.nwb_file; + elseif size(varargin,2)==1 + nwb_file = []; + eventsName = varargin{1}; + elseif size(varargin,2)==2 + parse(p,varargin{:}) + nwb_file = p.Results.nwb_file; + elseif size(varargin,2)==0 + nwb_file = []; + end + + %% let's check that there is an appropriate NWB file + if isempty(nwb_file) + %disp('No nwb_file given, so we look for a *nwb file in the current directory') + d = dir('*nwb'); + if length(d) > 1 % we assume one .nwb file or this should break + error('there is more than one .nwb file in this directory?'); + elseif length(d) == 0 + d = dir('*nwb'); + if isempty(d) + error('could not find an nwb file..') + end + end + nwb_file = fullfile(d.folder, d.name); + end + + nwb2 = nwbRead(nwb_file); + + % Check if an events field exists in the dataset + try + events_exist_here = ~isempty(nwb2.stimulus_presentation); + if ~events_exist_here + disp('No events in this .nwb file') + return + + else + all_event_keys = keys(nwb2.stimulus_presentation); + disp(' ') + disp('The following event types are present in this dataset') + disp('------------------------------------------------') + for iEvent = 1:length(all_event_keys) + disp(all_event_keys{iEvent}) + end + disp(' ') + end + catch + disp('No events in this .nwb file') + return + end + + + % Check if a specific event was called to be loaded. If not, display a + % pop-up list with the available events for selection + if ~exist('eventsName','var') + [iEvent, ~] = listdlg('PromptString','Which event type would you like to load?',... + 'ListString',all_event_keys,'SelectionMode','single'); + else + if sum(ismember(all_event_keys, eventsName))>0 + iEvent = find(ismember(all_event_keys, eventsName)); + end + end + + events = struct; + + events.timestamps = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).timestamps.load;% neuroscope compatible matrix with 1-2 columns - [starts stops] (in seconds) + events.detectorinfo.detectorname = 'N/A';% substructure with information about the detection method (fields below) + events.detectorinfo.detectionparms = 'N/A';% parameters used for detection + events.detectorinfo.detectiondate = 'N/A';% date of detection + events.detectorinfo.detectionintervals = 'N/A';% [start stop] pairs of intervals used for detection (optional) +% events.detectorinfo.detectionchannel = ;% channel used for detection (optional) + +% % Optional +% events.amplitudes = 1;% [Nx1 matrix] +% events.frequencies = 1;% [Nx1 matrix] +% events.durations = 1;% [Nx1 matrix] + + + % I add here the rest of the parameters that are saved on the nwb file + events.starting_time_unit = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).starting_time_unit; + events.timestamps_interval = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).timestamps_interval; + events.timestamps_unit = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).timestamps_unit; + events.comments = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).comments; + events.control = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).control; + events.control_description = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).control_description; + events.data = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).data; + events.data_conversion = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).data_conversion; + events.data_resolution = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).data_resolution; + events.data_unit = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).data_unit; + events.description = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).description; + events.starting_time = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).starting_time; + events.starting_time_rate = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).starting_time_rate; + events.help = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).help; + + +end + diff --git a/NWB/get_SessionInfoFromNWB.m b/NWB/get_SessionInfoFromNWB.m new file mode 100644 index 00000000..b5cd4f19 --- /dev/null +++ b/NWB/get_SessionInfoFromNWB.m @@ -0,0 +1,181 @@ +function sessionInfo = get_SessionInfoFromNWB(nwb2) + +%% get_SessionInfoFromNWB populates a variable similar to sessionInfo standards. +% Not all fields are filled (most are commented out). +% Only those that were needed so far in the functions tested. +% More can be filled on demand. + + + + +%% Get the type of the recording (its key will be used to get info from the nwb file) +% There will be 2 logical values: +% RawDataPresent +% LFPDataPresent + +% The keys are not finalized yet from the NWB format so some updating will +% be needed. + + +%% Populate fields + + +% First check if there are raw data present + +try + all_raw_keys = keys(nwb2.acquisition); + + for iKey = 1:length(all_raw_keys) + if ismember(all_raw_keys{iKey}, {'ECoG','bla bla bla'}) %%%%%%%% ADD MORE HERE, DON'T KNOW WHAT THE STANDARD FORMATS ARE + iRawDataKey = iKey; + RawDataPresent = 1; + else + RawDataPresent = 0; + end + end + % nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries.get('all_lfp').data +catch + RawDataPresent = 0; +end + + + + +try + % Check if the data is in LFP format + all_lfp_keys = keys(nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries); + + for iKey = 1:length(all_lfp_keys) + if ismember(all_lfp_keys{iKey}, {'lfp','bla bla bla'}) %%%%%%%% ADD MORE HERE, DON'T KNOW WHAT THE STANDARD FORMATS ARE + iLFPDataKey = iKey; + LFPDataPresent = 1; + break % Once you find the data don't look for other keys/trouble + else + LFPDataPresent = 0; + end + end +catch + LFPDataPresent = 0; +end + + +if ~RawDataPresent && ~LFPDataPresent + error 'There is no data in this .nwb - Maybe check if the Keys are labeled correctly' +end + + + +%% Creare SessionInfo +sessionInfo = struct; + +if RawDataPresent + sessionInfo.nChannels = nwb2.acquisition.get(all_raw_keys{iRawDataKey}).data.dims(2); + sessionInfo.samples_NWB = nwb2.acquisition.get(all_raw_keys{iRawDataKey}).data.dims(1); + sessionInfo.rates.wideband = nwb2.acquisition.get(all_raw_keys{iRawDataKey}).starting_time_rate; + sessionInfo.rates.lfp = 1250; %1250 - DEFAULT - CHECK THIS + sessionInfo.lfpSampleRate = 1250; %1250 - DEFAULT - CHECK THIS % tH lfpSampleRate bypasses + +elseif LFPDataPresent + sessionInfo.nChannels = nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries.get(all_lfp_keys{iLFPDataKey}).data.dims(2); + sessionInfo.samples_NWB = nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries.get(all_lfp_keys{iLFPDataKey}).data.dims(1); + sessionInfo.rates.wideband = nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries.get(all_lfp_keys{iLFPDataKey}).starting_time_rate; + sessionInfo.rates.lfp = nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries.get(all_lfp_keys{iLFPDataKey}).starting_time_rate; %I assign the LFP sampling rate that was already used. Not sure yet if a different value than 1250 Hz will cause problems + sessionInfo.lfpSampleRate = nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries.get(all_lfp_keys{iLFPDataKey}).starting_time_rate; %I assign the LFP sampling rate that was already used. Not sure yet if a different value than 1250 Hz will cause problems + + if sessionInfo.rates.wideband <1250 + warning 'Something weird will happen. Not thoroughly tested yet' + end +end + + + +% Get the ElectrodeGroups +% There is a lot of redundancy on this one + +% electrode_description = nwb2.general_extracellular_ephys_electrodes.vectordata.get('electrode_description').data; +% filtering = nwb2.general_extracellular_ephys_electrodes.vectordata.get('filtering').data; +% group = nwb2.general_extracellular_ephys_electrodes.vectordata.get('group').data; +% group_name = nwb2.general_extracellular_ephys_electrodes.vectordata.get('group_name').data; +% imp = nwb2.general_extracellular_ephys_electrodes.vectordata.get('imp').data.load; +location = nwb2.general_extracellular_ephys_electrodes.vectordata.get('location').data; +% shank = nwb2.general_extracellular_ephys_electrodes.vectordata.get('shank').data.load; +% x = nwb2.general_extracellular_ephys_electrodes.vectordata.get('x').data.load; +% y = nwb2.general_extracellular_ephys_electrodes.vectordata.get('y').data.load; +% z = nwb2.general_extracellular_ephys_electrodes.vectordata.get('z').data.load; + + + + + +% nGroups = sum(unique(shank)>0); % -1 doesn't belong to a Shank Group +nGroups = length(unique(nwb2.general_extracellular_ephys_electrodes.vectordata.get('group_name').data)); +uniqueGroupNames = unique(nwb2.general_extracellular_ephys_electrodes.vectordata.get('group_name').data); +groupNames = nwb2.general_extracellular_ephys_electrodes.vectordata.get('group_name').data; + + +sessionInfo.spikeGroups.nGroups = nGroups; +sessionInfo.spikeGroups.nSamples = ones(1,sessionInfo.spikeGroups.nGroups)*32; % The file I found had 32 here + + +id = nwb2.general_extracellular_ephys_electrodes.vectordata.get('amp_channel').data.load; + + +sessionInfo.spikeGroups.groups = cell(1,sessionInfo.spikeGroups.nGroups); +for iGroup = 1:sessionInfo.spikeGroups.nGroups + sessionInfo.spikeGroups.groups{iGroup} = id(ismember(groupNames, uniqueGroupNames{iGroup}))'; + + % Redundant + for iChannel = 1:length(id(ismember(groupNames, uniqueGroupNames{iGroup}))') + sessionInfo.ElecGp{1,iGroup}.channel{1,iChannel} = sessionInfo.spikeGroups.groups{iGroup}(iChannel); + end + + % Redundant + sessionInfo.SpkGrps(iGroup).Channels = id(ismember(groupNames, uniqueGroupNames{iGroup}))'; + sessionInfo.SpkGrps(iGroup).nSamples = 32; + sessionInfo.SpkGrps(iGroup).PeakSample = 16; + sessionInfo.SpkGrps(iGroup).nFeatures = 3; + + +end + + +% Add Region Info +for iChannel = 1:sessionInfo.nChannels + sessionInfo.region{iChannel} = location{iChannel}; +end + + + + +% Get the rest of the Info +sessionInfo.nBits = 16; % ASSUMING THAT NWB GIVES INT16 PRECISION +sessionInfo.rates.video = 0; +sessionInfo.FileName = nwb2.identifier;%%%%%%%%%%%%%%%% no extension - I DON'T USE THE FILENAME OF THE NWB HERE JUST IN CASE SOMEONE CHANGED IT. +% sessionInfo.SampleTime = 50; % 50 no idea +sessionInfo.nElecGps = nGroups; % 13 +% sessionInfo.ElecGp = []; % 1x13 cell (struct with 1x12 cell inside) +% sessionInfo.HiPassFreq = ;% probably the one from the LFP conversion +% sessionInfo.Date = +% sessionInfo.VoltageRange = % 20 +% sessionInfo.Amplification = % 1000 +% sessionInfo.Offset = 0; +% sessionInfo.AnatGrps = +% sessionInfo.spikeGroups.groups = {1:32}; +% sessionInfo.spikeGroups.nSamples = 1; % I ADDED THIS FOR bz_GetSpikes +sessionInfo.channels = 0:sessionInfo.nChannels-1 ;% 1x128 % starts from 0 THIS IS USED IN THE bz_GetLFP in the end to assign channels to regions +% sessionInfo.lfpChans = +% sessionInfo.thetaChans = +% sessionInfo.region = % cell 1x128 +% sessionInfo.depth = 1; % 3324 - single value +% sessionInfo.ca1 = 116; % 116 - single value +% sessionInfo.ca3 = [] % [] +% sessionInfo.ls = []; +% sessionInfo.animal = 'MONKEY'% string +% sessionInfo.refChan = 112; % single value + +sessionInfo.useRaw = RawDataPresent; % I add this to choose which data to convert to .lfp files + + + +end + diff --git a/NWB/testing_BUZSAKI_NWB_functions_fidelity.m b/NWB/testing_BUZSAKI_NWB_functions_fidelity.m new file mode 100644 index 00000000..ffd3f4e5 --- /dev/null +++ b/NWB/testing_BUZSAKI_NWB_functions_fidelity.m @@ -0,0 +1,130 @@ + + +%% Test BUZSAKI FUNCTIONS and NWB FUNCTIONS + +% For the NWB functions, the nwb file_path should be added as an input +% (no need to have only one .nwb in the same folder) +nwb_file = 'C:\Users\McGill\Documents\GitHub\matnwb\Nas\YutaMouse41\YutaMouse41\YutaMouse41.nwb'; + + +% If the nwb_file is not specified, the functions will search for a unique +% .nwb within the current directory + +% The current Buzcode functions and the NWB functions are differentiated by +% the _NWB suffix. + + +% Konstantinos Nasiotis 2019 + + +%% GET_LFP + +% channel ID 5 (= # 6) +tic +lfp = bz_GetLFP(5); +disp(['Loading a single channel with Buzsaki code takes: ' num2str(toc)]) % 18.0 seconds for a 40369 seconds recording at 1250 Hz + +tic +lfp_NWB = bz_GetLFP_NWB(5, 'nwb_file', nwb_file); +disp(['Loading a single channel with NWB code takes: ' num2str(toc)]) % 4.3 seconds for a 40369 seconds recording at 1250 Hz + +figure(1);plot(lfp.data) +figure(2);plot(lfp_NWB.data) + + +% channel ID 5 and 10 (= #6 and #11) +tic +lfp = bz_GetLFP([5 10]); +disp(['Loading a single channel with Buzsaki code takes: ' num2str(toc)]) % 18 seconds for a 40369 seconds recording at 1250 Hz + +tic +lfp_NWB = bz_GetLFP_NWB([5 10], 'nwb_file', nwb_file); +disp(['Loading a single channel with NWB code takes: ' num2str(toc)]) % 4.3 seconds for a 40369 seconds recording at 1250 Hz + +figure(1);subplot(1,2,1);plot(lfp.data(:,1)); subplot(1,2,2);plot(lfp.data(:,2)) +figure(2);subplot(1,2,1);plot(lfp_NWB.data(:,1)); subplot(1,2,2);plot(lfp_NWB.data(:,2)) + + + +% channel ID 5 (= # 6), from 0 to 120 seconds +lfp = bz_GetLFP(5,'restrict',[0 120]); +lfp_NWB = bz_GetLFP_NWB(5,'restrict',[0 120]); + +figure(1);plot(lfp.data) +figure(2);plot(lfp_NWB.data) + + + +% same, plus from 240.2 to 265.23 seconds +lfp = bz_GetLFP(5,'restrict',[0 120;240.2 265.23]); +lfp_NWB = bz_GetLFP_NWB(5, 'restrict',[0 120;240.2 265.23],'nwb_file', nwb_file); + + +figure(1);subplot(1,2,1); plot(lfp(1).data) ; subplot(1,2,2); plot(lfp(2).data) +figure(2);subplot(1,2,1); plot(lfp_NWB(1).data) ; subplot(1,2,2); plot(lfp_NWB(2).data) + + + +% Downsample recording +lfp = bz_GetLFP (5, 'restrict',[0 120], 'downsample', 2); +lfp_NWB = bz_GetLFP_NWB(5, 'restrict',[0 120], 'downsample', 2); + +figure(1);plot(lfp.data) +figure(2);plot(lfp_NWB.data) + + + +%% GET SPIKES + +% Select ALL Neurons +spikes = bz_GetSpikes('UID',[2 4 7]); +spikes_NWB = bz_GetSpikes_NWB('UID',[2 4 7]); +spikes_NWB = bz_GetSpikes_NWB('nwb_file', nwb_file, 'UID',[2 4 7]); + +% Select Neurons: #2, #4, #7 +spikes = bz_GetSpikes('UID',[2 4 7]); +spikes_NWB = bz_GetSpikes_NWB('nwb_file', nwb_file, 'UID',[2 4 7]); + +% Select Neurons from spikeGroups: #2, #4 +spikes = bz_GetSpikes('spikeGroups', [2,4]); +spikes_NWB = bz_GetSpikes_NWB('nwb_file', nwb_file, 'spikeGroups', [2,4]); + +% Select Neurons from region: 'unknown' +spikes = bz_GetSpikes('region', 'unknown'); +spikes_NWB = bz_GetSpikes_NWB('nwb_file', nwb_file, 'region', 'unknown'); + +% Loads and saves all spiking data into buzcode format .cellinfo. struct +spikes = bz_GetSpikes('saveMat',true); +spikes_NWB = bz_GetSpikes_NWB('nwb_file', nwb_file, 'saveMat',true); + + + +%% Load behavior + +% Test that it works +behavior = bz_LoadBehavior; +behavior = bz_LoadBehavior_NWB; +behavior = bz_LoadBehavior_NWB('nwb_file', nwb_file); +behavior = bz_LoadBehavior_NWB('nwb_file', nwb_file, 'EightMazePosition_norm_spatial_series'); +behavior = bz_LoadBehavior_NWB('EightMazePosition_norm_spatial_series'); + + +%% Load events +basePath = 'C:\Users\McGill\Documents\GitHub\matnwb\Nas\YutaMouse41\YutaMouse41'; +events = bz_LoadEvents(basePath, 'YutaMouse41.PulseStim_0V_10021ms_LD0.events'); + + +events = bz_LoadEvents_NWB; +events = bz_LoadEvents_NWB('nwb_file', nwb_file); +events = bz_LoadEvents_NWB('nwb_file', nwb_file, 'PulseStim_5V_77777ms_LD12'); +events = bz_LoadEvents_NWB('PulseStim_5V_77777ms_LD12'); + + +%% bz_firingMap1D + +[firingMaps] = bz_firingMap1D(spikes_NWB,behavior,4,'savemat',false); +% let's look at some ratemaps... +subplot(3,2,2) +imagesc(squeeze(firingMaps.rateMaps{1}(96,:,:))) +ylabel('trial #') +xlabel('position') From 9635a60323d5ce90c22b12cb432bb3334262a9f8 Mon Sep 17 00:00:00 2001 From: Konstantinos Date: Fri, 26 Apr 2019 10:44:33 -0400 Subject: [PATCH 2/5] Moved NWB into IO --- {NWB => io/NWB}/bz_GetLFP_NWB.m | 0 {NWB => io/NWB}/bz_GetSpikes_NWB.m | 0 {NWB => io/NWB}/bz_LoadBehavior_NWB.m | 0 {NWB => io/NWB}/bz_LoadEvents_NWB.m | 0 {NWB => io/NWB}/get_SessionInfoFromNWB.m | 0 {NWB => io/NWB}/testing_BUZSAKI_NWB_functions_fidelity.m | 0 6 files changed, 0 insertions(+), 0 deletions(-) rename {NWB => io/NWB}/bz_GetLFP_NWB.m (100%) rename {NWB => io/NWB}/bz_GetSpikes_NWB.m (100%) rename {NWB => io/NWB}/bz_LoadBehavior_NWB.m (100%) rename {NWB => io/NWB}/bz_LoadEvents_NWB.m (100%) rename {NWB => io/NWB}/get_SessionInfoFromNWB.m (100%) rename {NWB => io/NWB}/testing_BUZSAKI_NWB_functions_fidelity.m (100%) diff --git a/NWB/bz_GetLFP_NWB.m b/io/NWB/bz_GetLFP_NWB.m similarity index 100% rename from NWB/bz_GetLFP_NWB.m rename to io/NWB/bz_GetLFP_NWB.m diff --git a/NWB/bz_GetSpikes_NWB.m b/io/NWB/bz_GetSpikes_NWB.m similarity index 100% rename from NWB/bz_GetSpikes_NWB.m rename to io/NWB/bz_GetSpikes_NWB.m diff --git a/NWB/bz_LoadBehavior_NWB.m b/io/NWB/bz_LoadBehavior_NWB.m similarity index 100% rename from NWB/bz_LoadBehavior_NWB.m rename to io/NWB/bz_LoadBehavior_NWB.m diff --git a/NWB/bz_LoadEvents_NWB.m b/io/NWB/bz_LoadEvents_NWB.m similarity index 100% rename from NWB/bz_LoadEvents_NWB.m rename to io/NWB/bz_LoadEvents_NWB.m diff --git a/NWB/get_SessionInfoFromNWB.m b/io/NWB/get_SessionInfoFromNWB.m similarity index 100% rename from NWB/get_SessionInfoFromNWB.m rename to io/NWB/get_SessionInfoFromNWB.m diff --git a/NWB/testing_BUZSAKI_NWB_functions_fidelity.m b/io/NWB/testing_BUZSAKI_NWB_functions_fidelity.m similarity index 100% rename from NWB/testing_BUZSAKI_NWB_functions_fidelity.m rename to io/NWB/testing_BUZSAKI_NWB_functions_fidelity.m From 5fbdbad27d45506933d11a83ff0d4d247f2ac103 Mon Sep 17 00:00:00 2001 From: Konstantinos Date: Sun, 28 Apr 2019 01:28:54 -0400 Subject: [PATCH 3/5] Converter from Neuroscope2NWB The converter works in a convenient modular way. Example calls are in YutaMouse41_toNWB.m Added an xml importer since I had trouble with the default from Buzcode --- io/NWB/Neuroscope2NWB.m | 666 +++++++++++++++++++++++++++++++++++++ io/NWB/YutaMouse41_toNWB.m | 41 +++ io/NWB/xml2struct.m | 183 ++++++++++ 3 files changed, 890 insertions(+) create mode 100644 io/NWB/Neuroscope2NWB.m create mode 100644 io/NWB/YutaMouse41_toNWB.m create mode 100644 io/NWB/xml2struct.m diff --git a/io/NWB/Neuroscope2NWB.m b/io/NWB/Neuroscope2NWB.m new file mode 100644 index 00000000..96e16276 --- /dev/null +++ b/io/NWB/Neuroscope2NWB.m @@ -0,0 +1,666 @@ +classdef Neuroscope2NWB + +% This function wraps all electrophysiological, behavioral and analyzed +% data into a single NWB file. + +% It was tested on the YutaMouse41-150903 dataset: +% https://buzsakilab.nyumc.org/datasets/SenzaiY/YutaMouse41/YutaMouse41-150903/ + +% Only the folder needs to be specified. It assumes that all Neuroscope, +% .eeg and .mat files exist within the specified folder. + + +% Konstantinos Nasiotis 2019 + + + methods(Static) + + function xml = GetXMLInfo(folder_path) + %% Enter the folder and run everything from in there (easier paths). + % When the converter is done, go back to the initial folder + + [previous_paths, name] = fileparts(folder_path); + +% current_folder = pwd; +% cd (folder_path) + + %% This uses an .xml importer downloaded from MathWorks - File Exchange + % https://www.mathworks.com/matlabcentral/fileexchange/28518-xml2struct + % The loadxml from the Buzcode repo gave errors + + all_files_in_folder = dir(folder_path); + + iXML = []; + for iFile = 1:length(all_files_in_folder) + if strfind(all_files_in_folder(iFile).name,'.xml') + iXML = [iXML iFile]; + end + end + if isempty(iXML) + error 'There are no .xml files in this folder' + elseif length(iXML)>1 + error 'There is more than one .xml in this folder' + end + + xml = xml2struct([folder_path filesep all_files_in_folder(iXML).name]); + xml = xml.parameters; + xml.folder_path = folder_path; + xml.name = name; + end + + + + function nwb = GeneralInfo(xml) + + + %% General Info + nwb_version = '2.0b'; + + session_start_time = datetime(xml.generalInfo.date.Text, ... + 'Format', 'yyyy-MM-dd''T''HH:mm:ssZZ', ... + 'TimeZone', 'local'); + timestamps_reference_time = datetime(xml.generalInfo.date.Text, ... + 'Format', 'yyyy-MM-dd''T''HH:mm:ssZZ', ... + 'TimeZone', 'local'); + + file_create_date = datetime(datestr(clock), ... + 'Format', 'yyyy-MM-dd''T''HH:mm:ssZZ', ... + 'TimeZone', 'local'); + + + nwb = nwbfile( ... + 'session_description' , 'Mouse in open exploration and theta maze', ... + 'identifier' , xml.name, ... + 'session_start_time' , session_start_time,... + 'file_create_date' , file_create_date,... + 'general_experimenter' , xml.generalInfo.experimenters.Text,... + 'general_session_id' , xml.name,... + 'general_institution' , 'NYU' ,... + 'general_lab' , 'Buzsaki',... + 'subject' , 'YutaMouse',... + 'general_related_publications' , 'DOI:10.1016/j.neuron.2016.12.011',... + 'timestamps_reference_time' , session_start_time); + + nwb.general_subject = types.core.Subject( ... + 'description', 'mouse 5', 'genotype', 'POMC-Cre::Arch', 'age', '9 months', ... + 'sex', 'M', 'subject_id', xml.name, 'species', 'Mus musculus'); + end + + + function nwb = getElectrodeInfo (xml,nwb) + %% Get the electrodes' info + + nShanks = length(xml.spikeDetection.channelGroups.group); + + groups = xml.spikeDetection.channelGroups.group; % Use this for simplicity + + all_shank_channels = cell(nShanks,1); % This will hold the channel numbers that belong in each shank + + % Initialize variables + x = []; + y = []; + z = []; + imp = []; + location = []; + shank = []; + group_name = []; + group_object_view = []; + filtering = []; + shank_channel = []; + amp_channel_id = []; + + device_name = 'implant'; + nwb.general_devices.set(device_name, types.core.Device()); + device_link = types.untyped.SoftLink(['/general/devices/' device_name]); + + for iGroup = 1:nShanks + for iChannel = 1:length(groups{iGroup}.channels.channel) + all_shank_channels{iGroup} = [all_shank_channels{iGroup} str2double(groups{iGroup}.channels.channel{iChannel}.Text)]; + shank_channel = [shank_channel; iChannel-1]; + amp_channel_id = [amp_channel_id; str2double(groups{iGroup}.channels.channel{iChannel}.Text)]; + shank = [shank; iGroup]; + group_name = [group_name; ['shank' num2str(iGroup)]]; + group_object_view = [group_object_view; types.untyped.ObjectView(['/general/extracellular_ephys/' ['shank' num2str(iGroup)]])]; + + if ~isfield(groups{iGroup}.channels.channel{iChannel},'position') + x = [x; NaN]; + y = [y; NaN]; + z = [z; NaN]; + end + if ~isfield(groups{iGroup}.channels.channel{iChannel},'imp') + imp = [imp; NaN]; + end + if ~isfield(groups{iGroup}.channels.channel{iChannel},'location') + location{end+1,1} = 'unknown'; + end + if ~isfield(groups{iGroup}.channels.channel{iChannel},'filtering') + filtering = [filtering; NaN]; + end + end + nwb.general_extracellular_ephys.set(['shank' num2str(iGroup)], ... + types.core.ElectrodeGroup( ... + 'description', ['electrode group for shank' num2str(iGroup)], ... + 'location', 'unknown', ... + 'device', device_link)); + end + + variables = {'x'; 'y'; 'z'; 'imp'; 'location'; 'filtering'; 'group'; 'group_name'; 'shank'; 'shank_channel'; 'amp_channel'}; + + % In order to insert string to a table, they need to be converted to a cell + % first (e.g. location(iElectrode)) + for iElectrode = 1:length(x) + if iElectrode == 1 + tbl = table(x(iElectrode),y(iElectrode),z(iElectrode),imp(iElectrode),{location{iElectrode}},filtering(iElectrode),group_object_view(iElectrode),{group_name(iElectrode,:)},shank(iElectrode),shank_channel(iElectrode),amp_channel_id(iElectrode),... + 'VariableNames', variables); + else + tbl = [tbl; {x(iElectrode),y(iElectrode),z(iElectrode),imp(iElectrode),{location{iElectrode}},filtering(iElectrode),group_object_view(iElectrode),{group_name(iElectrode,:)},shank(iElectrode),shank_channel(iElectrode),amp_channel_id(iElectrode)}]; + end + end + + % add the |DynamicTable| object to the NWB file in + % /general/extracellular_ephys/electrodes + electrode_table = util.table2nwb(tbl, 'metadata about extracellular electrodes'); + nwb.general_extracellular_ephys_electrodes = electrode_table; + + end + + + + + function nwb = getUnitsInfo(xml, nwb) + + %% Add the units info (copied from bz_GetSpikes) + + getWaveforms = 1; % Set this to true if you want to add waveforms on the NWB file + + + spikes.samplingRate = str2double(xml.acquisitionSystem.samplingRate.Text); + + + disp('loading spikes from clu/res/spk files..') + % find res/clu/fet/spk files here + cluFiles = dir([xml.folder_path filesep '*.clu*']); + resFiles = dir([xml.folder_path filesep '*.res*']); + if any(getWaveforms) + spkFiles = dir([xml.folder_path filesep '*.spk*']); + end + + % remove *temp*, *autosave*, and *.clu.str files/directories + tempFiles = zeros(length(cluFiles),1); + for i = 1:length(cluFiles) + dummy = strsplit(cluFiles(i).name, '.'); % Check whether the component after the last dot is a number or not. If not, exclude the file/dir. + if ~isempty(findstr('temp',cluFiles(i).name)) | ~isempty(findstr('autosave',cluFiles(i).name)) | isempty(str2num(dummy{length(dummy)})) | find(contains(dummy, 'clu')) ~= length(dummy)-1 + tempFiles(i) = 1; + end + end + cluFiles(tempFiles==1)=[]; + tempFiles = zeros(length(resFiles),1); + for i = 1:length(resFiles) + if ~isempty(findstr('temp',resFiles(i).name)) | ~isempty(findstr('autosave',resFiles(i).name)) + tempFiles(i) = 1; + end + end + if any(getWaveforms) + resFiles(tempFiles==1)=[]; + tempFiles = zeros(length(spkFiles),1); + for i = 1:length(spkFiles) + if ~isempty(findstr('temp',spkFiles(i).name)) | ~isempty(findstr('autosave',spkFiles(i).name)) + tempFiles(i) = 1; + end + end + spkFiles(tempFiles==1)=[]; + end + + if isempty(cluFiles) + disp('no clu files found...') + spikes = []; + return + end + + + % ensures we load in sequential order (forces compatibility with FMAT + % ordering) + for i = 1:length(cluFiles) + temp = strsplit(cluFiles(i).name,'.'); + shanks(i) = str2num(temp{length(temp)}); + end + [shanks ind] = sort(shanks); + cluFiles = cluFiles(ind); %Bug here if there are any files x.clu.x that are not your desired clus + resFiles = resFiles(ind); + if any(getWaveforms) + spkFiles = spkFiles(ind); + end + + % check if there are matching #'s of files + if length(cluFiles) ~= length(resFiles) && length(cluFiles) ~= length(spkFiles) + error('found an incorrect number of res/clu/spk files...') + end + + % use the .clu files to get spike ID's and generate UID and spikeGroup + % use the .res files to get spike times + count = 1; + + ecephys = types.core.ProcessingModule; + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % This section is copied from the ElectrodesInfo + nShanks = length(xml.spikeDetection.channelGroups.group); + groups = xml.spikeDetection.channelGroups.group; % Use this for simplicity + all_shank_channels = cell(nShanks,1); % This will hold the channel numbers that belong in each shank + shank = []; + group_object_view = []; + + for iGroup = 1:nShanks + % Get all_shank_channls again for iGroup = 1:nShanks + for iChannel = 1:length(groups{iGroup}.channels.channel) + all_shank_channels{iGroup} = [all_shank_channels{iGroup} str2double(groups{iGroup}.channels.channel{iChannel}.Text)]; + shank = [shank iGroup]; + group_object_view = [group_object_view; types.untyped.ObjectView(['/general/extracellular_ephys/' ['shank' num2str(iGroup)]])]; + end + end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + for iShank=1:length(cluFiles) + disp(['working on ' cluFiles(iShank).name]) + + temp = strsplit(cluFiles(iShank).name,'.'); + shankID = str2num(temp{length(temp)}); %shankID is the spikegroup number + clu = load(fullfile(xml.folder_path,cluFiles(iShank).name)); + clu = clu(2:end); % toss the first sample to match res/spk files + res = load(fullfile(xml.folder_path,resFiles(iShank).name)); + spkGrpChans = all_shank_channels{iShank}; + + if any(getWaveforms) && sum(clu)>0 %bug fix if no clusters + nSamples = str2double(xml.spikeDetection.channelGroups.group{iShank}.nSamples.Text); + % load waveforms + chansPerSpikeGrp = length(all_shank_channels{iShank}); + fid = fopen(fullfile(xml.folder_path,spkFiles(iShank).name),'r'); + wav = fread(fid,[1 inf],'int16=>int16'); + try %bug in some spk files... wrong number of samples? + wav = reshape(wav,chansPerSpikeGrp,nSamples,[]); + catch + if strcmp(getWaveforms,'force') + wav = nan(chansPerSpikeGrp,nSamples,length(clu)); + display([spkFiles(iShank).name,' error.']) + else + error(['something is wrong with ',spkFiles(iShank).name,... + ' Use ''getWaveforms'', false to skip waveforms or ',... + '''getWaveforms'', ''force'' to write nans on bad shanks.']) + end + end + wav = permute(wav,[3 1 2]); + + + + + %% Get the DynamicTableRegion field for each shank + + electrodes_field = types.core.DynamicTableRegion('table',types.untyped.ObjectView('/general/extracellular_ephys/electrodes'),'description',['shank' num2str(iShank) ' region'],'data',nwb.general_extracellular_ephys_electrodes.id.data(find(shank == iShank)')); + SpikeEventSeries = types.core.SpikeEventSeries('data', wav, 'electrodes', electrodes_field, 'timestamps', res./ spikes.samplingRate); + + %% This section assigns the spike-waveforms in the .NWB + ecephys.nwbdatainterface.set(['SpikeEventSeries' num2str(iShank)],SpikeEventSeries); + + end + + + cells = unique(clu); + % remove MUA and NOISE clusters... + cells(cells==0) = []; + cells(cells==1) = []; % consider adding MUA as another input argument...? + + + for c = 1:length(cells) + spikes.UID(count) = count; % this only works if all shanks are loaded... how do we optimize this? + ind = find(clu == cells(c)); + spikes.times{count} = res(ind) ./ spikes.samplingRate; + spikes.shankID(count) = shankID; + spikes.cluID(count) = cells(c); + + %Waveforms + if any(getWaveforms) + wvforms = squeeze(mean(wav(ind,:,:)))-mean(mean(mean(wav(ind,:,:)))); % mean subtract to account for slower (theta) trends + if prod(size(wvforms))==length(wvforms)%in single-channel groups wvforms will squeeze too much and will have amplitude on D1 rather than D2 + wvforms = wvforms';%fix here + end + for t = 1:size(wvforms,1) + [a(t) b(t)] = max(abs(wvforms(t,:))); + end + [aa bb] = max(a,[],2); + spikes.rawWaveform{count} = wvforms(bb,:); + spikes.maxWaveformCh(count) = spkGrpChans(bb); % Use this in Brainstorm + % %Regions (needs waveform peak) + % if isfield(xml,'region') %if there is regions field in your metadata + % spikes.region{count} = 'unknown'; + % elseif isfield(xml,'units') %if no regions, but unit region from xml via Loadparamteres + % %Find the xml Unit that matches group/cluster + % unitnum = cellfun(@(X,Y) X==spikes.shankID(count) && Y==spikes.cluID(count),... + % {sessionInfo.Units(:).spikegroup},{sessionInfo.Units(:).cluster}); + % if sum(unitnum) == 0 + % display(['xml Missing Unit - spikegroup: ',... + % num2str(spikes.shankID(count)),' cluster: ',... + % num2str(spikes.cluID(count))]) + % spikes.region{count} = 'missingxml'; + % else %possible future bug: two xml units with same group/clu... + % spikes.region{count} = sessionInfo.Units(unitnum).structure; + % end + % end + clear a aa b bb + end + + count = count + 1; + + end + + ecephys.description = 'intermediate data from extracellular electrophysiology recordings, e.g., LFP'; + nwb.processing.set('ecephys', ecephys); + end + + + % Serialize spiketimes and cluIDs + spike_times = []; + spike_times_index = []; + + current_index = 0; + for iNeuron = 1:length(spikes.UID) + spike_times = [spike_times ; spikes.times{iNeuron}]; + spike_times_index = [spike_times_index; int64(length(spikes.times{iNeuron})+current_index)]; + current_index = spike_times_index(end); + end + + + % electrode_group - Assigns the group_object_view that was defined above at + % the electrodes, to specific neurons - I need to find how each neuron is + % assigned to a shank + electrode_group = []; + shank_that_neurons_belongs_to = zeros(length(spikes.UID),1); + for iNeuron = 1:length(spikes.UID) + shank_that_neurons_belongs_to(iNeuron) = str2double(xml.units.unit{iNeuron}.group.Text); + first_electrode_in_shank = find(shank == shank_that_neurons_belongs_to(iNeuron)); + first_electrode_in_shank = first_electrode_in_shank(1); + electrode_group = [electrode_group; group_object_view(first_electrode_in_shank)]; + end + + electrode_group = types.core.VectorData('data', electrode_group, 'description','the electrode group that each spike unit came from'); + + % Initialize the fields needed + spike_times = types.core.VectorData ('data', spike_times, 'description', 'the spike times for each unit'); + spike_times_index = types.core.VectorIndex ('data', spike_times_index, 'target', types.untyped.ObjectView('/units/spike_times')); % The ObjectView links the indices to the spike times + id = types.core.ElementIdentifiers('data', [0:length(xml.units.unit)-1]'); + + + % FOR THE VECTORDATA I NEED FILE: DG_all_6__UnitFeatureSummary_add (PROBABLY - ACCORDING TO THE CONVERTER) + + + % First, instantiate the table, listing all of the columns that will be + % added and us the |'id'| argument to indicate the number of rows. Ifa + % value is indexed, only the column name is included, not the index. For + % instance, |'spike_times_index'| is not added to the array. + nwb.units = types.core.Units( ... + 'electrode_group', electrode_group, 'electrodes', [], 'electrodes_index', [], 'obs_intervals', [], 'obs_intervals_index', [], ... + 'spike_times', spike_times, 'spike_times_index', spike_times_index, 'waveform_mean', [], 'waveform_sd', [], ... + 'colnames', {'shank_id'; 'spike_times'; 'electrode_group'; 'cell_type'; 'global_id'; 'max_electrode'}, ... + 'description', 'Generated from Neuroscope2NWB', 'id', id, 'vectordata', [], 'vectorindex', []); + + end + + + + function nwb = getEvents(xml,nwb) + %% Add events: nwb2.stimulus_presentation + + eventFiles = dir([xml.folder_path filesep '*.evt']); + + for iFile = 1:length(eventFiles) + + if ~strcmp(eventFiles(iFile).name,'YutaMouse41-150903.DS1.ch0.evt') && ~strcmp(eventFiles(iFile).name,'YutaMouse41-150903.DS2.ch0.evt') % Ignore those two files + + events = LoadEvents(eventFiles(iFile).name); % Load Events is a function from the Buzcode - THIS DOESN'T LOAD ANYTHING FOR 'PulseStim_0V_10021ms_LD0' - maybe because it has a single entrty??? + + if ~isempty(events.time) && ~isempty(events.description) + AnnotationSeries = types.core.AnnotationSeries('data',events.description,'timestamps',events.time); + nwb.stimulus_presentation.set(events.description{1}, AnnotationSeries); + end + end + + end + end + + + + function nwb = getBehavior(xml,nwb) + + %% Add behavioral data: nwb2.processing.get('behavior').nwbdatainterface + + behavior = types.core.ProcessingModule; + + behavioralFiles = dir([xml.folder_path filesep '*.position']); + + time_epochs = repmat(struct('label','','start_times',0,'stop_times',0),length(behavioralFiles),1); + + for iFile = 1:length(behavioralFiles) + + % The label of the behavior + behavioral_Label = strsplit(behavioralFiles(iFile).name,'__'); + behavioral_Label = erase(behavioral_Label{2},'.mat'); + + position_signals = load(behavioralFiles(iFile).name); + + % Some behavioral signals might have more than one signal in them + field_names = fieldnames(position_signals); + + the_position_field_NWB = types.core.Position; + + for iField = 1:length(field_names) + + behavioral_timestamps = position_signals.(field_names{iField})(:,1); + position_coordinates = position_signals.(field_names{iField})(:,2:end); + spatial_series = types.core.SpatialSeries('data', position_coordinates, 'timestamps', behavioral_timestamps, 'reference_frame', 'unknown', 'data_conversion', 1); + + the_position_field_NWB.spatialseries.set(field_names{iField}, spatial_series); + end + + behavior.nwbdatainterface.set(behavioral_Label,the_position_field_NWB); + + time_epochs(iFile).start_times = behavioral_timestamps(1,1); + time_epochs(iFile).stop_times = behavioral_timestamps(end,1); + time_epochs(iFile).label = behavioral_Label; + + end + + behavior.description = 'Behavioral signals'; + nwb.processing.set('behavior', behavior); + + end + + + + function nwb = getElectrophysiology(xml,nwb) + + lfpFile = dir([xml.folder_path filesep '*.eeg']); + + if length(lfpFile)>1 + error('More than one .eeg files are present here. Weird') + end + + + % Get the samples number, based on the size of the file + % Check for the precision that samples are saved first + + hdr.nBits = str2double(xml.acquisitionSystem.nBits.Text); + hdr.nChannels = str2double(xml.acquisitionSystem.nChannels.Text); + hdr.sRateOrig = str2double(xml.acquisitionSystem.samplingRate.Text); + hdr.Gain = str2double(xml.acquisitionSystem.amplification.Text); + hdr.sRateLfp = str2double(xml.fieldPotentials.lfpSamplingRate.Text); + + + % Get data type + switch lower(hdr.nBits) + case 16; + hdr.byteSize = 2; + hdr.byteFormat = 'int16'; + case 32; + hdr.byteSize = 4; + hdr.byteFormat = 'int32'; + end + % Guess the number of time points based on the file size + dirInfo = dir(lfpFile.name); + hdr.nSamples = floor(dirInfo.bytes ./ (hdr.nChannels * hdr.byteSize)); + + + + + + lfp_data = bz_LoadBinary(lfpFile.name, 'duration',Inf, 'frequency',hdr.sRateLfp,'nchannels',hdr.nChannels, 'channels', [1:hdr.nChannels]); % nSamples x 64 + + + % If the electrode Information has not already bwwn filled, + % do it now + if isempty(nwb.general_extracellular_ephys_electrodes) + nwb = Neuroscope2NWB.getElectrodeInfo(xml,nwb); + end + + + electrodes_field = types.core.DynamicTableRegion('table',types.untyped.ObjectView('/general/extracellular_ephys/electrodes'),'description','electrode table reference','data',nwb.general_extracellular_ephys_electrodes.id.data); + + lfp = types.core.ElectricalSeries('data', lfp_data', 'electrodes',electrodes_field, 'description', 'lfp signal for all shank electrodes', 'starting_time', 0, 'starting_time_rate', hdr.sRateLfp); + % I TRANSPOSED THE MATRIX HERE TO BE COMPATIBLE WITH THE FUNCTION BZ_GET_LFP + + LFP = types.core.LFP; + LFP.electricalseries.set('lfp',lfp); + + % Check if the ecephys field is already created + if ~ismember(keys(nwb.processing),'ecephys') + ecephys = types.core.ProcessingModule('description', ''); + ecephys.description = 'intermediate data from extracellular electrophysiology recordings, e.g., LFP'; + nwb.processing.set('ecephys', ecephys); + end + + nwb.processing.get('ecephys').nwbdatainterface.set('LFP', LFP); + + end + + + function nwb = getEpochs(xml,nwb) + + + %% Add Epochs + % The epochs are based on the separate behavioral files + % Each separate behavioral file = 1 epoch + + behavioralFiles = dir([xml.folder_path filesep '*.position']); + + time_epochs = repmat(struct('label','','start_times',0,'stop_times',0),length(behavioralFiles),1); + + for iFile = 1:length(behavioralFiles) + + for iField = 1:length(field_names) + behavioral_timestamps = position_signals.(field_names{iField})(:,1); + time_epochs(iFile).start_times = behavioral_timestamps(1,1); + time_epochs(iFile).stop_times = behavioral_timestamps(end,1); + time_epochs(iFile).label = behavioral_Label; + end + + end + + id_epochs = types.core.ElementIdentifiers('data',int64(0:length(time_epochs)-1)'); + + start_time_epochs = types.core.VectorData('description','Starting timepoint of Each Epoch','data',[time_epochs.start_times]'); + stop_time_epochs = types.core.VectorData('description','Ending timepoint of Each Epoch','data',[time_epochs.stop_times]'); + + + vectordata_epochs_label = types.core.VectorData('description','Label of Epoch','data',{time_epochs.label}'); + vectordata_epochs = types.untyped.Set('label', vectordata_epochs_label); + + + %%%%%% THE VECTORDATA IS IGNORED IN HERE - THE LABELS ARE NOT SAVED - + %%%%%% CHECK AGAIN + %%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + intervals_epochs = types.core.TimeIntervals('start_time',start_time_epochs,'stop_time',stop_time_epochs,... + 'colnames',{'start_time';'stop_time';'label'},... + 'description','experimental epochs','id',id_epochs, 'vectordata',vectordata_epochs); + + nwb.intervals_epochs = intervals_epochs; + + + end + + + + + function nwb = getTrials(xml,nwb) + + + %% Add Trials + + % This file holds a matrix with the trial info + trialsFile = dir([xml.folder_path filesep '*Run.mat']); + trialsInfo = load(trialsFile.name); + the_field = fieldnames(trialsInfo); + trialsInfo = trialsInfo.(the_field{1}); + + % This file holds a cell array with the labels for the matrix above (...) + runFile = dir([xml.folder_path filesep '*RunInfo.mat']); + runInfo = load(runFile.name); + the_field = fieldnames(runInfo); + runInfo = runInfo.(the_field{1}); + + + + start_time_trials = types.core.VectorData('description','Starting timepoint of Each Trial','data',trialsInfo(:,1)); + stop_time_trials = types.core.VectorData('description','Ending timepoint of Each Trial','data',trialsInfo(:,2)); + + id_trials = types.core.ElementIdentifiers('data',int64(0:size(trialsInfo,1)-1)'); + + conditions_trials = cell(size(trialsInfo,1),1); + for iTrial = 1:size(trialsInfo,1) + conditions_trials{iTrial} = runInfo{find(trialsInfo(iTrial,3:4))+2}; + end + + vectordata_trials = types.untyped.Set('both_visit', trialsInfo(:,7), 'condition', conditions_trials, 'error_run', trialsInfo(:,5), 'stim_run',trialsInfo(:,6)); + + intervals_trials = types.core.TimeIntervals('start_time',start_time_trials,'stop_time',stop_time_trials,... + 'colnames',{'start_time';'stop_time';'error_run';'stim_run';'both_visit';'condition'},... + 'description','experimental trials','id',id_trials, 'vectordata',vectordata_trials); + + nwb.intervals_trials = intervals_trials; + + end + + + + function nwb = special_YutaMouse_recordings(xml,nwb) + + %% Add raw recordings + + % value taken from Yuta's spreadsheet + + % HOW ABOUT POSITION0 - POSITION1 CHANNELS??? + + hdr.nChannels = str2double(xml.acquisitionSystem.nChannels.Text); + hdr.sRateLfp = str2double(xml.fieldPotentials.lfpSamplingRate.Text); + + lfpFile = dir([xml.folder_path filesep '*.eeg']); + + special_electrode_labels = {'ch_wait','ch_arm','ch_solL','ch_solR','ch_dig1','ch_dig2','ch_entL','ch_entR','ch_SsolL','ch_SsolR'}; + special_electrode_indices = [79,78,76,77,65,68,72,71,73,70]; + + acquisition = types.untyped.Set; + for iSpecialElectrode = 1:length(special_electrode_labels) + + special_Electrode_data = bz_LoadBinary(lfpFile.name, 'duration',Inf, 'frequency',hdr.sRateLfp,'nchannels',hdr.nChannels, 'channels', special_electrode_indices(iSpecialElectrode)); + single_Electrode = types.core.TimeSeries('description','environmental electrode recorded inline with neural data','data',special_Electrode_data,'starting_time', 0, 'starting_time_rate', hdr.sRateLfp, 'data_unit','V'); + acquisition.set(special_electrode_labels{iSpecialElectrode}, single_Electrode); + end + + nwb.acquisition = acquisition; + + end + + + end + +end diff --git a/io/NWB/YutaMouse41_toNWB.m b/io/NWB/YutaMouse41_toNWB.m new file mode 100644 index 00000000..0642da23 --- /dev/null +++ b/io/NWB/YutaMouse41_toNWB.m @@ -0,0 +1,41 @@ +%% Example conversion to NWB of the YutaMouse41-150903 dataset + +% Folder that contains all the files - This is the only input needed +folder_path = 'F:\NWBtoBuzcode\YutaMouse41-150903'; + +%% Start Adding fields to NWB + +% Get info from the xml file +xml = Neuroscope2NWB.GetXMLInfo(folder_path); + +% Add general info to the NWB file +nwb = Neuroscope2NWB.GeneralInfo(xml); + +% Add electrode info +nwb = Neuroscope2NWB.getElectrodeInfo(xml, nwb); + +% Add units info - By default, the spike waveforms are added to the file +nwb = Neuroscope2NWB.getUnitsInfo(xml,nwb); + +% Add stimulation events +nwb = Neuroscope2NWB.getEvents(xml,nwb); + +% Add behavioral info/channels +nwb = Neuroscope2NWB.getBehavior(xml,nwb); + +% Add electrophysiological channels +nwb = Neuroscope2NWB.getElectrophysiology(xml, nwb); + +% Add epochs +nwb = Neuroscope2NWB.getEpochs(xml,nwb); + +% Add trials +nwb = Neuroscope2NWB.getTrials(xml,nwb); + +% Add channels based on the Yuta spreadsheet +nwb = Neuroscope2NWB.special_YutaMouse_recordings(xml,nwb); + + +%% Export to nwb +nwbExport(nwb, 'YutaMouse41_converted.nwb') + diff --git a/io/NWB/xml2struct.m b/io/NWB/xml2struct.m new file mode 100644 index 00000000..bd0412d6 --- /dev/null +++ b/io/NWB/xml2struct.m @@ -0,0 +1,183 @@ +function [ s ] = xml2struct( file ) +%Convert xml file into a MATLAB structure +% [ s ] = xml2struct( file ) +% +% A file containing: +% +% Some text +% Some more text +% Even more text +% +% +% Will produce: +% s.XMLname.Attributes.attrib1 = "Some value"; +% s.XMLname.Element.Text = "Some text"; +% s.XMLname.DifferentElement{1}.Attributes.attrib2 = "2"; +% s.XMLname.DifferentElement{1}.Text = "Some more text"; +% s.XMLname.DifferentElement{2}.Attributes.attrib3 = "2"; +% s.XMLname.DifferentElement{2}.Attributes.attrib4 = "1"; +% s.XMLname.DifferentElement{2}.Text = "Even more text"; +% +% Please note that the following characters are substituted +% '-' by '_dash_', ':' by '_colon_' and '.' by '_dot_' +% +% Written by W. Falkena, ASTI, TUDelft, 21-08-2010 +% Attribute parsing speed increased by 40% by A. Wanner, 14-6-2011 +% Added CDATA support by I. Smirnov, 20-3-2012 +% +% Modified by X. Mo, University of Wisconsin, 12-5-2012 + + if (nargin < 1) + clc; + help xml2struct + return + end + + if isa(file, 'org.apache.xerces.dom.DeferredDocumentImpl') || isa(file, 'org.apache.xerces.dom.DeferredElementImpl') + % input is a java xml object + xDoc = file; + else + %check for existance + if (exist(file,'file') == 0) + %Perhaps the xml extension was omitted from the file name. Add the + %extension and try again. + if (isempty(strfind(file,'.xml'))) + file = [file '.xml']; + end + + if (exist(file,'file') == 0) + error(['The file ' file ' could not be found']); + end + end + %read the xml file + xDoc = xmlread(file); + end + + %parse xDoc into a MATLAB structure + s = parseChildNodes(xDoc); + +end + +% ----- Subfunction parseChildNodes ----- +function [children,ptext,textflag] = parseChildNodes(theNode) + % Recurse over node children. + children = struct; + ptext = struct; textflag = 'Text'; + if hasChildNodes(theNode) + childNodes = getChildNodes(theNode); + numChildNodes = getLength(childNodes); + + for count = 1:numChildNodes + theChild = item(childNodes,count-1); + [text,name,attr,childs,textflag] = getNodeData(theChild); + + if (~strcmp(name,'#text') && ~strcmp(name,'#comment') && ~strcmp(name,'#cdata_dash_section')) + %XML allows the same elements to be defined multiple times, + %put each in a different cell + if (isfield(children,name)) + if (~iscell(children.(name))) + %put existsing element into cell format + children.(name) = {children.(name)}; + end + index = length(children.(name))+1; + %add new element + children.(name){index} = childs; + if(~isempty(fieldnames(text))) + children.(name){index} = text; + end + if(~isempty(attr)) + children.(name){index}.('Attributes') = attr; + end + else + %add previously unknown (new) element to the structure + children.(name) = childs; + if(~isempty(text) && ~isempty(fieldnames(text))) + children.(name) = text; + end + if(~isempty(attr)) + children.(name).('Attributes') = attr; + end + end + else + ptextflag = 'Text'; + if (strcmp(name, '#cdata_dash_section')) + ptextflag = 'CDATA'; + elseif (strcmp(name, '#comment')) + ptextflag = 'Comment'; + end + + %this is the text in an element (i.e., the parentNode) + if (~isempty(regexprep(text.(textflag),'[\s]*',''))) + if (~isfield(ptext,ptextflag) || isempty(ptext.(ptextflag))) + ptext.(ptextflag) = text.(textflag); + else + %what to do when element data is as follows: + %Text More text + + %put the text in different cells: + % if (~iscell(ptext)) ptext = {ptext}; end + % ptext{length(ptext)+1} = text; + + %just append the text + ptext.(ptextflag) = [ptext.(ptextflag) text.(textflag)]; + end + end + end + + end + end +end + +% ----- Subfunction getNodeData ----- +function [text,name,attr,childs,textflag] = getNodeData(theNode) + % Create structure of node info. + + %make sure name is allowed as structure name + name = toCharArray(getNodeName(theNode))'; + name = strrep(name, '-', '_dash_'); + name = strrep(name, ':', '_colon_'); + name = strrep(name, '.', '_dot_'); + + attr = parseAttributes(theNode); + if (isempty(fieldnames(attr))) + attr = []; + end + + %parse child nodes + [childs,text,textflag] = parseChildNodes(theNode); + + if (isempty(fieldnames(childs)) && isempty(fieldnames(text))) + %get the data of any childless nodes + % faster than if any(strcmp(methods(theNode), 'getData')) + % no need to try-catch (?) + % faster than text = char(getData(theNode)); + text.(textflag) = toCharArray(getTextContent(theNode))'; + end + +end + +% ----- Subfunction parseAttributes ----- +function attributes = parseAttributes(theNode) + % Create attributes structure. + + attributes = struct; + if hasAttributes(theNode) + theAttributes = getAttributes(theNode); + numAttributes = getLength(theAttributes); + + for count = 1:numAttributes + %attrib = item(theAttributes,count-1); + %attr_name = regexprep(char(getName(attrib)),'[-:.]','_'); + %attributes.(attr_name) = char(getValue(attrib)); + + %Suggestion of Adrian Wanner + str = toCharArray(toString(item(theAttributes,count-1)))'; + k = strfind(str,'='); + attr_name = str(1:(k(1)-1)); + attr_name = strrep(attr_name, '-', '_dash_'); + attr_name = strrep(attr_name, ':', '_colon_'); + attr_name = strrep(attr_name, '.', '_dot_'); + attributes.(attr_name) = str((k(1)+2):(end-1)); + end + end +end \ No newline at end of file From 7240875b66ee0d6691a06f76865279f0b5dfb639 Mon Sep 17 00:00:00 2001 From: Konstantinos Date: Mon, 6 May 2019 11:46:41 -0400 Subject: [PATCH 4/5] Removed Neuroscope2NWB Converter from the PR --- io/NWB/Neuroscope2NWB.m | 666 ---------------------------------------- 1 file changed, 666 deletions(-) delete mode 100644 io/NWB/Neuroscope2NWB.m diff --git a/io/NWB/Neuroscope2NWB.m b/io/NWB/Neuroscope2NWB.m deleted file mode 100644 index 96e16276..00000000 --- a/io/NWB/Neuroscope2NWB.m +++ /dev/null @@ -1,666 +0,0 @@ -classdef Neuroscope2NWB - -% This function wraps all electrophysiological, behavioral and analyzed -% data into a single NWB file. - -% It was tested on the YutaMouse41-150903 dataset: -% https://buzsakilab.nyumc.org/datasets/SenzaiY/YutaMouse41/YutaMouse41-150903/ - -% Only the folder needs to be specified. It assumes that all Neuroscope, -% .eeg and .mat files exist within the specified folder. - - -% Konstantinos Nasiotis 2019 - - - methods(Static) - - function xml = GetXMLInfo(folder_path) - %% Enter the folder and run everything from in there (easier paths). - % When the converter is done, go back to the initial folder - - [previous_paths, name] = fileparts(folder_path); - -% current_folder = pwd; -% cd (folder_path) - - %% This uses an .xml importer downloaded from MathWorks - File Exchange - % https://www.mathworks.com/matlabcentral/fileexchange/28518-xml2struct - % The loadxml from the Buzcode repo gave errors - - all_files_in_folder = dir(folder_path); - - iXML = []; - for iFile = 1:length(all_files_in_folder) - if strfind(all_files_in_folder(iFile).name,'.xml') - iXML = [iXML iFile]; - end - end - if isempty(iXML) - error 'There are no .xml files in this folder' - elseif length(iXML)>1 - error 'There is more than one .xml in this folder' - end - - xml = xml2struct([folder_path filesep all_files_in_folder(iXML).name]); - xml = xml.parameters; - xml.folder_path = folder_path; - xml.name = name; - end - - - - function nwb = GeneralInfo(xml) - - - %% General Info - nwb_version = '2.0b'; - - session_start_time = datetime(xml.generalInfo.date.Text, ... - 'Format', 'yyyy-MM-dd''T''HH:mm:ssZZ', ... - 'TimeZone', 'local'); - timestamps_reference_time = datetime(xml.generalInfo.date.Text, ... - 'Format', 'yyyy-MM-dd''T''HH:mm:ssZZ', ... - 'TimeZone', 'local'); - - file_create_date = datetime(datestr(clock), ... - 'Format', 'yyyy-MM-dd''T''HH:mm:ssZZ', ... - 'TimeZone', 'local'); - - - nwb = nwbfile( ... - 'session_description' , 'Mouse in open exploration and theta maze', ... - 'identifier' , xml.name, ... - 'session_start_time' , session_start_time,... - 'file_create_date' , file_create_date,... - 'general_experimenter' , xml.generalInfo.experimenters.Text,... - 'general_session_id' , xml.name,... - 'general_institution' , 'NYU' ,... - 'general_lab' , 'Buzsaki',... - 'subject' , 'YutaMouse',... - 'general_related_publications' , 'DOI:10.1016/j.neuron.2016.12.011',... - 'timestamps_reference_time' , session_start_time); - - nwb.general_subject = types.core.Subject( ... - 'description', 'mouse 5', 'genotype', 'POMC-Cre::Arch', 'age', '9 months', ... - 'sex', 'M', 'subject_id', xml.name, 'species', 'Mus musculus'); - end - - - function nwb = getElectrodeInfo (xml,nwb) - %% Get the electrodes' info - - nShanks = length(xml.spikeDetection.channelGroups.group); - - groups = xml.spikeDetection.channelGroups.group; % Use this for simplicity - - all_shank_channels = cell(nShanks,1); % This will hold the channel numbers that belong in each shank - - % Initialize variables - x = []; - y = []; - z = []; - imp = []; - location = []; - shank = []; - group_name = []; - group_object_view = []; - filtering = []; - shank_channel = []; - amp_channel_id = []; - - device_name = 'implant'; - nwb.general_devices.set(device_name, types.core.Device()); - device_link = types.untyped.SoftLink(['/general/devices/' device_name]); - - for iGroup = 1:nShanks - for iChannel = 1:length(groups{iGroup}.channels.channel) - all_shank_channels{iGroup} = [all_shank_channels{iGroup} str2double(groups{iGroup}.channels.channel{iChannel}.Text)]; - shank_channel = [shank_channel; iChannel-1]; - amp_channel_id = [amp_channel_id; str2double(groups{iGroup}.channels.channel{iChannel}.Text)]; - shank = [shank; iGroup]; - group_name = [group_name; ['shank' num2str(iGroup)]]; - group_object_view = [group_object_view; types.untyped.ObjectView(['/general/extracellular_ephys/' ['shank' num2str(iGroup)]])]; - - if ~isfield(groups{iGroup}.channels.channel{iChannel},'position') - x = [x; NaN]; - y = [y; NaN]; - z = [z; NaN]; - end - if ~isfield(groups{iGroup}.channels.channel{iChannel},'imp') - imp = [imp; NaN]; - end - if ~isfield(groups{iGroup}.channels.channel{iChannel},'location') - location{end+1,1} = 'unknown'; - end - if ~isfield(groups{iGroup}.channels.channel{iChannel},'filtering') - filtering = [filtering; NaN]; - end - end - nwb.general_extracellular_ephys.set(['shank' num2str(iGroup)], ... - types.core.ElectrodeGroup( ... - 'description', ['electrode group for shank' num2str(iGroup)], ... - 'location', 'unknown', ... - 'device', device_link)); - end - - variables = {'x'; 'y'; 'z'; 'imp'; 'location'; 'filtering'; 'group'; 'group_name'; 'shank'; 'shank_channel'; 'amp_channel'}; - - % In order to insert string to a table, they need to be converted to a cell - % first (e.g. location(iElectrode)) - for iElectrode = 1:length(x) - if iElectrode == 1 - tbl = table(x(iElectrode),y(iElectrode),z(iElectrode),imp(iElectrode),{location{iElectrode}},filtering(iElectrode),group_object_view(iElectrode),{group_name(iElectrode,:)},shank(iElectrode),shank_channel(iElectrode),amp_channel_id(iElectrode),... - 'VariableNames', variables); - else - tbl = [tbl; {x(iElectrode),y(iElectrode),z(iElectrode),imp(iElectrode),{location{iElectrode}},filtering(iElectrode),group_object_view(iElectrode),{group_name(iElectrode,:)},shank(iElectrode),shank_channel(iElectrode),amp_channel_id(iElectrode)}]; - end - end - - % add the |DynamicTable| object to the NWB file in - % /general/extracellular_ephys/electrodes - electrode_table = util.table2nwb(tbl, 'metadata about extracellular electrodes'); - nwb.general_extracellular_ephys_electrodes = electrode_table; - - end - - - - - function nwb = getUnitsInfo(xml, nwb) - - %% Add the units info (copied from bz_GetSpikes) - - getWaveforms = 1; % Set this to true if you want to add waveforms on the NWB file - - - spikes.samplingRate = str2double(xml.acquisitionSystem.samplingRate.Text); - - - disp('loading spikes from clu/res/spk files..') - % find res/clu/fet/spk files here - cluFiles = dir([xml.folder_path filesep '*.clu*']); - resFiles = dir([xml.folder_path filesep '*.res*']); - if any(getWaveforms) - spkFiles = dir([xml.folder_path filesep '*.spk*']); - end - - % remove *temp*, *autosave*, and *.clu.str files/directories - tempFiles = zeros(length(cluFiles),1); - for i = 1:length(cluFiles) - dummy = strsplit(cluFiles(i).name, '.'); % Check whether the component after the last dot is a number or not. If not, exclude the file/dir. - if ~isempty(findstr('temp',cluFiles(i).name)) | ~isempty(findstr('autosave',cluFiles(i).name)) | isempty(str2num(dummy{length(dummy)})) | find(contains(dummy, 'clu')) ~= length(dummy)-1 - tempFiles(i) = 1; - end - end - cluFiles(tempFiles==1)=[]; - tempFiles = zeros(length(resFiles),1); - for i = 1:length(resFiles) - if ~isempty(findstr('temp',resFiles(i).name)) | ~isempty(findstr('autosave',resFiles(i).name)) - tempFiles(i) = 1; - end - end - if any(getWaveforms) - resFiles(tempFiles==1)=[]; - tempFiles = zeros(length(spkFiles),1); - for i = 1:length(spkFiles) - if ~isempty(findstr('temp',spkFiles(i).name)) | ~isempty(findstr('autosave',spkFiles(i).name)) - tempFiles(i) = 1; - end - end - spkFiles(tempFiles==1)=[]; - end - - if isempty(cluFiles) - disp('no clu files found...') - spikes = []; - return - end - - - % ensures we load in sequential order (forces compatibility with FMAT - % ordering) - for i = 1:length(cluFiles) - temp = strsplit(cluFiles(i).name,'.'); - shanks(i) = str2num(temp{length(temp)}); - end - [shanks ind] = sort(shanks); - cluFiles = cluFiles(ind); %Bug here if there are any files x.clu.x that are not your desired clus - resFiles = resFiles(ind); - if any(getWaveforms) - spkFiles = spkFiles(ind); - end - - % check if there are matching #'s of files - if length(cluFiles) ~= length(resFiles) && length(cluFiles) ~= length(spkFiles) - error('found an incorrect number of res/clu/spk files...') - end - - % use the .clu files to get spike ID's and generate UID and spikeGroup - % use the .res files to get spike times - count = 1; - - ecephys = types.core.ProcessingModule; - - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % This section is copied from the ElectrodesInfo - nShanks = length(xml.spikeDetection.channelGroups.group); - groups = xml.spikeDetection.channelGroups.group; % Use this for simplicity - all_shank_channels = cell(nShanks,1); % This will hold the channel numbers that belong in each shank - shank = []; - group_object_view = []; - - for iGroup = 1:nShanks - % Get all_shank_channls again for iGroup = 1:nShanks - for iChannel = 1:length(groups{iGroup}.channels.channel) - all_shank_channels{iGroup} = [all_shank_channels{iGroup} str2double(groups{iGroup}.channels.channel{iChannel}.Text)]; - shank = [shank iGroup]; - group_object_view = [group_object_view; types.untyped.ObjectView(['/general/extracellular_ephys/' ['shank' num2str(iGroup)]])]; - end - end - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - - for iShank=1:length(cluFiles) - disp(['working on ' cluFiles(iShank).name]) - - temp = strsplit(cluFiles(iShank).name,'.'); - shankID = str2num(temp{length(temp)}); %shankID is the spikegroup number - clu = load(fullfile(xml.folder_path,cluFiles(iShank).name)); - clu = clu(2:end); % toss the first sample to match res/spk files - res = load(fullfile(xml.folder_path,resFiles(iShank).name)); - spkGrpChans = all_shank_channels{iShank}; - - if any(getWaveforms) && sum(clu)>0 %bug fix if no clusters - nSamples = str2double(xml.spikeDetection.channelGroups.group{iShank}.nSamples.Text); - % load waveforms - chansPerSpikeGrp = length(all_shank_channels{iShank}); - fid = fopen(fullfile(xml.folder_path,spkFiles(iShank).name),'r'); - wav = fread(fid,[1 inf],'int16=>int16'); - try %bug in some spk files... wrong number of samples? - wav = reshape(wav,chansPerSpikeGrp,nSamples,[]); - catch - if strcmp(getWaveforms,'force') - wav = nan(chansPerSpikeGrp,nSamples,length(clu)); - display([spkFiles(iShank).name,' error.']) - else - error(['something is wrong with ',spkFiles(iShank).name,... - ' Use ''getWaveforms'', false to skip waveforms or ',... - '''getWaveforms'', ''force'' to write nans on bad shanks.']) - end - end - wav = permute(wav,[3 1 2]); - - - - - %% Get the DynamicTableRegion field for each shank - - electrodes_field = types.core.DynamicTableRegion('table',types.untyped.ObjectView('/general/extracellular_ephys/electrodes'),'description',['shank' num2str(iShank) ' region'],'data',nwb.general_extracellular_ephys_electrodes.id.data(find(shank == iShank)')); - SpikeEventSeries = types.core.SpikeEventSeries('data', wav, 'electrodes', electrodes_field, 'timestamps', res./ spikes.samplingRate); - - %% This section assigns the spike-waveforms in the .NWB - ecephys.nwbdatainterface.set(['SpikeEventSeries' num2str(iShank)],SpikeEventSeries); - - end - - - cells = unique(clu); - % remove MUA and NOISE clusters... - cells(cells==0) = []; - cells(cells==1) = []; % consider adding MUA as another input argument...? - - - for c = 1:length(cells) - spikes.UID(count) = count; % this only works if all shanks are loaded... how do we optimize this? - ind = find(clu == cells(c)); - spikes.times{count} = res(ind) ./ spikes.samplingRate; - spikes.shankID(count) = shankID; - spikes.cluID(count) = cells(c); - - %Waveforms - if any(getWaveforms) - wvforms = squeeze(mean(wav(ind,:,:)))-mean(mean(mean(wav(ind,:,:)))); % mean subtract to account for slower (theta) trends - if prod(size(wvforms))==length(wvforms)%in single-channel groups wvforms will squeeze too much and will have amplitude on D1 rather than D2 - wvforms = wvforms';%fix here - end - for t = 1:size(wvforms,1) - [a(t) b(t)] = max(abs(wvforms(t,:))); - end - [aa bb] = max(a,[],2); - spikes.rawWaveform{count} = wvforms(bb,:); - spikes.maxWaveformCh(count) = spkGrpChans(bb); % Use this in Brainstorm - % %Regions (needs waveform peak) - % if isfield(xml,'region') %if there is regions field in your metadata - % spikes.region{count} = 'unknown'; - % elseif isfield(xml,'units') %if no regions, but unit region from xml via Loadparamteres - % %Find the xml Unit that matches group/cluster - % unitnum = cellfun(@(X,Y) X==spikes.shankID(count) && Y==spikes.cluID(count),... - % {sessionInfo.Units(:).spikegroup},{sessionInfo.Units(:).cluster}); - % if sum(unitnum) == 0 - % display(['xml Missing Unit - spikegroup: ',... - % num2str(spikes.shankID(count)),' cluster: ',... - % num2str(spikes.cluID(count))]) - % spikes.region{count} = 'missingxml'; - % else %possible future bug: two xml units with same group/clu... - % spikes.region{count} = sessionInfo.Units(unitnum).structure; - % end - % end - clear a aa b bb - end - - count = count + 1; - - end - - ecephys.description = 'intermediate data from extracellular electrophysiology recordings, e.g., LFP'; - nwb.processing.set('ecephys', ecephys); - end - - - % Serialize spiketimes and cluIDs - spike_times = []; - spike_times_index = []; - - current_index = 0; - for iNeuron = 1:length(spikes.UID) - spike_times = [spike_times ; spikes.times{iNeuron}]; - spike_times_index = [spike_times_index; int64(length(spikes.times{iNeuron})+current_index)]; - current_index = spike_times_index(end); - end - - - % electrode_group - Assigns the group_object_view that was defined above at - % the electrodes, to specific neurons - I need to find how each neuron is - % assigned to a shank - electrode_group = []; - shank_that_neurons_belongs_to = zeros(length(spikes.UID),1); - for iNeuron = 1:length(spikes.UID) - shank_that_neurons_belongs_to(iNeuron) = str2double(xml.units.unit{iNeuron}.group.Text); - first_electrode_in_shank = find(shank == shank_that_neurons_belongs_to(iNeuron)); - first_electrode_in_shank = first_electrode_in_shank(1); - electrode_group = [electrode_group; group_object_view(first_electrode_in_shank)]; - end - - electrode_group = types.core.VectorData('data', electrode_group, 'description','the electrode group that each spike unit came from'); - - % Initialize the fields needed - spike_times = types.core.VectorData ('data', spike_times, 'description', 'the spike times for each unit'); - spike_times_index = types.core.VectorIndex ('data', spike_times_index, 'target', types.untyped.ObjectView('/units/spike_times')); % The ObjectView links the indices to the spike times - id = types.core.ElementIdentifiers('data', [0:length(xml.units.unit)-1]'); - - - % FOR THE VECTORDATA I NEED FILE: DG_all_6__UnitFeatureSummary_add (PROBABLY - ACCORDING TO THE CONVERTER) - - - % First, instantiate the table, listing all of the columns that will be - % added and us the |'id'| argument to indicate the number of rows. Ifa - % value is indexed, only the column name is included, not the index. For - % instance, |'spike_times_index'| is not added to the array. - nwb.units = types.core.Units( ... - 'electrode_group', electrode_group, 'electrodes', [], 'electrodes_index', [], 'obs_intervals', [], 'obs_intervals_index', [], ... - 'spike_times', spike_times, 'spike_times_index', spike_times_index, 'waveform_mean', [], 'waveform_sd', [], ... - 'colnames', {'shank_id'; 'spike_times'; 'electrode_group'; 'cell_type'; 'global_id'; 'max_electrode'}, ... - 'description', 'Generated from Neuroscope2NWB', 'id', id, 'vectordata', [], 'vectorindex', []); - - end - - - - function nwb = getEvents(xml,nwb) - %% Add events: nwb2.stimulus_presentation - - eventFiles = dir([xml.folder_path filesep '*.evt']); - - for iFile = 1:length(eventFiles) - - if ~strcmp(eventFiles(iFile).name,'YutaMouse41-150903.DS1.ch0.evt') && ~strcmp(eventFiles(iFile).name,'YutaMouse41-150903.DS2.ch0.evt') % Ignore those two files - - events = LoadEvents(eventFiles(iFile).name); % Load Events is a function from the Buzcode - THIS DOESN'T LOAD ANYTHING FOR 'PulseStim_0V_10021ms_LD0' - maybe because it has a single entrty??? - - if ~isempty(events.time) && ~isempty(events.description) - AnnotationSeries = types.core.AnnotationSeries('data',events.description,'timestamps',events.time); - nwb.stimulus_presentation.set(events.description{1}, AnnotationSeries); - end - end - - end - end - - - - function nwb = getBehavior(xml,nwb) - - %% Add behavioral data: nwb2.processing.get('behavior').nwbdatainterface - - behavior = types.core.ProcessingModule; - - behavioralFiles = dir([xml.folder_path filesep '*.position']); - - time_epochs = repmat(struct('label','','start_times',0,'stop_times',0),length(behavioralFiles),1); - - for iFile = 1:length(behavioralFiles) - - % The label of the behavior - behavioral_Label = strsplit(behavioralFiles(iFile).name,'__'); - behavioral_Label = erase(behavioral_Label{2},'.mat'); - - position_signals = load(behavioralFiles(iFile).name); - - % Some behavioral signals might have more than one signal in them - field_names = fieldnames(position_signals); - - the_position_field_NWB = types.core.Position; - - for iField = 1:length(field_names) - - behavioral_timestamps = position_signals.(field_names{iField})(:,1); - position_coordinates = position_signals.(field_names{iField})(:,2:end); - spatial_series = types.core.SpatialSeries('data', position_coordinates, 'timestamps', behavioral_timestamps, 'reference_frame', 'unknown', 'data_conversion', 1); - - the_position_field_NWB.spatialseries.set(field_names{iField}, spatial_series); - end - - behavior.nwbdatainterface.set(behavioral_Label,the_position_field_NWB); - - time_epochs(iFile).start_times = behavioral_timestamps(1,1); - time_epochs(iFile).stop_times = behavioral_timestamps(end,1); - time_epochs(iFile).label = behavioral_Label; - - end - - behavior.description = 'Behavioral signals'; - nwb.processing.set('behavior', behavior); - - end - - - - function nwb = getElectrophysiology(xml,nwb) - - lfpFile = dir([xml.folder_path filesep '*.eeg']); - - if length(lfpFile)>1 - error('More than one .eeg files are present here. Weird') - end - - - % Get the samples number, based on the size of the file - % Check for the precision that samples are saved first - - hdr.nBits = str2double(xml.acquisitionSystem.nBits.Text); - hdr.nChannels = str2double(xml.acquisitionSystem.nChannels.Text); - hdr.sRateOrig = str2double(xml.acquisitionSystem.samplingRate.Text); - hdr.Gain = str2double(xml.acquisitionSystem.amplification.Text); - hdr.sRateLfp = str2double(xml.fieldPotentials.lfpSamplingRate.Text); - - - % Get data type - switch lower(hdr.nBits) - case 16; - hdr.byteSize = 2; - hdr.byteFormat = 'int16'; - case 32; - hdr.byteSize = 4; - hdr.byteFormat = 'int32'; - end - % Guess the number of time points based on the file size - dirInfo = dir(lfpFile.name); - hdr.nSamples = floor(dirInfo.bytes ./ (hdr.nChannels * hdr.byteSize)); - - - - - - lfp_data = bz_LoadBinary(lfpFile.name, 'duration',Inf, 'frequency',hdr.sRateLfp,'nchannels',hdr.nChannels, 'channels', [1:hdr.nChannels]); % nSamples x 64 - - - % If the electrode Information has not already bwwn filled, - % do it now - if isempty(nwb.general_extracellular_ephys_electrodes) - nwb = Neuroscope2NWB.getElectrodeInfo(xml,nwb); - end - - - electrodes_field = types.core.DynamicTableRegion('table',types.untyped.ObjectView('/general/extracellular_ephys/electrodes'),'description','electrode table reference','data',nwb.general_extracellular_ephys_electrodes.id.data); - - lfp = types.core.ElectricalSeries('data', lfp_data', 'electrodes',electrodes_field, 'description', 'lfp signal for all shank electrodes', 'starting_time', 0, 'starting_time_rate', hdr.sRateLfp); - % I TRANSPOSED THE MATRIX HERE TO BE COMPATIBLE WITH THE FUNCTION BZ_GET_LFP - - LFP = types.core.LFP; - LFP.electricalseries.set('lfp',lfp); - - % Check if the ecephys field is already created - if ~ismember(keys(nwb.processing),'ecephys') - ecephys = types.core.ProcessingModule('description', ''); - ecephys.description = 'intermediate data from extracellular electrophysiology recordings, e.g., LFP'; - nwb.processing.set('ecephys', ecephys); - end - - nwb.processing.get('ecephys').nwbdatainterface.set('LFP', LFP); - - end - - - function nwb = getEpochs(xml,nwb) - - - %% Add Epochs - % The epochs are based on the separate behavioral files - % Each separate behavioral file = 1 epoch - - behavioralFiles = dir([xml.folder_path filesep '*.position']); - - time_epochs = repmat(struct('label','','start_times',0,'stop_times',0),length(behavioralFiles),1); - - for iFile = 1:length(behavioralFiles) - - for iField = 1:length(field_names) - behavioral_timestamps = position_signals.(field_names{iField})(:,1); - time_epochs(iFile).start_times = behavioral_timestamps(1,1); - time_epochs(iFile).stop_times = behavioral_timestamps(end,1); - time_epochs(iFile).label = behavioral_Label; - end - - end - - id_epochs = types.core.ElementIdentifiers('data',int64(0:length(time_epochs)-1)'); - - start_time_epochs = types.core.VectorData('description','Starting timepoint of Each Epoch','data',[time_epochs.start_times]'); - stop_time_epochs = types.core.VectorData('description','Ending timepoint of Each Epoch','data',[time_epochs.stop_times]'); - - - vectordata_epochs_label = types.core.VectorData('description','Label of Epoch','data',{time_epochs.label}'); - vectordata_epochs = types.untyped.Set('label', vectordata_epochs_label); - - - %%%%%% THE VECTORDATA IS IGNORED IN HERE - THE LABELS ARE NOT SAVED - - %%%%%% CHECK AGAIN - %%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - intervals_epochs = types.core.TimeIntervals('start_time',start_time_epochs,'stop_time',stop_time_epochs,... - 'colnames',{'start_time';'stop_time';'label'},... - 'description','experimental epochs','id',id_epochs, 'vectordata',vectordata_epochs); - - nwb.intervals_epochs = intervals_epochs; - - - end - - - - - function nwb = getTrials(xml,nwb) - - - %% Add Trials - - % This file holds a matrix with the trial info - trialsFile = dir([xml.folder_path filesep '*Run.mat']); - trialsInfo = load(trialsFile.name); - the_field = fieldnames(trialsInfo); - trialsInfo = trialsInfo.(the_field{1}); - - % This file holds a cell array with the labels for the matrix above (...) - runFile = dir([xml.folder_path filesep '*RunInfo.mat']); - runInfo = load(runFile.name); - the_field = fieldnames(runInfo); - runInfo = runInfo.(the_field{1}); - - - - start_time_trials = types.core.VectorData('description','Starting timepoint of Each Trial','data',trialsInfo(:,1)); - stop_time_trials = types.core.VectorData('description','Ending timepoint of Each Trial','data',trialsInfo(:,2)); - - id_trials = types.core.ElementIdentifiers('data',int64(0:size(trialsInfo,1)-1)'); - - conditions_trials = cell(size(trialsInfo,1),1); - for iTrial = 1:size(trialsInfo,1) - conditions_trials{iTrial} = runInfo{find(trialsInfo(iTrial,3:4))+2}; - end - - vectordata_trials = types.untyped.Set('both_visit', trialsInfo(:,7), 'condition', conditions_trials, 'error_run', trialsInfo(:,5), 'stim_run',trialsInfo(:,6)); - - intervals_trials = types.core.TimeIntervals('start_time',start_time_trials,'stop_time',stop_time_trials,... - 'colnames',{'start_time';'stop_time';'error_run';'stim_run';'both_visit';'condition'},... - 'description','experimental trials','id',id_trials, 'vectordata',vectordata_trials); - - nwb.intervals_trials = intervals_trials; - - end - - - - function nwb = special_YutaMouse_recordings(xml,nwb) - - %% Add raw recordings - - % value taken from Yuta's spreadsheet - - % HOW ABOUT POSITION0 - POSITION1 CHANNELS??? - - hdr.nChannels = str2double(xml.acquisitionSystem.nChannels.Text); - hdr.sRateLfp = str2double(xml.fieldPotentials.lfpSamplingRate.Text); - - lfpFile = dir([xml.folder_path filesep '*.eeg']); - - special_electrode_labels = {'ch_wait','ch_arm','ch_solL','ch_solR','ch_dig1','ch_dig2','ch_entL','ch_entR','ch_SsolL','ch_SsolR'}; - special_electrode_indices = [79,78,76,77,65,68,72,71,73,70]; - - acquisition = types.untyped.Set; - for iSpecialElectrode = 1:length(special_electrode_labels) - - special_Electrode_data = bz_LoadBinary(lfpFile.name, 'duration',Inf, 'frequency',hdr.sRateLfp,'nchannels',hdr.nChannels, 'channels', special_electrode_indices(iSpecialElectrode)); - single_Electrode = types.core.TimeSeries('description','environmental electrode recorded inline with neural data','data',special_Electrode_data,'starting_time', 0, 'starting_time_rate', hdr.sRateLfp, 'data_unit','V'); - acquisition.set(special_electrode_labels{iSpecialElectrode}, single_Electrode); - end - - nwb.acquisition = acquisition; - - end - - - end - -end From 2cd4cb7d3f661fc2d640b3fe682bd1fbe3a774de Mon Sep 17 00:00:00 2001 From: Konstantinos Date: Wed, 22 May 2019 00:52:46 -0400 Subject: [PATCH 5/5] Merged Buzcode and NWB loading functions - Separated add2NWB functions according to standard and non-standard files --- io/NWB/GeneralInfo.m | 58 ++ io/NWB/GetXMLInfo.m | 27 + io/NWB/Neuroscope/addUnitsInfo_Neuroscope.m | 249 ++++++++ io/NWB/Yuta/YutaMouse41_toNWB.m | 45 ++ io/NWB/Yuta/addBehavior_Yuta.m | 53 ++ io/NWB/Yuta/addEpochs_Yuta.m | 60 ++ io/NWB/Yuta/addEvents_Yuta.m | 13 + io/NWB/Yuta/addSpecial_YutaMouse_recordings.m | 22 + io/NWB/Yuta/addTrials_Yuta.m | 53 ++ io/NWB/YutaMouse41_toNWB.m | 41 -- io/NWB/addBehavior.m | 77 +++ io/NWB/addElectrodeInfo.m | 88 +++ io/NWB/addElectrophysiology.m | 84 +++ io/NWB/addEvents.m | 84 +++ io/NWB/addTrials.m | 152 +++++ io/NWB/addUnitsInfo.m | 71 +++ io/NWB/bz_GetLFP_NWB.m | 201 ------ io/NWB/bz_GetSpikes_NWB.m | 245 -------- io/NWB/bz_LoadBehavior_NWB.m | 436 ------------- io/NWB/bz_LoadEvents_NWB.m | 121 ---- ...script_to_add_fields_from_Buzcode_format.m | 42 ++ io/NWB/get_SessionInfoFromNWB.m | 63 +- io/bz_GetLFP.m | 90 ++- io/bz_GetSpikes.m | 583 +++++++++++------- io/bz_LoadBehavior.m | 227 ++++++- io/bz_LoadCellinfo.m | 106 +++- io/bz_LoadEvents.m | 126 +++- 27 files changed, 2095 insertions(+), 1322 deletions(-) create mode 100644 io/NWB/GeneralInfo.m create mode 100644 io/NWB/GetXMLInfo.m create mode 100644 io/NWB/Neuroscope/addUnitsInfo_Neuroscope.m create mode 100644 io/NWB/Yuta/YutaMouse41_toNWB.m create mode 100644 io/NWB/Yuta/addBehavior_Yuta.m create mode 100644 io/NWB/Yuta/addEpochs_Yuta.m create mode 100644 io/NWB/Yuta/addEvents_Yuta.m create mode 100644 io/NWB/Yuta/addSpecial_YutaMouse_recordings.m create mode 100644 io/NWB/Yuta/addTrials_Yuta.m delete mode 100644 io/NWB/YutaMouse41_toNWB.m create mode 100644 io/NWB/addBehavior.m create mode 100644 io/NWB/addElectrodeInfo.m create mode 100644 io/NWB/addElectrophysiology.m create mode 100644 io/NWB/addEvents.m create mode 100644 io/NWB/addTrials.m create mode 100644 io/NWB/addUnitsInfo.m delete mode 100644 io/NWB/bz_GetLFP_NWB.m delete mode 100644 io/NWB/bz_GetSpikes_NWB.m delete mode 100644 io/NWB/bz_LoadBehavior_NWB.m delete mode 100644 io/NWB/bz_LoadEvents_NWB.m create mode 100644 io/NWB/example_script_to_add_fields_from_Buzcode_format.m diff --git a/io/NWB/GeneralInfo.m b/io/NWB/GeneralInfo.m new file mode 100644 index 00000000..7041511d --- /dev/null +++ b/io/NWB/GeneralInfo.m @@ -0,0 +1,58 @@ +function nwb = GeneralInfo(xml) + % Adds info in: nwb.general_subject + % Konstantinos Nasiotis 2019 + + %% General Info + nwb_version = '2.0b'; + + session_start_time = datetime(xml.generalInfo.date.Text, ... + 'Format', 'yyyy-MM-dd''T''HH:mm:ssZZ', ... + 'TimeZone', 'local'); + timestamps_reference_time = datetime(xml.generalInfo.date.Text, ... + 'Format', 'yyyy-MM-dd''T''HH:mm:ssZZ', ... + 'TimeZone', 'local'); + + + %% Check for a .lfp or a .eeg file. Use the creation date of that file + % to store in the NWB + [ff, basename] = fileparts(xml.folder_path); + lfpFile = dir([xml.folder_path filesep basename '*.lfp']); + if length(lfpFile)>1 + disp('More than one .eeg files are present here. No Electrophysiology signals were added') + return + elseif length(lfpFile)==0 + lfpFile = dir([xml.folder_path filesep basename '*.eeg']); + if length(lfpFile)>1 + disp('More than one .lfp files are present here. No Electrophysiology signals were added') + return + elseif length(lfpFile)==0 + disp('No .eeg or .lfp files are present in the selected directory. No Electrophysiology signals were added') + return + end + end + + + file_create_date = datetime(lfpFile.date, ... + 'Format', 'yyyy-MM-dd''T''HH:mm:ssZZ', ... + 'TimeZone', 'local'); + + %% + nwb = nwbfile( ... + 'session_description' , 'Mouse in open exploration and theta maze', ... + 'identifier' , xml.name, ... + 'session_start_time' , session_start_time,... + 'file_create_date' , file_create_date,... + 'general_experimenter' , xml.generalInfo.experimenters.Text,... + 'general_session_id' , xml.name,... + 'general_institution' , 'NYU' ,... + 'general_lab' , 'Buzsaki',... + 'subject' , 'YutaMouse',... + 'general_related_publications' , 'DOI:10.1016/j.neuron.2016.12.011',... + 'timestamps_reference_time' , session_start_time); + + nwb.general_subject = types.core.Subject( ... + 'description', 'mouse 5', 'genotype', 'POMC-Cre::Arch', 'age', '9 months', ... + 'sex', 'M', 'subject_id', xml.name, 'species', 'Mus musculus'); + + disp('General info added..') +end \ No newline at end of file diff --git a/io/NWB/GetXMLInfo.m b/io/NWB/GetXMLInfo.m new file mode 100644 index 00000000..1c8825ab --- /dev/null +++ b/io/NWB/GetXMLInfo.m @@ -0,0 +1,27 @@ +function xml = GetXMLInfo(folder_path) + + %% This uses an .xml importer downloaded from MathWorks - File Exchange + % https://www.mathworks.com/matlabcentral/fileexchange/28518-xml2struct + % The loadxml from the Buzcode repo gave errors + + [previous_path, name] = fileparts(folder_path); + + all_files_in_folder = dir(folder_path); + + iXML = []; + for iFile = 1:length(all_files_in_folder) + if strfind(all_files_in_folder(iFile).name,'.xml') + iXML = [iXML iFile]; + end + end + if isempty(iXML) + error 'There are no .xml files in this folder' + elseif length(iXML)>1 + error 'There is more than one .xml in this folder' + end + + xml = xml2struct([folder_path filesep all_files_in_folder(iXML).name]); + xml = xml.parameters; + xml.folder_path = folder_path; + xml.name = name; +end \ No newline at end of file diff --git a/io/NWB/Neuroscope/addUnitsInfo_Neuroscope.m b/io/NWB/Neuroscope/addUnitsInfo_Neuroscope.m new file mode 100644 index 00000000..ce71df61 --- /dev/null +++ b/io/NWB/Neuroscope/addUnitsInfo_Neuroscope.m @@ -0,0 +1,249 @@ + function nwb = addUnitsInfo_Neuroscope(xml, nwb) + %% Add the units info (copied from bz_GetSpikes) + % Adds unit info in: nwb.units + + % This code takes the unit information from the .clu, .res, .spk files + + + %% Get unit information + getWaveforms = 1; % Set this to true if you want to add waveforms on the NWB file + + + spikes.samplingRate = str2double(xml.acquisitionSystem.samplingRate.Text); + + + disp('loading spikes from clu/res/spk files..') + % find res/clu/fet/spk files here + cluFiles = dir([xml.folder_path filesep '*.clu*']); + resFiles = dir([xml.folder_path filesep '*.res*']); + if any(getWaveforms) + spkFiles = dir([xml.folder_path filesep '*.spk*']); + end + + % remove *temp*, *autosave*, and *.clu.str files/directories + tempFiles = zeros(length(cluFiles),1); + for i = 1:length(cluFiles) + dummy = strsplit(cluFiles(i).name, '.'); % Check whether the component after the last dot is a number or not. If not, exclude the file/dir. + if ~isempty(findstr('temp',cluFiles(i).name)) | ~isempty(findstr('autosave',cluFiles(i).name)) | isempty(str2num(dummy{length(dummy)})) | find(contains(dummy, 'clu')) ~= length(dummy)-1 + tempFiles(i) = 1; + end + end + cluFiles(tempFiles==1)=[]; + tempFiles = zeros(length(resFiles),1); + for i = 1:length(resFiles) + if ~isempty(findstr('temp',resFiles(i).name)) | ~isempty(findstr('autosave',resFiles(i).name)) + tempFiles(i) = 1; + end + end + if any(getWaveforms) + resFiles(tempFiles==1)=[]; + tempFiles = zeros(length(spkFiles),1); + for i = 1:length(spkFiles) + if ~isempty(findstr('temp',spkFiles(i).name)) | ~isempty(findstr('autosave',spkFiles(i).name)) + tempFiles(i) = 1; + end + end + spkFiles(tempFiles==1)=[]; + end + + if isempty(cluFiles) + disp('no clu files found...') + spikes = []; + return + end + + + % ensures we load in sequential order (forces compatibility with FMAT + % ordering) + for i = 1:length(cluFiles) + temp = strsplit(cluFiles(i).name,'.'); + shanks(i) = str2num(temp{length(temp)}); + end + [shanks ind] = sort(shanks); + cluFiles = cluFiles(ind); %Bug here if there are any files x.clu.x that are not your desired clus + resFiles = resFiles(ind); + if any(getWaveforms) + spkFiles = spkFiles(ind); + end + + % check if there are matching #'s of files + if length(cluFiles) ~= length(resFiles) && length(cluFiles) ~= length(spkFiles) + error('found an incorrect number of res/clu/spk files...') + end + + % use the .clu files to get spike ID's and generate UID and spikeGroup + % use the .res files to get spike times + count = 1; + + ecephys = types.core.ProcessingModule; + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % This section is copied from the ElectrodesInfo + nShanks = length(xml.spikeDetection.channelGroups.group); + groups = xml.spikeDetection.channelGroups.group; % Use this for simplicity + all_shank_channels = cell(nShanks,1); % This will hold the channel numbers that belong in each shank + shank = []; + group_object_view = []; + + for iGroup = 1:nShanks + % Get all_shank_channls again for iGroup = 1:nShanks + for iChannel = 1:length(groups{iGroup}.channels.channel) + all_shank_channels{iGroup} = [all_shank_channels{iGroup} str2double(groups{iGroup}.channels.channel{iChannel}.Text)]; + shank = [shank iGroup]; + group_object_view = [group_object_view; types.untyped.ObjectView(['/general/extracellular_ephys/' ['shank' num2str(iGroup)]])]; + end + end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + for iShank=1:length(cluFiles) + disp(['working on ' cluFiles(iShank).name]) + + temp = strsplit(cluFiles(iShank).name,'.'); + shankID = str2num(temp{length(temp)}); %shankID is the spikegroup number + clu = load(fullfile(xml.folder_path,cluFiles(iShank).name)); + clu = clu(2:end); % toss the first sample to match res/spk files + res = load(fullfile(xml.folder_path,resFiles(iShank).name)); + spkGrpChans = all_shank_channels{iShank}; + + if any(getWaveforms) && sum(clu)>0 %bug fix if no clusters + nSamples = str2double(xml.spikeDetection.channelGroups.group{iShank}.nSamples.Text); + % load waveforms + chansPerSpikeGrp = length(all_shank_channels{iShank}); + fid = fopen(fullfile(xml.folder_path,spkFiles(iShank).name),'r'); + wav = fread(fid,[1 inf],'int16=>int16'); + try %bug in some spk files... wrong number of samples? + wav = reshape(wav,chansPerSpikeGrp,nSamples,[]); + catch + if strcmp(getWaveforms,'force') + wav = nan(chansPerSpikeGrp,nSamples,length(clu)); + display([spkFiles(iShank).name,' error.']) + else + error(['something is wrong with ',spkFiles(iShank).name,... + ' Use ''getWaveforms'', false to skip waveforms or ',... + '''getWaveforms'', ''force'' to write nans on bad shanks.']) + end + end + wav = permute(wav,[3 1 2]); + + %% Get the DynamicTableRegion field for each shank + + % First check if the electrodes field has been filled + if isempty(nwb.general_extracellular_ephys_electrodes) + nwb = Neuroscope2NWB.getElectrodeInfo(xml, nwb); + end + + electrodes_field = types.core.DynamicTableRegion('table',types.untyped.ObjectView('/general/extracellular_ephys/electrodes'),'description',['shank' num2str(iShank) ' region'],'data',nwb.general_extracellular_ephys_electrodes.id.data(find(shank == iShank)')); + SpikeEventSeries = types.core.SpikeEventSeries('data', wav, 'electrodes', electrodes_field, 'timestamps', res./ spikes.samplingRate); + + %% This section assigns the spike-waveforms in the .NWB + ecephys.nwbdatainterface.set(['SpikeEventSeries' num2str(iShank)],SpikeEventSeries); + + end + + + cells = unique(clu); + % remove MUA and NOISE clusters... + cells(cells==0) = []; + cells(cells==1) = []; % consider adding MUA as another input argument...? + + + for c = 1:length(cells) + spikes.UID(count) = count; % this only works if all shanks are loaded... how do we optimize this? + ind = find(clu == cells(c)); + spikes.times{count} = res(ind) ./ spikes.samplingRate; + spikes.shankID(count) = shankID; + spikes.cluID(count) = cells(c); + + %Waveforms + if any(getWaveforms) + wvforms = squeeze(mean(wav(ind,:,:)))-mean(mean(mean(wav(ind,:,:)))); % mean subtract to account for slower (theta) trends + if prod(size(wvforms))==length(wvforms)%in single-channel groups wvforms will squeeze too much and will have amplitude on D1 rather than D2 + wvforms = wvforms';%fix here + end + for t = 1:size(wvforms,1) + [a(t) b(t)] = max(abs(wvforms(t,:))); + end + [aa bb] = max(a,[],2); + spikes.rawWaveform{count} = wvforms(bb,:); + spikes.maxWaveformCh(count) = spkGrpChans(bb); % Use this in Brainstorm + % %Regions (needs waveform peak) + % if isfield(xml,'region') %if there is regions field in your metadata + % spikes.region{count} = 'unknown'; + % elseif isfield(xml,'units') %if no regions, but unit region from xml via Loadparamteres + % %Find the xml Unit that matches group/cluster + % unitnum = cellfun(@(X,Y) X==spikes.shankID(count) && Y==spikes.cluID(count),... + % {sessionInfo.Units(:).spikegroup},{sessionInfo.Units(:).cluster}); + % if sum(unitnum) == 0 + % display(['xml Missing Unit - spikegroup: ',... + % num2str(spikes.shankID(count)),' cluster: ',... + % num2str(spikes.cluID(count))]) + % spikes.region{count} = 'missingxml'; + % else %possible future bug: two xml units with same group/clu... + % spikes.region{count} = sessionInfo.Units(unitnum).structure; + % end + % end + clear a aa b bb + end + + count = count + 1; + + end + + ecephys.description = 'intermediate data from extracellular electrophysiology recordings, e.g., LFP'; + nwb.processing.set('ecephys', ecephys); + end + + + % Serialize spiketimes and cluIDs + spike_times = []; + spike_times_index = []; + + current_index = 0; + for iNeuron = 1:length(spikes.UID) + spike_times = [spike_times ; spikes.times{iNeuron}]; + spike_times_index = [spike_times_index; int64(length(spikes.times{iNeuron})+current_index)]; + current_index = spike_times_index(end); + end + + + % electrode_group - Assigns the group_object_view that was defined above at + % the electrodes, to specific neurons - I need to find how each neuron is + % assigned to a shank + electrode_group = []; + shank_that_neurons_belongs_to = zeros(length(spikes.UID),1); + for iNeuron = 1:length(spikes.UID) + shank_that_neurons_belongs_to(iNeuron) = str2double(xml.units.unit{iNeuron}.group.Text); + first_electrode_in_shank = find(shank == shank_that_neurons_belongs_to(iNeuron)); + first_electrode_in_shank = first_electrode_in_shank(1); + electrode_group = [electrode_group; group_object_view(first_electrode_in_shank)]; + end + + electrode_group = types.core.VectorData('data', electrode_group, 'description','the electrode group that each spike unit came from'); + + % Initialize the fields needed + spike_times = types.core.VectorData ('data', spike_times, 'description', 'the spike times for each unit'); + spike_times_index = types.core.VectorIndex ('data', spike_times_index, 'target', types.untyped.ObjectView('/units/spike_times')); % The ObjectView links the indices to the spike times + id = types.core.ElementIdentifiers('data', [0:length(xml.units.unit)-1]'); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% THIS GAVE AN ERROR WHEN ASSIGNING CELL ARRAY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + waveform_mean = types.core.VectorData('data', spikes.rawWaveform, 'description', 'The mean Waveform for each unit'); + waveform_mean = []; + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + %% Fill the units fields + nwb.units = types.core.Units( ... + 'electrode_group', electrode_group, 'electrodes', [], 'electrodes_index', [], 'obs_intervals', [], 'obs_intervals_index', [], ... + 'spike_times', spike_times, 'spike_times_index', spike_times_index, 'waveform_mean', waveform_mean, 'waveform_sd', [], ... + 'colnames', {'shank_id'; 'spike_times'; 'electrode_group'; 'cell_type'; 'global_id'; 'max_electrode'}, ... + 'description', 'Generated from Neuroscope2NWB', 'id', id, 'vectorindex', []); + + %% Extra Unit Info + % FOR THE VECTORDATA, IDEALLY I NEED FILE: DG_all_6__UnitFeatureSummary_add (ACCORDING TO BEN'S CONVERTER - THIS HOLDS INFO ABOUT THE CELL_TYPE, GLOBAL_ID) + nwb.units.vectordata.set('cluID', types.core.VectorData('description', 'cluster ID', 'data', spikes.cluID)); + nwb.units.vectordata.set('maxWaveformCh', types.core.VectorData('description', 'The electrode where each unit showed maximum Waveform', 'data', spikes.maxWaveformCh)); + + disp('Spikes info added..') +end \ No newline at end of file diff --git a/io/NWB/Yuta/YutaMouse41_toNWB.m b/io/NWB/Yuta/YutaMouse41_toNWB.m new file mode 100644 index 00000000..108af3bd --- /dev/null +++ b/io/NWB/Yuta/YutaMouse41_toNWB.m @@ -0,0 +1,45 @@ +%% Example conversion to NWB of the YutaMouse41-150903 dataset + +% Folder that contains all the files - This is the only input needed +folder_path = 'F:\NWBtoBuzcode\YutaMouse41-150903'; +% folder_path = 'C:\Users\McGill\Documents\GitHub\buzcode\tutorials\exampleDataStructs\20170505_396um_0um_merge'; + + + +%% Start Adding fields to NWB + +% Get info from the xml file +xml = GetXMLInfo(folder_path); + +% Add general info to the NWB file +nwb = GeneralInfo(xml); + +% Add electrode info +nwb = addElectrodeInfo(xml, nwb); + +% Add units info - By default, the spike waveforms are added to the file +nwb = addUnitsInfo_Neuroscope(xml,nwb); + +% Add stimulation events +nwb = addEvents_Yuta(xml,nwb); + +% Add behavioral info/channels +nwb = addBehavior_Yuta(xml,nwb); + +% Add electrophysiological channels +nwb = addElectrophysiology(xml, nwb); + +% Add epochs +nwb = addEpochs_Yuta(xml,nwb); + +% Add trials +nwb = addTrials_Yuta(xml,nwb); + +% Add channels based on the Yuta spreadsheet +nwb = addSpecial_YutaMouse_recordings(xml,nwb); + + +%% Export to nwb +% nwbExport(nwb, 'YutaMouse41_converted.nwb') +nwbExport(nwb, 'F:\NWBtoBuzcode\YutaMouse41-150903\YutaMouse41.nwb') + diff --git a/io/NWB/Yuta/addBehavior_Yuta.m b/io/NWB/Yuta/addBehavior_Yuta.m new file mode 100644 index 00000000..a6a7ca7f --- /dev/null +++ b/io/NWB/Yuta/addBehavior_Yuta.m @@ -0,0 +1,53 @@ +function nwb = addBehavior_Yuta(xml,nwb) + + %% Add behavioral data: nwb2.processing.get('behavior').nwbdatainterface + % Check for YutaMouse behavioral files (*position*) + + behavioralFiles = dir(fullfile(xml.folder_path, '*position*')); + + if ~isempty(behavioralFiles) + + behavior_NWB = types.core.ProcessingModule; + + + time_epochs = repmat(struct('label','','start_times',0,'stop_times',0),length(behavioralFiles),1); + + for iFile = 1:length(behavioralFiles) + + % The label of the behavior + behavioral_Label = strsplit(behavioralFiles(iFile).name,'__'); + behavioral_Label = erase(behavioral_Label{2},'.mat'); + + position_signals = load(fullfile(behavioralFiles(iFile).folder, behavioralFiles(iFile).name)); + + % Some behavioral signals might have more than one signal in them + field_names = fieldnames(position_signals); + + the_position_field_NWB = types.core.Position; + + for iField = 1:length(field_names) + + behavioral_timestamps = position_signals.(field_names{iField})(:,1); + position_coordinates = position_signals.(field_names{iField})(:,2:end); + spatial_series = types.core.SpatialSeries('data', position_coordinates, 'timestamps', behavioral_timestamps, 'reference_frame', 'unknown', 'data_conversion', 1); + + the_position_field_NWB.spatialseries.set(field_names{iField}, spatial_series); + end + + behavior_NWB.nwbdatainterface.set(behavioral_Label,the_position_field_NWB); + + time_epochs(iFile).start_times = behavioral_timestamps(1,1); + time_epochs(iFile).stop_times = behavioral_timestamps(end,1); + time_epochs(iFile).label = behavioral_Label; + end + + behavior_NWB.description = 'Behavioral signals'; + nwb.processing.set('behavior', behavior_NWB); + + disp('Behavioral info added..') + + else + disp('No Yuta behavioral signals present in this directory (*position*)') + return + end +end \ No newline at end of file diff --git a/io/NWB/Yuta/addEpochs_Yuta.m b/io/NWB/Yuta/addEpochs_Yuta.m new file mode 100644 index 00000000..5ab71929 --- /dev/null +++ b/io/NWB/Yuta/addEpochs_Yuta.m @@ -0,0 +1,60 @@ +function nwb = addEpochs_Yuta(xml,nwb) + + %% Add Epochs + % The epochs are based on separate behavioral files + % Each separate behavioral file = 1 epoch + + behavioralFiles = dir([xml.folder_path filesep '*Position*']); + + if length(behavioralFiles) == 0 + disp('There are no *Position* files in this folder. No Epochs will be added') + return + end + + + intervals_epochs = types.core.TimeIntervals(); + id_epochs = types.core.ElementIdentifiers('data',1:length(behavioralFiles)); + + start_time_epochs = zeros(length(behavioralFiles),1); % Start time + stop_time_epochs = zeros(length(behavioralFiles),1); % Stop time + colnames = cell(length(behavioralFiles),1); % Labels + + for iFile = 1:length(behavioralFiles) + + % The label of the behavior + behavioral_Label = strsplit(behavioralFiles(iFile).name,'__'); + behavioral_Label = erase(behavioral_Label{2},'.mat'); + + position_signals = load(fullfile(behavioralFiles(iFile).folder,behavioralFiles(iFile).name)); + + % Some behavioral signal files might have more than one + % type of signals in them (e.g. twhl_linearized, twhl_norm) + field_names = fieldnames(position_signals); + + % NO NEED TO KEEP THE DATA - COMMENTING OUT +% the_position_field_NWB = types.untyped.Set(); +% for iField = 1:length(field_names) +% the_position_field_NWB.set(field_names{iField}, types.core.SpatialSeries('description', [field_names{iField} ' position signals from ' behavioral_Label ' epoch'], 'data', position_signals.(field_names{iField})(:,2:end))); +% end +% intervals_epochs.vectordata.set(behavioral_Label,types.core.VectorData('description', ['Position signals from ' behavioral_Label ' epoch'], 'data', the_position_field_NWB)); + + start_time_epochs(iFile) = position_signals.(field_names{1})(1,1); + stop_time_epochs(iFile) = position_signals.(field_names{1})(end,1); + colnames{iFile} = behavioral_Label; + + end + + % Transform doubles to VectorData + start_time_epochs = types.core.VectorData('description','Starting timestamp of epochs', 'data', start_time_epochs); + stop_time_epochs = types.core.VectorData('description','Stopping timestamp of epochs', 'data', stop_time_epochs); + + intervals_epochs.start_time = start_time_epochs; + intervals_epochs.stop_time = stop_time_epochs; + intervals_epochs.description = 'Experimental epochs'; + intervals_epochs.id = id_epochs; + intervals_epochs.colnames = colnames; + + nwb.intervals_epochs = intervals_epochs; + + disp('Epoch info added..') +end \ No newline at end of file diff --git a/io/NWB/Yuta/addEvents_Yuta.m b/io/NWB/Yuta/addEvents_Yuta.m new file mode 100644 index 00000000..e5ef2051 --- /dev/null +++ b/io/NWB/Yuta/addEvents_Yuta.m @@ -0,0 +1,13 @@ +function nwb = addEvents_Yuta(xml,nwb) + %% Add events: nwb2.stimulus_presentation + eventFiles = dir([xml.folder_path filesep '*.evt']); + + for iFile = 1:length(eventFiles) + events = LoadEvents(fullfile(eventFiles(iFile).folder,eventFiles(iFile).name)); % LoadEvents is a function from the Buzcode - THIS DOESN'T LOAD ANYTHING FOR 'PulseStim_0V_10021ms_LD0' - maybe because it has a single entrty??? + if ~isempty(events.time) && ~isempty(events.description) + AnnotationSeries = types.core.AnnotationSeries('data',events.description,'timestamps',events.time); + nwb.stimulus_presentation.set(events.description{1}, AnnotationSeries); + end + end + disp('Events added..') +end \ No newline at end of file diff --git a/io/NWB/Yuta/addSpecial_YutaMouse_recordings.m b/io/NWB/Yuta/addSpecial_YutaMouse_recordings.m new file mode 100644 index 00000000..21486788 --- /dev/null +++ b/io/NWB/Yuta/addSpecial_YutaMouse_recordings.m @@ -0,0 +1,22 @@ +function nwb = addSpecial_YutaMouse_recordings(xml,nwb) + %% Add raw recordings in: nwb.acquisition + + % values taken from Yuta's spreadsheet + % HOW ABOUT POSITION0 - POSITION1 CHANNELS??? + + hdr.nChannels = str2double(xml.acquisitionSystem.nChannels.Text); + hdr.sRateLfp = str2double(xml.fieldPotentials.lfpSamplingRate.Text); + + lfpFile = dir([xml.folder_path filesep '*.eeg']); + + special_electrode_labels = {'ch_wait','ch_arm','ch_solL','ch_solR','ch_dig1','ch_dig2','ch_entL','ch_entR','ch_SsolL','ch_SsolR'}; + special_electrode_indices = [79,78,76,77,65,68,72,71,73,70]; + + for iSpecialElectrode = 1:length(special_electrode_labels) + special_Electrode_data = bz_LoadBinary(fullfile(lfpFile.folder,lfpFile.name), 'duration',Inf, 'frequency',hdr.sRateLfp,'nchannels',hdr.nChannels, 'channels', special_electrode_indices(iSpecialElectrode)); + single_Electrode = types.core.TimeSeries('description','environmental electrode recorded inline with neural data','data',special_Electrode_data,'starting_time', 0, 'starting_time_rate', hdr.sRateLfp, 'data_unit','V'); + nwb.acquisition.set(special_electrode_labels{iSpecialElectrode}, single_Electrode); + end + + disp('Special YutaMouse channel info added..') +end \ No newline at end of file diff --git a/io/NWB/Yuta/addTrials_Yuta.m b/io/NWB/Yuta/addTrials_Yuta.m new file mode 100644 index 00000000..318786dd --- /dev/null +++ b/io/NWB/Yuta/addTrials_Yuta.m @@ -0,0 +1,53 @@ +function nwb = addTrials_Yuta(xml,nwb) + %% THIS SECTION LOADS THE INFO FROM *RUN.MAT - YUTA - NON STANDARD FORMAT + + % NOTE: Assumes the presence of a single *Run.mat and *RunInfo.mat file. If more are present, the code needs to change to store everything in: nwb.processing + + %% Add Trials + + % This file holds a matrix with the trial info + trialsFile = dir([xml.folder_path filesep '*Run.mat']); + + if ~isempty(trialsFile) + + trialsInfo = load(fullfile(trialsFile.folder, trialsFile.name)); + the_field = fieldnames(trialsInfo); + trialsInfo = trialsInfo.(the_field{1}); + + colname = the_field{1}; + + % This file holds a cell array with the labels for the matrix above (...) + runFile = dir([xml.folder_path filesep '*RunInfo.mat']); + runInfo = load(fullfile(runFile.folder,runFile.name)); + the_field = fieldnames(runInfo); + runInfo = runInfo.(the_field{1}); + + + start_time_trials = types.core.VectorData('description','Starting timepoint of Each Trial','data',trialsInfo(:,1)); + stop_time_trials = types.core.VectorData('description','Ending timepoint of Each Trial','data',trialsInfo(:,2)); + + id_trials = types.core.ElementIdentifiers('data',int64(0:size(trialsInfo,1)-1)'); + + conditions_trials = cell(size(trialsInfo,1),1); + for iTrial = 1:size(trialsInfo,1) + conditions_trials{iTrial} = runInfo{find(trialsInfo(iTrial,3:4))+2}; + end + + + nwb.intervals_trials = types.core.TimeIntervals('start_time',start_time_trials,'stop_time',stop_time_trials,... + 'colnames',colname,... + 'description',['experimental trials from ' colname 'RunInfo.mat file'],'id',id_trials); + + nwb.intervals_trials.vectordata.set('both_visit', types.core.VectorData('description', 'Both visit condition', 'data', trialsInfo(:,7))); + nwb.intervals_trials.vectordata.set('condition', types.core.VectorData('description', 'Condition Label', 'data', conditions_trials)); + nwb.intervals_trials.vectordata.set('error_run', types.core.VectorData('description', 'Error run', 'data', trialsInfo(:,5))); + nwb.intervals_trials.vectordata.set('stim_run', types.core.VectorData('description', 'Stimulation run condition', 'data', trialsInfo(:,6))); + + disp('Yuta Trial info added..') + + else + disp('No Yuta Trial info present in the selected directory..') + end + + +end \ No newline at end of file diff --git a/io/NWB/YutaMouse41_toNWB.m b/io/NWB/YutaMouse41_toNWB.m deleted file mode 100644 index 0642da23..00000000 --- a/io/NWB/YutaMouse41_toNWB.m +++ /dev/null @@ -1,41 +0,0 @@ -%% Example conversion to NWB of the YutaMouse41-150903 dataset - -% Folder that contains all the files - This is the only input needed -folder_path = 'F:\NWBtoBuzcode\YutaMouse41-150903'; - -%% Start Adding fields to NWB - -% Get info from the xml file -xml = Neuroscope2NWB.GetXMLInfo(folder_path); - -% Add general info to the NWB file -nwb = Neuroscope2NWB.GeneralInfo(xml); - -% Add electrode info -nwb = Neuroscope2NWB.getElectrodeInfo(xml, nwb); - -% Add units info - By default, the spike waveforms are added to the file -nwb = Neuroscope2NWB.getUnitsInfo(xml,nwb); - -% Add stimulation events -nwb = Neuroscope2NWB.getEvents(xml,nwb); - -% Add behavioral info/channels -nwb = Neuroscope2NWB.getBehavior(xml,nwb); - -% Add electrophysiological channels -nwb = Neuroscope2NWB.getElectrophysiology(xml, nwb); - -% Add epochs -nwb = Neuroscope2NWB.getEpochs(xml,nwb); - -% Add trials -nwb = Neuroscope2NWB.getTrials(xml,nwb); - -% Add channels based on the Yuta spreadsheet -nwb = Neuroscope2NWB.special_YutaMouse_recordings(xml,nwb); - - -%% Export to nwb -nwbExport(nwb, 'YutaMouse41_converted.nwb') - diff --git a/io/NWB/addBehavior.m b/io/NWB/addBehavior.m new file mode 100644 index 00000000..801ca3de --- /dev/null +++ b/io/NWB/addBehavior.m @@ -0,0 +1,77 @@ +function nwb = addBehavior(xml,nwb) + %% This function adds Behavior information on an nwb file - COMBINE IT WITH addTrials in order to use bz_LoadBehavior + % It is based on the Buzcode tutorial Behavior file: 20170505_396um_0um_merge.track.behavior.mat + % and the Buzcode wiki: https://github.com/buzsakilab/buzcode/wiki/Data-Formatting-Standards#behavior + + % Konstantinos Nasiotis 2019 + %% Add behavioral data: nwb2.processing.get('behavior').nwbdatainterface + behavior_NWB = types.core.ProcessingModule; + + behavioralFiles = dir([xml.folder_path filesep '*behavior.mat']); + + if length(behavioralFiles) ~= 0 + + for iFile = 1:length(behavioralFiles) + + % The label of the behavior + behavioral_Label = erase(behavioralFiles(iFile).name,'.behavior.mat'); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + behavioral_Label = strsplit(behavioral_Label,'.'); % '20170505_396um_0um_merge' 'track' + behavioral_Label = behavioral_Label{2}; % This section is hardcoded, maybe improve. + % I assumed here that the standardized way of saving behavior files is: experimentName.BehaviorLabel.behavior.mat + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % Load the file + behaviorstruct = load(behavioralFiles(iFile).name); % The example I have has the variable "behavior" saved + behavior = behaviorstruct.behavior; + + + + behavioral_timestamps = behavior.timestamps; + + + behavioral_signals_NWB = types.core.Position; + behavior_field_names = fieldnames(behaviorstruct.behavior)'; + + %% First add the fields that contain channels + iChannelTypeFields = find(ismember(behavior_field_names,{'position','orientation','pupil'})); + + for iField = iChannelTypeFields + + channel_fields = fieldnames(behavior.(behavior_field_names{iField})); + behavioral_channels = zeros(length(behavior.(behavior_field_names{iField}).(channel_fields{1})),length(channel_fields)); + + for iChannel = 1:length(channel_fields) + behavioral_channels(:,iChannel) = behavior.(behavior_field_names{iField}).(channel_fields{iChannel}); + end + + spatial_series = types.core.SpatialSeries('data', behavioral_channels, 'timestamps', behavioral_timestamps,'data_unit', behavior.units); + behavioral_signals_NWB.spatialseries.set(behavior_field_names{iField}, spatial_series); + end + behavior_NWB.nwbdatainterface.set(behavioral_Label,behavioral_signals_NWB); + + + %% Add behaviorinfo field - THIS SHOULDN'T BE SPATIALSERIES - HOWEVER IF IT IS NOT, THE IMPORTING FAILS + behaviorinfoField = find(ismember(behavior_field_names,{'behaviorinfo'})); + + behaviorInfo_fields = fieldnames(behavior.(behavior_field_names{behaviorinfoField})); + + ErrorPerMarkerSignal = behavior.(behavior_field_names{behaviorinfoField}).errorPerMarker; + + spatial_series = types.core.SpatialSeries('data', ErrorPerMarkerSignal, 'timestamps', behavioral_timestamps,'data_unit', behavior.units, ... + 'comments', 'The data field represent the errorPerMarker vector', 'description', behavior.(behavior_field_names{behaviorinfoField}).description, 'control_description', behavior.(behavior_field_names{behaviorinfoField}).acquisitionsystem); + behavioral_signals_NWB.spatialseries.set(behavior_field_names{behaviorinfoField}, spatial_series); + behavior_NWB.nwbdatainterface.set(behavioral_Label,behavioral_signals_NWB); + + end + + behavior_NWB.description = ['Behavioral signals from ' behavioral_Label 'behavior.mat file']; + nwb.processing.set('behavior', behavior_NWB); + disp('Behavioral info added..') + + else + disp('No behavior.mat file present in this directory') + return + end +end \ No newline at end of file diff --git a/io/NWB/addElectrodeInfo.m b/io/NWB/addElectrodeInfo.m new file mode 100644 index 00000000..caecde09 --- /dev/null +++ b/io/NWB/addElectrodeInfo.m @@ -0,0 +1,88 @@ + function nwb = addElectrodeInfo (xml,nwb) + % Adds electrode info in: nwb.general_extracellular_ephys_electrodes + % Konstantinos Nasiotis 2019 + + %% + nShanks = length(xml.spikeDetection.channelGroups.group); + + groups = xml.spikeDetection.channelGroups.group; % Use this for simplicity + + all_shank_channels = cell(nShanks,1); % This will hold the channel numbers that belong in each shank + + % Initialize variables + x = []; + y = []; + z = []; + imp = []; + location = []; + shank = []; + group_name = []; + group_object_view = []; + filtering = []; + shank_channel = []; + amp_channel_id = []; + + device_name = 'implant'; + nwb.general_devices.set(device_name, types.core.Device()); + device_link = types.untyped.SoftLink(['/general/devices/' device_name]); + + for iGroup = 1:nShanks + for iChannel = 1:length(groups{iGroup}.channels.channel) + all_shank_channels{iGroup} = [all_shank_channels{iGroup} str2double(groups{iGroup}.channels.channel{iChannel}.Text)]; + shank_channel = [shank_channel; iChannel-1]; + amp_channel_id = [amp_channel_id; str2double(groups{iGroup}.channels.channel{iChannel}.Text)]; + shank = [shank; iGroup]; + + if nShanks > 9 && iGroup<10 + group_name = [group_name; 'shank' num2str(iGroup) ' ']; + else + group_name = [group_name; 'shank' num2str(iGroup)]; + end + + group_object_view = [group_object_view; types.untyped.ObjectView(['/general/extracellular_ephys/' ['shank' num2str(iGroup)]])]; + + if ~isfield(groups{iGroup}.channels.channel{iChannel},'position') + x = [x; NaN]; + y = [y; NaN]; + z = [z; NaN]; + end + if ~isfield(groups{iGroup}.channels.channel{iChannel},'imp') + imp = [imp; NaN]; + end + if ~isfield(groups{iGroup}.channels.channel{iChannel},'location') + location{end+1,1} = 'unknown'; + end + if ~isfield(groups{iGroup}.channels.channel{iChannel},'filtering') + filtering = [filtering; NaN]; + end + + end + nwb.general_extracellular_ephys.set(['shank' num2str(iGroup)], ... + types.core.ElectrodeGroup( ... + 'description', ['electrode group for shank' num2str(iGroup)], ... + 'location', 'unknown', ... + 'device', device_link)); + + end + + variables = {'x'; 'y'; 'z'; 'imp'; 'location'; 'filtering'; 'group'; 'group_name'; 'shank'; 'shank_channel'; 'amp_channel'}; + + % In order to insert string to a table, they need to be converted to a cell + % first (e.g. location(iElectrode)) + for iElectrode = 1:length(x) + if iElectrode == 1 + tbl = table(x(iElectrode),y(iElectrode),z(iElectrode),imp(iElectrode),{location{iElectrode}},filtering(iElectrode),group_object_view(iElectrode),{group_name(iElectrode,:)},shank(iElectrode),shank_channel(iElectrode),amp_channel_id(iElectrode),... + 'VariableNames', variables); + else + tbl = [tbl; {x(iElectrode),y(iElectrode),z(iElectrode),imp(iElectrode),{location{iElectrode}},filtering(iElectrode),group_object_view(iElectrode),{group_name(iElectrode,:)},shank(iElectrode),shank_channel(iElectrode),amp_channel_id(iElectrode)}]; + end + end + + % add the |DynamicTable| object to the NWB file in + % /general/extracellular_ephys/electrodes + electrode_table = util.table2nwb(tbl, 'metadata about extracellular electrodes'); + nwb.general_extracellular_ephys_electrodes = electrode_table; + + disp('Electrode info added..') + +end \ No newline at end of file diff --git a/io/NWB/addElectrophysiology.m b/io/NWB/addElectrophysiology.m new file mode 100644 index 00000000..468220d9 --- /dev/null +++ b/io/NWB/addElectrophysiology.m @@ -0,0 +1,84 @@ +function nwb = addElectrophysiology(xml,nwb) + %% Adds electrophysiology signals in: nwb.processing.get('ecephys').nwbdatainterface.get('LFP') + % This code was tested on the Buzcode tutorial dataset: 20170505_396um_0um_merge + % and the YutaMouse41-150903 + % Konstantinos Nasiotis 2019 + + %% + [ff basename] = fileparts(xml.folder_path); + + + % bz_LoadBinary + lfpFile = dir([xml.folder_path filesep basename '*.lfp']); + + if length(lfpFile)>1 + disp('More than one .eeg files are present here. No Electrophysiology signals were added') + return + elseif length(lfpFile)==0 + lfpFile = dir([xml.folder_path filesep basename '*.eeg']); + if length(lfpFile)>1 + disp('More than one .lfp files are present here. No Electrophysiology signals were added') + return + elseif length(lfpFile)==0 + disp('No .eeg or .lfp files are present in the selected directory. No Electrophysiology signals were added') + return + end + end + + + % Get the samples number, based on the size of the file + % Check for the precision that samples are saved first + + hdr.nBits = str2double(xml.acquisitionSystem.nBits.Text); + hdr.nChannels = str2double(xml.acquisitionSystem.nChannels.Text); + hdr.sRateOrig = str2double(xml.acquisitionSystem.samplingRate.Text); + hdr.Gain = str2double(xml.acquisitionSystem.amplification.Text); + hdr.sRateLfp = str2double(xml.fieldPotentials.lfpSamplingRate.Text); + + + % Get data type + switch lower(hdr.nBits) + case 16; + hdr.byteSize = 2; + hdr.byteFormat = 'int16'; + case 32; + hdr.byteSize = 4; + hdr.byteFormat = 'int32'; + end + % Guess the number of time points based on the file size + dirInfo = dir(fullfile(lfpFile.folder,lfpFile.name)); + hdr.nSamples = floor(dirInfo.bytes ./ (hdr.nChannels * hdr.byteSize)); + + lfp_data = bz_LoadBinary(fullfile(lfpFile.folder,lfpFile.name), 'duration',Inf, 'frequency',hdr.sRateLfp,'nchannels',hdr.nChannels, 'channels', [1:hdr.nChannels]); % nSamples x 64 + + % If the electrode Information has not already been filled, + % do it now + if isempty(nwb.general_extracellular_ephys_electrodes) + nwb = Neuroscope2NWB.getElectrodeInfo(xml,nwb); + end + + + electrodes_field = types.core.DynamicTableRegion('table',types.untyped.ObjectView('/general/extracellular_ephys/electrodes'),'description','electrode table reference','data',nwb.general_extracellular_ephys_electrodes.id.data); + + lfp = types.core.ElectricalSeries('data', lfp_data', 'electrodes',electrodes_field, 'description', 'lfp signal for all shank electrodes', 'starting_time', 0, 'starting_time_rate', hdr.sRateLfp); + % I TRANSPOSED THE MATRIX HERE TO BE COMPATIBLE WITH THE FUNCTION BZ_GET_LFP + + LFP = types.core.LFP; + LFP.electricalseries.set('lfp',lfp); + + % Check if the ecephys field is already created + if isempty(keys(nwb.processing)) + ecephys = types.core.ProcessingModule('description', ''); + ecephys.description = 'intermediate data from extracellular electrophysiology recordings, e.g., LFP'; + nwb.processing.set('ecephys', ecephys); + else + if ~ismember(keys(nwb.processing),'ecephys') + ecephys = types.core.ProcessingModule('description', ''); + ecephys.description = 'intermediate data from extracellular electrophysiology recordings, e.g., LFP'; + nwb.processing.set('ecephys', ecephys); + end + end + + nwb.processing.get('ecephys').nwbdatainterface.set('LFP', LFP); + disp('Electrophysiological signals added..') +end \ No newline at end of file diff --git a/io/NWB/addEvents.m b/io/NWB/addEvents.m new file mode 100644 index 00000000..632f2aaf --- /dev/null +++ b/io/NWB/addEvents.m @@ -0,0 +1,84 @@ +function nwb = addEvents(xml,nwb) + % Add events: nwb2.stimulus_presentation + % This code was tested on the Buzcode tutorial dataset: 20170505_396um_0um_merge + % Konstantinos Nasiotis 2019 + + % Info about IntervalSeries here: https://nwb-schema.readthedocs.io/en/latest/format.html#sec-intervalseries + + % All the timestamps enter in a vectorized way (even for 2 dimensional + % events with start and stop). + % the start and stop flags are indicated in the .data field by a + % positive (start) or negative (stop) value + + %% + eventFiles = dir([xml.folder_path filesep '*.events.mat']); + + if ~isempty(eventFiles) + all_events = types.core.ProcessingModule(); + + for iFile = 1:length(eventFiles) + events = load(fullfile(eventFiles(iFile).folder,eventFiles(iFile).name)); + + name = fields(events); name = name{1}; % Event name (e.g. ripples) + + % Get the timestamps + timestamps = events.(name).timestamps; + + if size(timestamps,2) == 2 % Start-stop stimulation - The events need to be saved as IntervalSeries + IntervalSeries = types.core.IntervalSeries(); + IntervalSeries.timestamps = [timestamps(:,1) ; timestamps(:,2)]; % I linearize the start and stop timestamps: First event start then event stop + IntervalSeries.data = [ones(size(timestamps,1),1) ; ones(size(timestamps,1),1)* (-1)]; % First event start then event stop: 1 signifies events start, -1 event stop + + all_events.nwbdatainterface.set(name,IntervalSeries); + + + else % Stimulation Onset only - The events need to be saved as AnnotationSeries + AnnotationSeries = types.core.AnnotationSeries('data', repmat({name},length(timestamps),1),'timestamps',timestamps); + all_events.nwbdatainterface.set(name,AnnotationSeries); + + end + + + + +% % % % % % % % % % % % % %% Add the detectorinfo information +% % % % % % % % % % % % % +% % % % % % % % % % % % % % This is no properly set. It simply saves the extra info from +% % % % % % % % % % % % % % the behavior.mat file within the nwb file. +% % % % % % % % % % % % % % Lab members should help in which information should be saved. +% % % % % % % % % % % % % +% % % % % % % % % % % % % detectorinfo = types.core.DynamicTable; +% % % % % % % % % % % % % +% % % % % % % % % % % % % detector_fields = fields(events.(name).detectorinfo); +% % % % % % % % % % % % % for iField = 1:length(detector_fields) +% % % % % % % % % % % % % +% % % % % % % % % % % % % detectorinfo.vectordata.set(detector_fields{iField}, types.core.VectorData('description','Name of the buzcode function that created the events', 'data', events.(name).detectorinfo.(detector_fields{iField}))); +% % % % % % % % % % % % % % +% % % % % % % % % % % % % % IntervalSeries.timestamps = [events.(name).detectorinfo.(detector_fields{iField})(:,1) ; events.(name).detectorinfo.(detector_fields{iField})(:,2)]; % First event start then event stop +% % % % % % % % % % % % % % IntervalSeries.data = [ones(size(events.(name).detectorinfo.(detector_fields{iField}),1),1) ; ones(size(events.(name).detectorinfo.(detector_fields{iField}),1),1)* (-1)]; % First event start then event stop - 1 signifies events start, -1 event stop +% % % % % % % % % % % % % % detectorinfo.set(detector_fields{iField}, IntervalSeries); +% % % % % % % % % % % % % % else +% % % % % % % % % % % % % % detectorinfo.set(detector_fields{iField}, events.(name).detectorinfo.(detector_fields{iField})); +% % % % % % % % % % % % % % end +% % % % % % % % % % % % % end +% % % % % % % % % % % % % detectorinfo.colnames = detector_fields; +% % % % % % % % % % % % % detectorinfo.description = ['Structure fields from ' name ' .behavior.mat file']; +% % % % % % % % % % % % % detectorinfo.id = types.core.ElementIdentifiers('data', 1); +% % % % % % % % % % % % % +% % % % % % % % % % % % % all_events.nwbdatainterface.set([name '_detectorinfo'],detectorinfo); + + + all_events.description = 'Events created from analysis functions'; + nwb.processing.set('events', all_events); + + + end + disp('Events added..') + else + disp('No *.events.mat files found') + end +end + + + + diff --git a/io/NWB/addTrials.m b/io/NWB/addTrials.m new file mode 100644 index 00000000..fe408177 --- /dev/null +++ b/io/NWB/addTrials.m @@ -0,0 +1,152 @@ +function nwb = addTrials(xml,nwb) + + % THIS SECTION ADDS THE TRIAL INFO FROM THE .BEHAVIOR.MAT FILES TO THE NWB FILE - COMBINE IT WITH addBehavior in order to use bz_LoadBehavior + % The trials part of Each file is stored in: nwb.processing.get('behavior').nwbdatainterface.get('trials_BEHAVIORNAME')) + % This is based on the Buzcode tutorial Behavior file: 20170505_396um_0um_merge.track.behavior.mat + % and the Buzcode wiki: https://github.com/buzsakilab/buzcode/wiki/Data-Formatting-Standards#behavior + + % Konstantinos Nasiotis 2019 + %% Fill the fields + behavioralFiles = dir([xml.folder_path filesep '*behavior.mat']); + + if ~isempty(behavioralFiles) + + intervals_trials = types.core.TimeIntervals(); + + for iFile = 1:length(behavioralFiles) + + % The label of the behavior + behavioral_Label = erase(behavioralFiles(iFile).name,'.behavior.mat'); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + behavioral_Label = strsplit(behavioral_Label,'.'); % '20170505_396um_0um_merge' 'track' + behavioral_Label = behavioral_Label{2}; % This section is hardcoded, maybe improve. + % I assumed here that the standardized way of saving behavior files is: experimentName.BehaviorLabel.behavior.mat + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % Load the file + behaviorstruct = load(behavioralFiles(iFile).name); % The example I have has the variable "behavior" saved + behavior = behaviorstruct.behavior; + + all_start_times = types.core.VectorData('description','Starting timepoint of Each Trial','data',behavior.events.trialIntervals(:,1)); + all_stop_times = types.core.VectorData('description','Ending timepoint of Each Trial','data',behavior.events.trialIntervals(:,2)); + + + + %% Add optional events field + if isfield(behavior, 'events') + + trials = behavior.events.trials; + nTrials = length(trials); + + % Find the maximum length of a trial + % All trials will be concatenated on a single matrix + % Smaller trials than the maximum will be filled with + % NaNs + + maxLength = 0; + for iTrial = 1:nTrials + if maxLength1 + disp('More than one spikes.cellinfo file is present. Exiting without adding spiking information') + return + elseif isempty(cellinfoFiles) + disp('No spikes.cellinfo.mat files present in the selected directory.') + return + end + + % If there a single spikes.cellinfo.mat file, load it + load(cellinfoFiles.name) + + %% Add spikes info + % Serialize spiketimes and cluIDs + spike_times = []; + spike_times_index = []; + + current_index = 0; + for iNeuron = 1:length(spikes.UID) + spike_times = [spike_times ; spikes.times{iNeuron}]; + spike_times_index = [spike_times_index; int64(length(spikes.times{iNeuron})+current_index)]; + current_index = spike_times_index(end); + end + + electrode_group = []; + shank_that_neurons_belongs_to = zeros(length(spikes.UID),1); + for iNeuron = 1:length(spikes.UID) + electrode_group = [electrode_group; types.untyped.ObjectView(['/general/extracellular_ephys/' ['shank' num2str(spikes.shankID(iNeuron))]])]; + end + + electrode_group = types.core.VectorData('data', electrode_group, 'description','the electrode group that each spike unit came from'); + + % Initialize the fields needed + spike_times = types.core.VectorData ('data', spike_times, 'description', 'the spike times for each unit'); + spike_times_index = types.core.VectorIndex ('data', spike_times_index, 'target', types.untyped.ObjectView('/units/spike_times')); % The ObjectView links the indices to the spike times + id = types.core.ElementIdentifiers('data', spikes.UID'); + + + waveform_mean = zeros(length(spikes.UID),length(spikes.rawWaveform{1})); + for iNeuron = 1:length(spikes.UID) + waveform_mean(iNeuron,:) = spikes.rawWaveform{iNeuron}; + end + waveform_mean = types.core.VectorData('data', waveform_mean, 'description', 'The mean Waveform for each unit. nNeurons x nSamples'); + + + + %% Fill the units fields + nwb.units = types.core.Units( ... + 'electrode_group', electrode_group, 'electrodes', [], 'electrodes_index', [], 'obs_intervals', [], 'obs_intervals_index', [], ... + 'spike_times', spike_times, 'spike_times_index', spike_times_index, 'waveform_mean', waveform_mean, 'waveform_sd', [], ... + 'colnames', {'shank_id'; 'spike_times'; 'electrode_group'; 'cell_type'; 'global_id'; 'max_electrode'}, ... + 'description', 'Generated from Buzcode2NWB - addUnitsInfo', 'id', id, 'vectorindex', []); + + %% Extra Unit Info + nwb.units.vectordata.set('shankID' , types.core.VectorData('description', 'Which shank each unit belongs to', 'data', spikes.shankID)); + nwb.units.vectordata.set('cluID' , types.core.VectorData('description', 'cluster ID', 'data', spikes.cluID)); + nwb.units.vectordata.set('maxWaveformCh', types.core.VectorData('description', 'The electrode where each unit showed maximum Waveform', 'data', spikes.maxWaveformCh)); + nwb.units.vectordata.set('region' , types.core.VectorData('description', 'The region where each neuron belongs to', 'data', spikes.region)); + + + disp('Spikes info added..') + + +end \ No newline at end of file diff --git a/io/NWB/bz_GetLFP_NWB.m b/io/NWB/bz_GetLFP_NWB.m deleted file mode 100644 index 28e75894..00000000 --- a/io/NWB/bz_GetLFP_NWB.m +++ /dev/null @@ -1,201 +0,0 @@ -function [lfp] = bz_GetLFP_NWB(varargin) -% bz_GetLFP - Get local field potentials. - -% Load local field potentials from NWB file - - -% USAGE -% -% [lfp] = bz_GetLFP(channels,) -% -% INPUTS -% -% channels(required) -must be first input, numeric -% list of channels to load (use keyword 'all' for all) -% channID is 0-indexing, a la neuroscope -% Name-value paired inputs: - -% 'intervals' -list of time intervals [0 10; 20 30] to read from -% the LFP file (default is [0 inf]) -% 'downsample' -factor to downsample the LFP (i.e. 'downsample',5 -% will load a 1250Hz .lfp file at 250Hz) -% -% OUTPUT -% -% lfp struct of lfp data. Can be a single struct or an array -% of structs for different intervals. lfp(1), lfp(2), -% etc for intervals(1,:), intervals(2,:), etc -% .data [Nt x Nd] matrix of the LFP data -% .timestamps [Nt x 1] vector of timestamps to match LFP data -% .interval [1 x 2] vector of start/stop times of LFP interval -% .channels [Nd X 1] vector of channel ID's -% .samplingRate LFP sampling rate [default = 1250] -% .duration duration, in seconds, of LFP interval - -% Copyright (C) 2004-2011 by Michaël Zugaro -% edited by David Tingley, 2017 - - -%% EXAMPLE CALLS -% lfp_NWB = bz_GetLFP_NWB(5); % Loads the entire signal from channel 5 (#6) from the current folder -% lfp_NWB = bz_GetLFP_NWB(5, 'nwb_file', nwb_file); % Loads the entire signal from channel 5 (#6) from the current folder from the specified NWB file -% lfp_NWB = bz_GetLFP_NWB([5 10], 'nwb_file', nwb_file); % Loads channels 5 (#6) and 10 (#11) -% lfp_NWB = bz_GetLFP_NWB(5, 'restrict',[0 120]); % Loads channel 5 (#6) between [0, 120] seconds -% lfp_NWB = bz_GetLFP_NWB(5, 'restrict',[0 120;240.2 265.23]); % Loads channel 5 (#6) between [0, 120] and [240.2 265.23] seconds -% lfp_NWB = bz_GetLFP_NWB(5, 'restrict',[0 120], 'downsample', 2); % Loads channel 5 (#6) between [0, 120] seconds and downsamples by a factor of 2 - - -% Added support for NWB: Konstantinos Nasiotis 2019 - - -%% Parse the inputs! - -channelsValidation = @(x) isnumeric(x) || strcmp(x,'all'); - -% parse args -p = inputParser; -addParameter(p,'nwb_file','',@isstr); -addRequired(p,'channels',channelsValidation) -addParameter(p,'basename','',@isstr) -addParameter(p,'intervals',[],@isnumeric) -addParameter(p,'restrict',[],@isnumeric) -addParameter(p,'basepath',pwd,@isstr); -addParameter(p,'downsample',1,@isnumeric); -addParameter(p,'saveMat',false,@islogical); -addParameter(p,'forceReload',false,@islogical); -addParameter(p,'noPrompts',false,@islogical); - -parse(p,varargin{:}) -nwb_file = p.Results.nwb_file; -channels = p.Results.channels; -downsamplefactor = p.Results.downsample; -basepath = p.Results.basepath; -noPrompts = p.Results.noPrompts; - -% doing this so you can use either 'intervals' or 'restrict' as parameters to do the same thing -intervals = p.Results.intervals; -restrict = p.Results.restrict; -if isempty(intervals) && isempty(restrict) % both empty - intervals = [0 Inf]; -elseif isempty(intervals) && ~isempty(restrict) % intervals empty, restrict isn't - intervals = restrict; -end - -[filepath,name,ext] = fileparts(nwb_file); - - -%% let's check that there is an appropriate NWB file -if isempty(nwb_file) - %disp('No nwb_file given, so we look for a *nwb file in the current directory') - d = dir('*nwb'); - if length(d) > 1 % we assume one .nwb file or this should break - error('there is more than one .nwb file in this directory?'); - elseif length(d) == 0 - d = dir('*nwb'); - if isempty(d) - error('could not find an nwb file..') - end - end - nwb_file = fullfile(d.folder, d.name); - - lfp.Filename = d.name; -else - lfp.Filename = [name ext]; -end - -%% things we can parse from sessionInfo or xml file - -sessionInfo = bz_getSessionInfo(basepath, 'noPrompts', noPrompts); - -try - samplingRate = sessionInfo.lfpSampleRate; -catch - samplingRate = sessionInfo.rates.lfp; % old ugliness we need to get rid of -end -samplingRateLFP = samplingRate./downsamplefactor; - -if mod(samplingRateLFP,1)~=0 - error('samplingRate/downsamplefactor must be an integer') -end -%% Channel load options -%Right now this assumes that all means channels 0:nunchannels-1 (neuroscope -%indexing), we could also add options for this to be select region or spike -%group from the xml... -if strcmp(channels,'all') - channels = sessionInfo.channels; -end - -%% get the data -disp('loading LFP file...') -nIntervals = size(intervals,1); - - -nwb2 = nwbRead(nwb_file); - - -try - % Check if the data is in LFP format - all_lfp_keys = keys(nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries); - - for iKey = 1:length(all_lfp_keys) - if ismember(all_lfp_keys{iKey}, {'lfp','bla bla bla'}) %%%%%%%% ADD MORE HERE, DON'T KNOW WHAT THE STANDARD FORMATS ARE - iLFPDataKey = iKey; - LFPDataPresent = 1; - break % Once you find the data don't look for other keys/trouble - else - LFPDataPresent = 0; - end - end -catch - LFPDataPresent = 0; -end - -if ~LFPDataPresent - error ('No LFP signals exist in this .nwb file') -end - - - - - -% returns lfp/bz format -for i = 1:nIntervals - lfp(i).duration = (intervals(i,2)-intervals(i,1)); - lfp(i).interval = [intervals(i,1) intervals(i,2)]; - - % Load data and put into struct - % we assume 0-indexing like neuroscope, but bz_LoadBinary uses 1-indexing to - % load.... - - if intervals(i,2) == Inf - use_this_bracket = sessionInfo.samples_NWB/sessionInfo.lfpSampleRate; - data_temp = zeros(length(channels),ceil((use_this_bracket - intervals(i,1))*samplingRateLFP)); - else - use_this_bracket = intervals(i,2); - data_temp = zeros(length(channels),floor((intervals(i,2) - intervals(i,1))*samplingRateLFP)); - end - - % This is not optimized yet - - for iChannel = 1:length(channels) - disp(['Now concatenating channel: ' num2str(channels(iChannel))]) - data_temp(iChannel,:) = nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries.get(all_lfp_keys{iLFPDataKey}).data.load([channels(iChannel)+1, intervals(i,1)*samplingRate+1],[1, downsamplefactor],[channels(iChannel)+1, use_this_bracket*samplingRate]); - end - - lfp(i).data = data_temp'; clear data_temp - lfp(i).timestamps = [lfp(i).interval(1):(1/samplingRateLFP):... - (lfp(i).interval(1)+(length(lfp(i).data)-1)/... - samplingRateLFP)]'; - lfp(i).channels = channels; - lfp(i).samplingRate = samplingRateLFP; - % check if duration is inf, and reset to actual duration... - if lfp(i).interval(2) == inf - lfp(i).interval(2) = length(lfp(i).timestamps)/lfp(i).samplingRate; - lfp(i).duration = (lfp(i).interval(i,2)-lfp(i).interval(i,1)); - end - - if isfield(sessionInfo,'region') && isfield(sessionInfo,'channels') - [~,~,regionidx] = intersect(lfp(i).channels,sessionInfo.channels,'stable'); - lfp(i).region = sessionInfo.region(regionidx); % match region order to channel order.. - end -end diff --git a/io/NWB/bz_GetSpikes_NWB.m b/io/NWB/bz_GetSpikes_NWB.m deleted file mode 100644 index 525674b3..00000000 --- a/io/NWB/bz_GetSpikes_NWB.m +++ /dev/null @@ -1,245 +0,0 @@ -function spikes = bz_GetSpikes_NWB(varargin) -% bz_getSpikes - Get spike timestamps. -% -% USAGE -% -% spikes = bz_getSpikes(varargin) -% -% INPUTS -% -% spikeGroups -vector subset of shank IDs to load (Default: all) -% region -string region ID to load neurons from specific region -% (requires sessionInfo file or units->structures in xml) -% UID -vector subset of UID's to load -% getWaveforms -logical (default=true) to load mean of raw waveform data -% forceReload -logical (default=false) to force loading from -% res/clu/spk files -% saveMat -logical (default=false) to save in buzcode format -% noPrompts -logical (default=false) to supress any user prompts - - -% spikes_NWB = bz_GetSpikes_NWB('UID',[2 4 7]); % Loads neurons [2, 4 7] from an NWB at the current directory -% spikes_NWB = bz_GetSpikes_NWB('nwb_file', nwb_file, 'UID',[2 4 7]); % Loads neurons [2, 4 7] from a specific NWB file -% spikes_NWB = bz_GetSpikes_NWB('nwb_file', nwb_file, 'spikeGroups', [2,4]); % Loads the neurons from a specific NWB file and spikeGroups [2,4] -% spikes_NWB = bz_GetSpikes_NWB('nwb_file', nwb_file, 'region', 'unknown'); % Loads the neurons from a specific NWB file and a specified region -% spikes_NWB = bz_GetSpikes_NWB('nwb_file', nwb_file, 'saveMat',true); % saves a cellInfo.mat file with all the Neurons - - -% Konstantinos Nasiotis 2019 - - - -%% TODO -disp(' 1. The template waveform should be filled by: nwb2.units.waveform_mean The example dataset didnt have any. I added a template') -disp(' 2. The maxWaveformCh should be added on the example dataset') -disp(' 3. CluID is probably filled by Kilosort - should be added on the example dataset') -disp(' 4. getWaveforms is not supported. Not sure nwb holds the spiking waveforms - Havent seen a field - UPDATE, THIS IS NOW SUPPORTED - CHECK nwb2.processing.get("ecephys").nwbdatainterface.get("SpikeEventSeries1")') - - - -%% Deal With Inputs -spikeGroupsValidation = @(x) assert(isnumeric(x) || strcmp(x,'all'),... - 'spikeGroups must be numeric or "all"'); - -p = inputParser; -addParameter(p,'nwb_file','',@isstr); -addParameter(p,'spikeGroups','all',spikeGroupsValidation); -addParameter(p,'region','',@isstr); % won't work without sessionInfodata -addParameter(p,'UID',[],@isvector); -addParameter(p,'basepath',pwd,@isstr); -addParameter(p,'getWaveforms',true,@islogical) -addParameter(p,'forceReload',false,@islogical); -addParameter(p,'saveMat',false,@islogical); -addParameter(p,'noPrompts',false,@islogical); - -parse(p,varargin{:}) - -%% Adding support only for spikeGroups, region and UID for now - -nwb_file = p.Results.nwb_file; -spikeGroups = p.Results.spikeGroups; -region = p.Results.region; -UID = p.Results.UID; -% basepath = p.Results.basepath; -% getWaveforms = p.Results.getWaveforms; -% forceReload = p.Results.forceReload; -saveMat = p.Results.saveMat; -% noPrompts = p.Results.noPrompts; - - - - -%% let's check that there is an appropriate NWB file -if isempty(nwb_file) - %disp('No nwb_file given, so we look for a *nwb file in the current directory') - d = dir('*nwb'); - if length(d) > 1 % we assume one .nwb file or this should break - error('there is more than one .nwb file in this directory?'); - elseif length(d) == 0 - d = dir('*nwb'); - if isempty(d) - error('could not find an nwb file..') - end - end - nwb_file = fullfile(d.folder, d.name); -end - - - -nwb2 = nwbRead(nwb_file); -[the_path, ~, ~] = fileparts(nwb_file); - - -% load([new_path_for_files filesep name '.sessionInfo.mat']) -sessionInfo = get_SessionInfoFromNWB(nwb2); - - - nNeurons = length(nwb2.units.id.data.load); - - -% % % % % % % % % % % % % % % % % % % % % % -% CHEK THIS nwb2.units.electrode_group.data(1).refresh(nwb2) -% % % % % % % % % % % % % % % % % % % % - - - % Assign neuron to region based on the region that its Shank belongs to - - % Get a single electrode that belongs in that shank - electrodes2Shank = nwb2.general_extracellular_ephys_electrodes.vectordata.get('group_name').data; - electrodes2Region = nwb2.general_extracellular_ephys_electrodes.vectordata.get('location').data; - - neurons2ShankNames = cell(1,nNeurons); - all_region = cell(1,nNeurons); - - for iNeuron = 1:length(nwb2.units.electrode_group.data) - path = nwb2.units.electrode_group.data(iNeuron).path; - inds = strfind(path, '/'); - neurons2ShankNames{iNeuron} = path(inds(end)+1:end); - - ElectrodesInThatShank = find(strcmp(electrodes2Shank,neurons2ShankNames{iNeuron})); - all_region{iNeuron} = electrodes2Region{ElectrodesInThatShank(1)}; - - end - - - uniqueShanks = unique(neurons2ShankNames,'legacy'); - shankID_of_selected_Neurons = zeros(1, nNeurons); - for iNeuron = 1:nNeurons - shankID_of_selected_Neurons(iNeuron) = find(strcmp(uniqueShanks, neurons2ShankNames{iNeuron})); - end - - - %% Make the check here of what will be loaded - - selected_neurons_UID = false(1,nNeurons); - selected_neurons_spikeGroups = false(1,nNeurons); - selected_neurons_region = false(1,nNeurons); - - % Check which neurons were selected - if ~isempty(UID) - selected_neurons_UID(UID) = true; - end - %Check which Shanks were selected - if ~isempty(spikeGroups) - for iShank = spikeGroups - selected_neurons_spikeGroups(find(shankID_of_selected_Neurons==iShank)) = true; - end - end - %Check which regions were selected - if ~isempty(region) - selected_neurons_region(find(contains(all_region,region))) = true; - end - - UID = find(selected_neurons_UID | selected_neurons_spikeGroups | selected_neurons_region); - - if isempty(UID) - warning(['Warning: The selection made didnt specify any subgroup of neurons. Loading all neurons!']) - UID = 1:nNeurons; - end - - - %% The first time that this function is loaded, make sure to save a .spikes.cellInfo.mat that contains information from all neurons - - if saveMat - UID_forSaveMAt = 1:nNeurons; - temp = get_the_spikes_from_selected_UIDs(UID_forSaveMAt, nwb2, sessionInfo, saveMat, all_region, neurons2ShankNames, the_path); clear temp - end - - spikes = get_the_spikes_from_selected_UIDs(UID, nwb2, sessionInfo, 0, all_region, neurons2ShankNames, the_path); - -end - - - - -function spikes = get_the_spikes_from_selected_UIDs(UID, nwb2, sessionInfo, saveMat, all_region, neurons2ShankNames, the_path) - - %% Get Spikes - spikes = struct; - spikes.samplingRate = sessionInfo.rates.wideband; - spikes.UID = UID; - - % The template waveform should be filled by: nwb2.units.waveform_mean - - template_Waveform = [21.2963012297339,21.2206998551635,21.5609060407305,... % This is from a template I used on Kilosort. - 22.8392565561944,28.2241362812803,30.1107342194246,... % I assign the same on every neuron. - 23.9595314702837,24.3822118826549,23.4681225355758,... % Check if I can get this from nwb - 18.0866792366068,11.9973321575690,-2.96486715514583,... - -110.074832790885,-186.628097395696,-240.085142069235,... - -183.947685024562,-143.562805299476,-104.992358564081,... - -3.29132763624548,22.3478476214865,44.7121087898714,... - 73.1141706455415,79.3615933259538,82.9182943568817,... - 75.0076414359195,70.1691534634109,65.2413184118645,... - 51.9698407486342,45.6296345630672,41.3822118826549,... - 37.1519713328267,31.5334146317958]; - - times = cell(1,length(UID)); - rawWaveform = cell(1,length(UID)); - spindices = []; - - entry = 0; - for iNeuron = UID - entry = entry+1; - if iNeuron == 1 - times_temp = nwb2.units.spike_times.data.load(1:sum(nwb2.units.spike_times_index.data.load(iNeuron))); - else - times_temp = nwb2.units.spike_times.data.load(sum(nwb2.units.spike_times_index.data.load(iNeuron-1))+1:sum(nwb2.units.spike_times_index.data.load(iNeuron))); - end - times{entry} = times_temp(times_temp~=0)'; - - rawWaveform{entry} = template_Waveform; % The template waveform should be filled by: nwb2.units.waveform_mean The example dataset didn't have any - spindices = [spindices ; times{entry} ones(length(times{entry}),1)*iNeuron]; - end - - % Spindices have to be sorted according to when each spike occured - [~,sortedIndices] = sort(spindices(:,1)); - spindices = spindices(sortedIndices,:); - - - uniqueShanks = unique(neurons2ShankNames,'legacy'); - shankID_of_selected_Neurons = zeros(1, length(UID)); - for iNeuron = 1:length(UID) - shankID_of_selected_Neurons(iNeuron) = find(strcmp(uniqueShanks, neurons2ShankNames{UID(iNeuron)})); - end - - - spikes.times = times; - spikes.shankID = shankID_of_selected_Neurons; - spikes.cluID = ones(1,length(UID))*(-1); % THESE ARE THE SPIKING TEMPLATES. THEY ARE FILLED FROM KILOSORT. I ADD A NEGATIVE VALUE TO SEE IF IT CAUSES AN ERROR SOMEWHERE - spikes.rawWaveform = rawWaveform; - spikes.maxWaveformCh = ones(1,length(UID))*(-1); % THESE ASSIGN THE MAXIMUM WAVEFORM TO A ACHANNEL. CHECK HOW TO ADD THIS. I ADD A NEGATIVE VALUE TO SEE IF IT CAUSES AN ERROR SOMEWHERE - spikes.sessionName = sessionInfo.FileName; - spikes.numcells = length(UID); - spikes.spindices = spindices; % This holds the timing of each spike, sorted, and the neuron it belongs to. - - spikes.region = all_region(UID); - - if saveMat - save(fullfile(the_path, [sessionInfo.FileName '.spikes.cellinfo.mat']),'spikes') - end -end - - - - - diff --git a/io/NWB/bz_LoadBehavior_NWB.m b/io/NWB/bz_LoadBehavior_NWB.m deleted file mode 100644 index bc48115a..00000000 --- a/io/NWB/bz_LoadBehavior_NWB.m +++ /dev/null @@ -1,436 +0,0 @@ -function behavior = bz_LoadBehavior_NWB(varargin) - -%% Load Behavior from NWB files following the Buzcode format -% behavior = bz_LoadBehavior_NWB; % Checks the current directory for a NWB file. Gives a list with the available behavior fields in the NWB file for the user to choose which one to load -% behavior = bz_LoadBehavior_NWB('nwb_file', nwb_file); % Checks the specified NWB file. Gives a list with the available behavior fields in the NWB file for the user to choose which one to load -% behavior = bz_LoadBehavior_NWB('nwb_file', nwb_file, 'EightMazePosition_norm_spatial_series'); % Loads the specific behavior from the specified NWB file -% behavior = bz_LoadBehavior_NWB('EightMazePosition_norm_spatial_series'); % Loads the specific behavior from a NWB file in the current directory - - -% The behaviors are loaded straight from the NWB file. - -% This code creates a file with the behavior selected the first time it's used. -% This is done in order to speed up the process of calling it later -% (.map takes a long time to compute if the trials are long). - - -% This code is based on the bz_LoadBehavior -% Konstantinos Nasiotis 2019 - - -%% TODO - -disp('1. Check if the BAD trials should be included in the computation of the .map field') - - - - -%% - p = inputParser; - addParameter(p,'nwb_file','',@isstr); - - if size(varargin,2)>2 - behaviorName = varargin{3}; - varargin = {varargin{1,1:2}}; - parse(p,varargin{:}) - nwb_file = p.Results.nwb_file; - elseif size(varargin,2)==1 - nwb_file = []; - behaviorName = varargin{1}; - elseif size(varargin,2)==2 - parse(p,varargin{:}) - nwb_file = p.Results.nwb_file; - behaviorName = []; - elseif size(varargin,2)==0 - nwb_file = []; - behaviorName = []; - end - - - %% let's check that there is an appropriate NWB file - if isempty(nwb_file) - %disp('No nwb_file given, so we look for a *nwb file in the current directory') - d = dir('*nwb'); - if length(d) > 1 % we assume one .nwb file or this should break - error('there is more than one .nwb file in this directory?'); - elseif length(d) == 0 - d = dir('*nwb'); - if isempty(d) - error('could not find an NWB file..') - end - end - nwb_file = fullfile(d.folder, d.name); - end - - - nwb2 = nwbRead(nwb_file); - - % Check if behavior fields exists in the dataset - try - allBehaviorKeys = keys(nwb2.processing.get('behavior').nwbdatainterface); - - behavior_exist_here = ~isempty(allBehaviorKeys); - if ~behavior_exist_here - disp('No behavior in this .nwb file') - return - - else - disp(' ') - disp('The following behavior types are present in this dataset') - disp('------------------------------------------------') - for iBehavior = 1:length(allBehaviorKeys) - disp(allBehaviorKeys{iBehavior}) - end - disp(' ') - end - catch - disp('No behavior in this .nwb file') - return - end - - %% Check if the behavior file has already been computed (it takes some time to fill the .mapping field, especially when the trials are long) - % If the behavior file is already computed, just load the file and return - - % FIRST CHECK - CHECK IF THE behaviorName REQUESTED EXISTS - - [the_path, name, ext] = fileparts(nwb_file); - - if exist(fullfile(the_path, [name '.' behaviorName '.behavior.mat']) , 'file') == 2 - load(fullfile(the_path, [name '.' behaviorName '.behavior.mat'])) - disp('Just loading the already computed behavior struct') - return - else - disp('A file with the specified Behavior hasnt already been created. Select the behavior from the list') - end - - %% If a the behavior file doesn't exist, create a new one - - % Check if a specific behavior was called to be loaded. If not, display a - % pop-up list with the available behaviors for selection - [iBehavior, ~] = listdlg('PromptString','Which behavior type would you like to load?',... - 'ListString',allBehaviorKeys,'SelectionMode','single'); - - % Fill the behavior - behavior = struct; % Initialize - - - % Select the key that is within the Behavior selection - allBehaviorSUBKeys = keys(nwb2.processing.get('behavior').nwbdatainterface.get(allBehaviorKeys(iBehavior)).spatialseries); - if length(allBehaviorSUBKeys)>1 - [iSubKeyBehavior, ~] = listdlg('PromptString','Multiple Behavior channel-selections within the selected behavior',... - 'ListString',allBehaviorSUBKeys,'SelectionMode','single'); - else - iSubKeyBehavior = 1; - end - - %% Check if the behavior file has already been computed (it takes some time to fill the .mapping field, especially when the trials are long) - % If the behavior file is already computed, just load the file and return - - % SECOND CHECK - HERE A SUBKEY HAS ALREADY BEEN SET FROM THE SELECTION - - if exist(fullfile(the_path, [name '.' allBehaviorSUBKeys{iSubKeyBehavior} '.behavior.mat']) , 'file') == 2 - load(fullfile(the_path, [name '.' allBehaviorSUBKeys{iSubKeyBehavior} '.behavior.mat'])) - disp('Just loading the already computed behavior struct') - return - end - - - %% Get the behavior info - disp('computing') - - selected_behavior = nwb2.processing.get('behavior').nwbdatainterface.get(allBehaviorKeys(iBehavior)).spatialseries.get(allBehaviorSUBKeys{iSubKeyBehavior}); % This is for easier reference through the code - - behavior.samplingRate = 1000 ./ nanmean(diff(selected_behavior.timestamps.load))./1000; - behavior.units = selected_behavior.data_unit; - - Info = allBehaviorKeys{iBehavior}; -% behavior.description = selected_behavior.description; - - % Check if timestamps and data have values. If not something weird - % is going on - if ~isempty(selected_behavior.timestamps) - behavior.timestamps = selected_behavior.timestamps.load; - else - behavior.timestamps = []; - warning(['Behavior: ' Info ' --- Timestamps are empty: weird']) - end - if ~isempty(selected_behavior.data) - - % Specifically for position signals, seperate to x and y - if ~isempty(strfind(allBehaviorKeys(iBehavior),'position')) - data = selected_behavior.data.load; - datasubstruct.x = data(1,:)'; - if size(data,1)>1 - datasubstruct.y = data(2,:)'; - end - if size(data,1)>2 - datasubstruct.z = data(3,:)'; - end - else - datasubstruct.data = selected_behavior.data.load; - - end - - - else - datasubstruct.data = []; - warning(['Behavior: ' Info ' --- Data is empty: weird']) - end - - - % position: .x, .y, and .z - % units: millimeters, centimeters, meters[default] - % orientation: .x, .y, .z, and .w - % rotationType: euler or quaternion[default] - % pupil: .x, .y, .diameter - - -% datasubstruct.reference_frame = nwb2.acquisition.get(all_raw_keys(allBehaviorKeys(iBehavior))).reference_frame; This seems to be present only on the position_sensor channels - - % Conside removing these - datasubstruct.starting_time_unit = selected_behavior.starting_time_unit; - datasubstruct.timestamps_interval = selected_behavior.timestamps_interval; - datasubstruct.comments = selected_behavior.comments; - datasubstruct.control = selected_behavior.control; - datasubstruct.control_description = selected_behavior.control_description; - datasubstruct.data_resolution = selected_behavior.data_resolution; - datasubstruct.starting_time = selected_behavior.starting_time; - datasubstruct.help = selected_behavior.help; - - - - %% Exclude special characters from the behavior name - BehaviorLabel = allBehaviorKeys{iBehavior}; - - clean_string = regexprep(BehaviorLabel,'[^a-zA-Z0-9_]',''); - if ~strcmp(clean_string, BehaviorLabel) - disp(['The variable name (' BehaviorLabel ') of the Behavior was changed to exclude special characters']) - BehaviorLabel = clean_string; - end - - %% If this is a position signal, rename to position to have compatibility with the Buzsaki functions - if ~isempty(strfind(allBehaviorKeys(iBehavior),'position')) - BehaviorLabel = 'position'; - end - - - behavior.(BehaviorLabel) = datasubstruct; - - % BehaviorInfo - behaviorinfo.description = selected_behavior.description; - behaviorinfo.acquisitionsystem = 'Fill Me'; - behaviorinfo.substructnames = {BehaviorLabel}; - behavior.behaviorinfo = behaviorinfo; - - - - %% Fill the events substructure - - nTrials = length(nwb2.intervals_trials.start_time.data.load); - uniqueConditions = unique(nwb2.intervals_trials.vectordata.get('condition').data); - - - all_conditions = nwb2.intervals_trials.vectordata.get('condition').data; - for iCondition = 1:length(all_conditions) - - events.trialConditions(iCondition) = find(strcmp(all_conditions{iCondition}, uniqueConditions));%1x221 double which unique condition each trialCondition belongs to - INDEX - - end - events.trialIntervals = [nwb2.intervals_trials.start_time.data.load nwb2.intervals_trials.stop_time.data.load]; % 221x2 double start-end of trial in seconds - events.conditionType = unique(nwb2.intervals_trials.vectordata.get('condition').data)'; % 1x10 cell (central, wheel ...) - -% nwb2.intervals_trials.id.data.load -% nwb2.intervals_trials.colnames -% nwb2.intervals_trials.vectordata -% nwb2.intervals_trials.vectordata.get('both_visit').data.load -% nwb2.intervals_trials.vectordata.get('condition').data -% nwb2.intervals_trials.vectordata.get('error_run').data.load -% nwb2.intervals_trials.vectordata.get('stim_run').data.load -% nwb2.intervals_trials.vectorindex - - trialTimestampBounds = [nwb2.intervals_trials.start_time.data.load nwb2.intervals_trials.stop_time.data.load]; - position_timestamps = selected_behavior.timestamps.load; - - for iTrial = 1:nTrials - - %% Load only the samples from each trial - % Find the indices of the timestamps that are selected - - [~, iPositionTimestamps] = histc(trialTimestampBounds(iTrial,:), position_timestamps); - - if iPositionTimestamps(1)~=0 && iPositionTimestamps(2)==0 - iPositionTimestamps(2) = length(position_timestamps); - end - - try - position = selected_behavior.data.load([1, iPositionTimestamps(1)], [Inf, iPositionTimestamps(2)]); - catch - error ('Error loading selected behavior file. Are there datapoints within the trials selection time-period?') - end - - % The boundaries struct will hold the position boundaries of each - % TYPE of CONDITION, for all the trials on each condition Type - if ~exist('boundaries') - if size(position,1)==1 - boundaries = repmat(struct('x_min',0,'x_max',0),length(uniqueConditions),1); - elseif size(position,1)==2 - boundaries = repmat(struct('x_min',0,'x_max',0,'y_min',0,'y_max',0),length(uniqueConditions),1); - elseif size(position,1)==3 - boundaries = repmat(struct('x_min',0,'x_max',0,'y_min',0,'y_max',0,'z_min',0,'z_max',0),length(uniqueConditions),1); - end - end - - -% clr = rand(1,3); -% plot(position(1,:),position(2,:),'.','color',clr); -% drawnow - - % Get the index of the condition that the trials belong to - this will - % be used for defining the spatial boundaries of the condition - iConditionOfTrial = find(strcmp(uniqueConditions, nwb2.intervals_trials.vectordata.get('condition').data(iTrial)))'; - - - trials{1,iTrial}.x = position(1,:)';% 608 x 1 - if size(position,1)>1 - trials{1,iTrial}.y = position(2,:)';% 608 x 1 - end - if size(position,1)>2 - trials{1,iTrial}.z = position(3,:)';% 608 x 1 - end - - trials{1,iTrial}.timestamps = position_timestamps(iPositionTimestamps(1):iPositionTimestamps(2));% 608 x 1 - trials{1,iTrial}.direction = 'FILL ME'; % 'clockwise' 'counterclockwise' - trials{1,iTrial}.type = nwb2.intervals_trials.vectordata.get('condition').data{iTrial}; % 'central alternation' - % trials{1,iTrial}.orientation = % 1x1 struct - % trials{1,iTrial}.errorPerMarker = 608 x 1 - - end - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % Map the trials to a 1x201 vector - bins = 200; - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - % normalize positions to template - c=1; - - uniqueConditions = unique((events.trialConditions)); - - for iCondition = 1:length(uniqueConditions) - - map{iCondition}=[]; - t_conc=[]; - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % MAYBE SELECT ONLY THE TRIALS THAT WEREN'T BAD HERE -% iTrialsInCondition = find(events.trialConditions == uniqueConditions(iCondition) & ~nwb2.intervals_trials.vectordata.get('error_run').data.load'); - iTrialsInCondition = find(events.trialConditions == uniqueConditions(iCondition)); -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - - - for iTrial = 1:length(iTrialsInCondition) - selected_trial = trials{iTrialsInCondition(iTrial)}; % I set this for faster typing - - if size(position,1)==1 - t_conc = [selected_trial.timestamps, selected_trial.x, 20*(selected_trial.timestamps - selected_trial.timestamps(1))]; % Check what this 20 is - elseif size(position,1)==2 - t_conc = [selected_trial.timestamps, selected_trial.x selected_trial.y, 20*(selected_trial.timestamps - selected_trial.timestamps(1))]; % Check what this 20 is - elseif size(position,1)==3 - t_conc = [selected_trial.timestamps, selected_trial.x selected_trial.y selected_trial.z, 20*(selected_trial.timestamps - selected_trial.timestamps(1))]; % Check what this 20 is - end - -% % % % % % templateVector = 1:nBins; -% % % % % % target = size(position,2); -% % % % % % n_ent = length(templateVector); -% % % % % % trials{1,iTrial}.mapping = round(interp1( 1:n_ent, templateVector, linspace(1, n_ent, target) ))'; - - disp('the computation below doesnt make the usage of this function practical for on the fly retrieval of the Behavior') - if length(t_conc)>=bins - while length(t_conc)>bins+1 - - di = pdist(t_conc); - s = squareform(di); - s(find(eye(size(s))))=nan; - [a b] = min(s(:)); - [coords blah] = find(s==a); - t_conc(coords(1),:) = (t_conc(coords(1),:)+t_conc(coords(2),:))./2; - t_conc(coords(2),:) = []; - end - t_conc_all(iTrial,:,:) = t_conc; - end - end - if length(iTrialsInCondition)>0 - map{iCondition} = squeeze(median(t_conc_all(:,:,:),1)); - end - - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % Assign to the right behavior field - if size(position,1)==1 % 1D - events.map{iCondition}.x = map{iCondition}(:,2); - elseif size(position,1)==2 % 2D - events.map{iCondition}.x = map{iCondition}(:,2); - events.map{iCondition}.y = map{iCondition}(:,3); - elseif size(position,1)==3 % 3D - events.map{iCondition}.x = map{iCondition}(:,2); - events.map{iCondition}.y = map{iCondition}(:,3); - events.map{iCondition}.z = map{iCondition}(:,4); - end - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - clear t_conc_all - - disp('finding mapping...') - for iTrial =1:length(iTrialsInCondition) % all trial types (rotations) - selected_trial = trials{iTrialsInCondition(iTrial)}; % I set this for faster typing - for iPoint = 1:length(selected_trial.timestamps) - - if size(position,1)==1 % 1D - [a b] = min(nansum(abs([selected_trial.timestamps(iPoint)-map{iCondition}(:,1),... % Timestamp - selected_trial.x(iPoint)-map{iCondition}(:,2),... % X_COORDINATE - (selected_trial.timestamps(iPoint)-selected_trial.timestamps(1))*50-map{iCondition}(:,1),... % penalty for time differences - (40*(iPoint./length(selected_trial.timestamps)*length(map{iCondition}) - (1:length(map{iCondition})))')])')); % penalty for order differences - - elseif size(position,1)==2 % 2D - [a b] = min(nansum(abs([selected_trial.timestamps(iPoint)-map{iCondition}(:,1),... % Timestamp - selected_trial.x(iPoint)-map{iCondition}(:,2),... % X_COORDINATE - selected_trial.y(iPoint)-map{iCondition}(:,3),... % Y_COORDINATE - (selected_trial.timestamps(iPoint)-selected_trial.timestamps(1))*50-map{iCondition}(:,1),... % penalty for time differences - (40*(iPoint./length(selected_trial.timestamps)*length(map{iCondition}) - (1:length(map{iCondition})))')])')); % penalty for order differences - - elseif size(position,1)==3 % 3D - [a b] = min(nansum(abs([selected_trial.timestamps(iPoint)-map{iCondition}(:,1),... % Timestamp - selected_trial.x(iPoint)-map{iCondition}(:,2),... % X_COORDINATE - selected_trial.y(iPoint)-map{iCondition}(:,3),... % Y_COORDINATE - selected_trial.z(iPoint)-map{iCondition}(:,4),... % Z_COORDINATE - (selected_trial.timestamps(iPoint)-selected_trial.timestamps(1))*50-map{iCondition}(:,1),... % penalty for time differences - (40*(iPoint./length(selected_trial.timestamps)*length(map{iCondition}) - (1:length(map{iCondition})))')])')); % penalty for order differences - end - - mapping{iCondition}{iTrial}(iPoint,:) = [map{iCondition}(b,1:end) b selected_trial.timestamps(iPoint)]; - - end - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % Assign to the right behavior field - trials{iTrialsInCondition(iTrial)}.mapping = mapping{iCondition}{iTrial}(:,end-1); - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - end - end - - events.trials = trials; - behavior.events = events; - - %% Check that the behavior structure meets buzcode standards - [isBehavior] = bz_isBehavior(behavior); - switch isBehavior - case false - warning('Your behavior structure does not meet buzcode standards. Sad.') - end - - %% Save the behavior files if they are not already saved - if exist(fullfile(the_path, [name '.' allBehaviorSUBKeys{iSubKeyBehavior} '.behavior.mat']) ,'file') ~=2 - save(fullfile(the_path, [name '.' allBehaviorSUBKeys{iSubKeyBehavior} '.behavior.mat']), 'behavior') - end - -end \ No newline at end of file diff --git a/io/NWB/bz_LoadEvents_NWB.m b/io/NWB/bz_LoadEvents_NWB.m deleted file mode 100644 index 76afc640..00000000 --- a/io/NWB/bz_LoadEvents_NWB.m +++ /dev/null @@ -1,121 +0,0 @@ -function events = bz_LoadEvents_NWB(varargin) - -%% Load events from NWB files following the Buzcode format - -% Example calls: -% events = bz_LoadEvents_NWB; % Checks the current directory for a NWB file. Gives a list with the available events in the NWB file for the user to choose which one to load -% events = bz_LoadEvents_NWB('nwb_file', nwb_file); % Checks the specified NWB file. Gives a list with the available events in the NWB file for the user to choose which one to load -% events = bz_LoadEvents_NWB('nwb_file', nwb_file, 'PulseStim_5V_77777ms_LD12'); % Loads the specified event type from the NWB file -% events = bz_LoadEvents_NWB('PulseStim_5V_77777ms_LD12'); % Checks the current directory for a NWB file. Loads the specified file - - -% The events are loaded straight from the NWB file. No extra files needed. - -% This code is based on the bz_LoadEvents -% Konstantinos Nasiotis 2019 - - -%% - p = inputParser; - addParameter(p,'nwb_file','',@isstr); - - if size(varargin,2)>2 - eventsName = varargin{3}; - varargin = {varargin{1,1:2}}; - parse(p,varargin{:}) - nwb_file = p.Results.nwb_file; - elseif size(varargin,2)==1 - nwb_file = []; - eventsName = varargin{1}; - elseif size(varargin,2)==2 - parse(p,varargin{:}) - nwb_file = p.Results.nwb_file; - elseif size(varargin,2)==0 - nwb_file = []; - end - - %% let's check that there is an appropriate NWB file - if isempty(nwb_file) - %disp('No nwb_file given, so we look for a *nwb file in the current directory') - d = dir('*nwb'); - if length(d) > 1 % we assume one .nwb file or this should break - error('there is more than one .nwb file in this directory?'); - elseif length(d) == 0 - d = dir('*nwb'); - if isempty(d) - error('could not find an nwb file..') - end - end - nwb_file = fullfile(d.folder, d.name); - end - - nwb2 = nwbRead(nwb_file); - - % Check if an events field exists in the dataset - try - events_exist_here = ~isempty(nwb2.stimulus_presentation); - if ~events_exist_here - disp('No events in this .nwb file') - return - - else - all_event_keys = keys(nwb2.stimulus_presentation); - disp(' ') - disp('The following event types are present in this dataset') - disp('------------------------------------------------') - for iEvent = 1:length(all_event_keys) - disp(all_event_keys{iEvent}) - end - disp(' ') - end - catch - disp('No events in this .nwb file') - return - end - - - % Check if a specific event was called to be loaded. If not, display a - % pop-up list with the available events for selection - if ~exist('eventsName','var') - [iEvent, ~] = listdlg('PromptString','Which event type would you like to load?',... - 'ListString',all_event_keys,'SelectionMode','single'); - else - if sum(ismember(all_event_keys, eventsName))>0 - iEvent = find(ismember(all_event_keys, eventsName)); - end - end - - events = struct; - - events.timestamps = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).timestamps.load;% neuroscope compatible matrix with 1-2 columns - [starts stops] (in seconds) - events.detectorinfo.detectorname = 'N/A';% substructure with information about the detection method (fields below) - events.detectorinfo.detectionparms = 'N/A';% parameters used for detection - events.detectorinfo.detectiondate = 'N/A';% date of detection - events.detectorinfo.detectionintervals = 'N/A';% [start stop] pairs of intervals used for detection (optional) -% events.detectorinfo.detectionchannel = ;% channel used for detection (optional) - -% % Optional -% events.amplitudes = 1;% [Nx1 matrix] -% events.frequencies = 1;% [Nx1 matrix] -% events.durations = 1;% [Nx1 matrix] - - - % I add here the rest of the parameters that are saved on the nwb file - events.starting_time_unit = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).starting_time_unit; - events.timestamps_interval = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).timestamps_interval; - events.timestamps_unit = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).timestamps_unit; - events.comments = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).comments; - events.control = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).control; - events.control_description = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).control_description; - events.data = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).data; - events.data_conversion = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).data_conversion; - events.data_resolution = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).data_resolution; - events.data_unit = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).data_unit; - events.description = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).description; - events.starting_time = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).starting_time; - events.starting_time_rate = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).starting_time_rate; - events.help = nwb2.stimulus_presentation.get(all_event_keys{iEvent}).help; - - -end - diff --git a/io/NWB/example_script_to_add_fields_from_Buzcode_format.m b/io/NWB/example_script_to_add_fields_from_Buzcode_format.m new file mode 100644 index 00000000..85e6b6d5 --- /dev/null +++ b/io/NWB/example_script_to_add_fields_from_Buzcode_format.m @@ -0,0 +1,42 @@ + + +% test standard formats +folder_path = 'C:\Users\McGill\Documents\GitHub\buzcode\tutorials\exampleDataStructs\20170505_396um_0um_merge'; + + +%% Start Adding fields to NWB + +% Get info from the rhd file +xml = GetXMLInfo(folder_path); + +% Add general info to the NWB file +nwb = GeneralInfo(xml); + +% Add electrode information +nwb = addElectrodeInfo(xml, nwb); + +% Add Behavior information +nwb = addBehavior(xml, nwb); + +% Add Behavior information +nwb = addTrials(xml, nwb); + +% Add events +nwb = addEvents(xml, nwb); + +% Add electrophysiological data (.lfp, .eeg) +nwb = addElectrophysiology(xml, nwb); + +% Add units info +nwb = addUnitsInfo(xml, nwb); + +%% Export to nwb +nwbExport(nwb, 'F:\NWBtoBuzcode\test_Buzcode_Standards.nwb') + + + + + + + + diff --git a/io/NWB/get_SessionInfoFromNWB.m b/io/NWB/get_SessionInfoFromNWB.m index b5cd4f19..75eacd9d 100644 --- a/io/NWB/get_SessionInfoFromNWB.m +++ b/io/NWB/get_SessionInfoFromNWB.m @@ -25,15 +25,19 @@ try all_raw_keys = keys(nwb2.acquisition); - for iKey = 1:length(all_raw_keys) - if ismember(all_raw_keys{iKey}, {'ECoG','bla bla bla'}) %%%%%%%% ADD MORE HERE, DON'T KNOW WHAT THE STANDARD FORMATS ARE - iRawDataKey = iKey; - RawDataPresent = 1; - else - RawDataPresent = 0; + if ~isempty(all_raw_keys) + for iKey = 1:length(all_raw_keys) + if ismember(all_raw_keys{iKey}, {'ECoG','bla bla bla'}) %%%%%%%% ADD MORE HERE, DON'T KNOW WHAT THE STANDARD FORMATS ARE + iRawDataKey = iKey; + RawDataPresent = 1; + else + RawDataPresent = 0; + end end + % nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries.get('all_lfp').data + else + RawDataPresent = 0; end - % nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries.get('all_lfp').data catch RawDataPresent = 0; end @@ -45,14 +49,18 @@ % Check if the data is in LFP format all_lfp_keys = keys(nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries); - for iKey = 1:length(all_lfp_keys) - if ismember(all_lfp_keys{iKey}, {'lfp','bla bla bla'}) %%%%%%%% ADD MORE HERE, DON'T KNOW WHAT THE STANDARD FORMATS ARE - iLFPDataKey = iKey; - LFPDataPresent = 1; - break % Once you find the data don't look for other keys/trouble - else - LFPDataPresent = 0; + if ~isempty(all_lfp_keys) + for iKey = 1:length(all_lfp_keys) + if ismember(all_lfp_keys{iKey}, {'lfp','bla bla bla'}) %%%%%%%% ADD MORE HERE, DON'T KNOW WHAT THE STANDARD FORMATS ARE + iLFPDataKey = iKey; + LFPDataPresent = 1; + break % Once you find the data don't look for other keys/trouble + else + LFPDataPresent = 0; + end end + else + LFPDataPresent = 0; end catch LFPDataPresent = 0; @@ -71,7 +79,7 @@ if RawDataPresent sessionInfo.nChannels = nwb2.acquisition.get(all_raw_keys{iRawDataKey}).data.dims(2); sessionInfo.samples_NWB = nwb2.acquisition.get(all_raw_keys{iRawDataKey}).data.dims(1); - sessionInfo.rates.wideband = nwb2.acquisition.get(all_raw_keys{iRawDataKey}).starting_time_rate; + sessionInfo.rates.wideband = nwb2.acquisition.get(all_raw_keys{iRawDataKey}).starting_time_rate; % THIS SHOULD BE 20,000 HZ FOR THE TUTORIAL sessionInfo.rates.lfp = 1250; %1250 - DEFAULT - CHECK THIS sessionInfo.lfpSampleRate = 1250; %1250 - DEFAULT - CHECK THIS % tH lfpSampleRate bypasses @@ -82,9 +90,8 @@ sessionInfo.rates.lfp = nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries.get(all_lfp_keys{iLFPDataKey}).starting_time_rate; %I assign the LFP sampling rate that was already used. Not sure yet if a different value than 1250 Hz will cause problems sessionInfo.lfpSampleRate = nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries.get(all_lfp_keys{iLFPDataKey}).starting_time_rate; %I assign the LFP sampling rate that was already used. Not sure yet if a different value than 1250 Hz will cause problems - if sessionInfo.rates.wideband <1250 - warning 'Something weird will happen. Not thoroughly tested yet' - end + sessionInfo.LFPDataPresent = LFPDataPresent; % This is added for easy retrieval from bz_GetLFP + sessionInfo.LFPDataKey = all_lfp_keys{iLFPDataKey}; % This is added for easy retrieval from bz_GetLFP end @@ -97,7 +104,7 @@ % group = nwb2.general_extracellular_ephys_electrodes.vectordata.get('group').data; % group_name = nwb2.general_extracellular_ephys_electrodes.vectordata.get('group_name').data; % imp = nwb2.general_extracellular_ephys_electrodes.vectordata.get('imp').data.load; -location = nwb2.general_extracellular_ephys_electrodes.vectordata.get('location').data; +location = nwb2.general_extracellular_ephys_electrodes.vectordata.get('location').data.load; % shank = nwb2.general_extracellular_ephys_electrodes.vectordata.get('shank').data.load; % x = nwb2.general_extracellular_ephys_electrodes.vectordata.get('x').data.load; % y = nwb2.general_extracellular_ephys_electrodes.vectordata.get('y').data.load; @@ -105,12 +112,10 @@ - - % nGroups = sum(unique(shank)>0); % -1 doesn't belong to a Shank Group -nGroups = length(unique(nwb2.general_extracellular_ephys_electrodes.vectordata.get('group_name').data)); -uniqueGroupNames = unique(nwb2.general_extracellular_ephys_electrodes.vectordata.get('group_name').data); -groupNames = nwb2.general_extracellular_ephys_electrodes.vectordata.get('group_name').data; +nGroups = length(unique(nwb2.general_extracellular_ephys_electrodes.vectordata.get('group_name').data.load)); +uniqueGroupNames = unique(nwb2.general_extracellular_ephys_electrodes.vectordata.get('group_name').data.load); +groupNames = nwb2.general_extracellular_ephys_electrodes.vectordata.get('group_name').data.load; sessionInfo.spikeGroups.nGroups = nGroups; @@ -140,15 +145,17 @@ % Add Region Info -for iChannel = 1:sessionInfo.nChannels +% If the recording has additional channels (behavioral channels - they wouldn't have a location field - +% I assume here that the numbering of the channels starts with the electrophysiological, and the behavioral are concatenated after) +% Thats why I use length(groupNames) here and not sessionInfo.nChannels + +for iChannel = 1:length(groupNames) sessionInfo.region{iChannel} = location{iChannel}; end - - % Get the rest of the Info -sessionInfo.nBits = 16; % ASSUMING THAT NWB GIVES INT16 PRECISION +sessionInfo.nBits = 16; % ASSUMING THAT NWB HAS SAVED DATA IN INT16 PRECISION sessionInfo.rates.video = 0; sessionInfo.FileName = nwb2.identifier;%%%%%%%%%%%%%%%% no extension - I DON'T USE THE FILENAME OF THE NWB HERE JUST IN CASE SOMEONE CHANGED IT. % sessionInfo.SampleTime = 50; % 50 no idea diff --git a/io/bz_GetLFP.m b/io/bz_GetLFP.m index 0f17ee72..7720c4b5 100755 --- a/io/bz_GetLFP.m +++ b/io/bz_GetLFP.m @@ -53,6 +53,8 @@ % Copyright (C) 2004-2011 by Michaël Zugaro % editied by David Tingley, 2017 % +% added support for NWB by Konstantinos Nasiotis, 2019 + % NOTES % -'select' option has been removed, it allowed switching between 0 and 1 % indexing. This should no longer be necessary with .lfp.mat structs @@ -102,7 +104,20 @@ elseif length(d) == 0 d = dir([basepath filesep '*eeg']); if isempty(d) - error('could not find an lfp/eeg file..') + % If no .lfp file is present, check for a .nwb + d = dir([basepath filesep '*.nwb']); + if length(d) > 1 % If more than one .nwb files exist in the directory, select which one to load from + warning('there is more than one .nwb file in this directory'); + % Check if a specific behavior was called to be loaded. If not, display a + % pop-up list with the available behaviors for selection + [iNWBFile, ~] = listdlg('PromptString','Which NWB file would you like to load?',... + 'ListString',{d.name},'SelectionMode','single'); + d = d(iNWBFile); + else + if isempty(d) + error('could not find an lfp/eeg/nwb file..') + end + end end end lfp.Filename = d.name; @@ -124,15 +139,47 @@ elseif length(d) == 0 d = dir([basepath filesep basename '.eeg']); if isempty(d) - error('could not find an lfp/eeg file..') + d = dir([basepath filesep basename '.nwb']); + if length(d) > 1 % we assume one .nwb file or this should break + warning('there is more than one .nwb file in this directory'); + % Check if a specific behavior was called to be loaded. If not, display a + % pop-up list with the available behaviors for selection + [iNWBFile, ~] = listdlg('PromptString','Which NWB file would you like to load?',... + 'ListString',{d.name},'SelectionMode','single'); + d = d(iNWBFile); + else + if isempty(d) + error('could not find an lfp/eeg/nwb file..') + end + end end end lfp.Filename = d.name; end + +%% Add a flag that shows if variables are loaded from a .nwb file +% nwb has everything saved in a single file + +if ~isempty(strfind(lfp.Filename,'nwb')) + nwb_loaded = 1; +else + nwb_loaded = 0; +end + %% things we can parse from sessionInfo or xml file -sessionInfo = bz_getSessionInfo(basepath, 'noPrompts', noPrompts); +if ~nwb_loaded + sessionInfo = bz_getSessionInfo(basepath, 'noPrompts', noPrompts); +else + % Load the .nwb info + nwb2 = nwbRead(lfp.Filename); + sessionInfo = get_SessionInfoFromNWB(nwb2); % The get_SessionInfoFromNWB loads only the needed fields for the functions tested. Maybe improve + + if ~sessionInfo.LFPDataPresent + error('No LFP data detected in this file. Make sure that there is a field within: nwb2.processing.get("ecephys").nwbdatainterface.get("LFP").electricalseries') + end +end try samplingRate = sessionInfo.lfpSampleRate; @@ -165,11 +212,38 @@ % Load data and put into struct % we assume 0-indexing like neuroscope, but bz_LoadBinary uses 1-indexing to % load.... - lfp(i).data = bz_LoadBinary([basepath filesep lfp.Filename],... - 'duration',double(lfp(i).duration),... - 'frequency',samplingRate,'nchannels',sessionInfo.nChannels,... - 'start',double(lfp(i).interval(1)),'channels',channels+1,... - 'downsample',downsamplefactor); + + if ~nwb_loaded + lfp(i).data = bz_LoadBinary([basepath filesep lfp.Filename],... + 'duration',double(lfp(i).duration),... + 'frequency',samplingRate,'nchannels',sessionInfo.nChannels,... + 'start',double(lfp(i).interval(1)),'channels',channels+1,... + 'downsample',downsamplefactor); + + else + %% NWB + % Memory preallocation is used for faster assignment + if intervals(i,2) == Inf + use_this_bracket = sessionInfo.samples_NWB/sessionInfo.lfpSampleRate; + data_temp = zeros(length(channels),ceil((use_this_bracket - intervals(i,1))*samplingRateLFP)); + else + use_this_bracket = intervals(i,2); + data_temp = zeros(length(channels),floor((intervals(i,2) - intervals(i,1))*samplingRateLFP)); + end + + % Differentiate between sequential loading of each channel and + % simultaneous loading of neighboring channels + if all(diff(channels) == 1) + data_temp = nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries.get(sessionInfo.LFPDataKey).data.load([channels(1)+1, intervals(i,1)*samplingRate+1],[1, downsamplefactor],[channels(end)+1, use_this_bracket*samplingRate]); + else + for iChannel = 1:length(channels) + disp(['Now concatenating channel: ' num2str(channels(iChannel))]) + data_temp(iChannel,:) = nwb2.processing.get('ecephys').nwbdatainterface.get('LFP').electricalseries.get(sessionInfo.LFPDataKey).data.load([channels(iChannel)+1, intervals(i,1)*samplingRate+1],[1, downsamplefactor],[channels(iChannel)+1, use_this_bracket*samplingRate]); + end + end + lfp(i).data = data_temp'; clear data_temp + end + lfp(i).timestamps = [lfp(i).interval(1):(1/samplingRateLFP):... (lfp(i).interval(1)+(length(lfp(i).data)-1)/... samplingRateLFP)]'; diff --git a/io/bz_GetSpikes.m b/io/bz_GetSpikes.m index 0d62b758..5c14bd44 100755 --- a/io/bz_GetSpikes.m +++ b/io/bz_GetSpikes.m @@ -60,6 +60,9 @@ % % % written by David Tingley, 2017 + +% Added NWB support by Konstantinos Nasiotis 2019 + %% Deal With Inputs spikeGroupsValidation = @(x) assert(isnumeric(x) || strcmp(x,'all'),... 'spikeGroups must be numeric or "all"'); @@ -87,263 +90,359 @@ noPrompts = p.Results.noPrompts; onlyLoad = p.Results.onlyLoad; +% First check if there is an .xml or a sessionInfo in the current path +try + [sessionInfo] = bz_getSessionInfo(basepath, 'noPrompts', noPrompts); + nwb_loaded = 0; + +catch + % If the sessionInfo can't be created, check for a .nwb + d = dir([basepath filesep '*.nwb']); + if length(d) > 1 % If more than one .nwb files exist in the directory, select which one to load from + warning('there is more than one .nwb file in this directory'); + [iNWBFile, ~] = listdlg('PromptString','Which NWB file would you like to load?',... + 'ListString',{d.name},'SelectionMode','single'); + d = d(iNWBFile); + + else + if isempty(d) + error('could not find an .xml or a .nwb in the current path') + end + end + nwb_file = d.name; + [the_path, ~, ~] = fileparts(nwb_file); + + % Load NWB info + nwb2 = nwbRead(nwb_file); + sessionInfo = get_SessionInfoFromNWB(nwb2); % The get_SessionInfoFromNWB loads only the needed fields for the functions tested. Maybe improve -[sessionInfo] = bz_getSessionInfo(basepath, 'noPrompts', noPrompts); + nwb_loaded = 1; +end spikes.samplingRate = sessionInfo.rates.wideband; nChannels = sessionInfo.nChannels; -%% if the cellinfo file exist and we don't want to re-load files -if exist([basepath filesep sessionInfo.FileName '.spikes.cellinfo.mat'],'file') && forceReload == false - disp('loading spikes from cellinfo file..') - load([basepath filesep sessionInfo.FileName '.spikes.cellinfo.mat']) - %Check that the spikes structure fits cellinfo requirements - [iscellinfo] = bz_isCellInfo(spikes); - switch iscellinfo - case false - warning(['The spikes structure in baseName.spikes.cellinfo.mat ',... - 'does not fit buzcode standards. Sad.']) - end - - %If regions have been added since creation... add them - if ~isfield(spikes,'region') & isfield(sessionInfo,'region') - for cc = 1:spikes.numcells - spikes.region{cc} = sessionInfo.region{spikes.maxWaveformCh(cc)==sessionInfo.channels}; +%% Get spikes +if ~nwb_loaded + + %% if the cellinfo file exist and we don't want to re-load files + if exist([basepath filesep sessionInfo.FileName '.spikes.cellinfo.mat'],'file') && forceReload == false + disp('loading spikes from cellinfo file..') + load([basepath filesep sessionInfo.FileName '.spikes.cellinfo.mat']) + %Check that the spikes structure fits cellinfo requirements + [iscellinfo] = bz_isCellInfo(spikes); + switch iscellinfo + case false + warning(['The spikes structure in baseName.spikes.cellinfo.mat ',... + 'does not fit buzcode standards. Sad.']) end - - if saveMat - save([basepath filesep sessionInfo.FileName '.spikes.cellinfo.mat'],'spikes') + + %If regions have been added since creation... add them + if ~isfield(spikes,'region') & isfield(sessionInfo,'region') + for cc = 1:spikes.numcells + spikes.region{cc} = sessionInfo.region{spikes.maxWaveformCh(cc)==sessionInfo.channels}; + end + + if saveMat + save([basepath filesep sessionInfo.FileName '.spikes.cellinfo.mat'],'spikes') + end end - end - -else % do the below then filter by inputs... (Load from clu/res/fet) - - if ~noPrompts & saveMat == 0 %Inform the user that they should save a file for later - savebutton = questdlg(['Would you like to save your spikes in ',... - sessionInfo.FileName,'.spikes.cellinfo.mat? ',... - 'This will save significant load time later.']); - if strcmp(savebutton,'Yes'); saveMat = true; end - end - -disp('loading spikes from clu/res/spk files..') -% find res/clu/fet/spk files here -cluFiles = dir([basepath filesep '*.clu*']); -resFiles = dir([basepath filesep '*.res.*']); -if any(getWaveforms) - spkFiles = dir([basepath filesep '*.spk*']); -end + else % do the below then filter by inputs... (Load from clu/res/fet) + + if ~noPrompts & saveMat == 0 %Inform the user that they should save a file for later + savebutton = questdlg(['Would you like to save your spikes in ',... + sessionInfo.FileName,'.spikes.cellinfo.mat? ',... + 'This will save significant load time later.']); + if strcmp(savebutton,'Yes'); saveMat = true; end + end -% remove *temp*, *autosave*, and *.clu.str files/directories -tempFiles = zeros(length(cluFiles),1); -for i = 1:length(cluFiles) - dummy = strsplit(cluFiles(i).name, '.'); % Check whether the component after the last dot is a number or not. If not, exclude the file/dir. - if ~isempty(findstr('temp',cluFiles(i).name)) | ~isempty(findstr('autosave',cluFiles(i).name)) | isempty(str2num(dummy{length(dummy)})) | find(contains(dummy, 'clu')) ~= length(dummy)-1 - tempFiles(i) = 1; + disp('loading spikes from clu/res/spk files..') + % find res/clu/fet/spk files here + cluFiles = dir([basepath filesep '*.clu*']); + resFiles = dir([basepath filesep '*.res.*']); + if any(getWaveforms) + spkFiles = dir([basepath filesep '*.spk*']); end -end -cluFiles(tempFiles==1)=[]; -tempFiles = zeros(length(resFiles),1); -for i = 1:length(resFiles) - if ~isempty(findstr('temp',resFiles(i).name)) | ~isempty(findstr('autosave',resFiles(i).name)) - tempFiles(i) = 1; + + + % remove *temp*, *autosave*, and *.clu.str files/directories + tempFiles = zeros(length(cluFiles),1); + for i = 1:length(cluFiles) + dummy = strsplit(cluFiles(i).name, '.'); % Check whether the component after the last dot is a number or not. If not, exclude the file/dir. + if ~isempty(findstr('temp',cluFiles(i).name)) | ~isempty(findstr('autosave',cluFiles(i).name)) | isempty(str2num(dummy{length(dummy)})) | find(contains(dummy, 'clu')) ~= length(dummy)-1 + tempFiles(i) = 1; + end end -end -if any(getWaveforms) - resFiles(tempFiles==1)=[]; - tempFiles = zeros(length(spkFiles),1); - for i = 1:length(spkFiles) - if ~isempty(findstr('temp',spkFiles(i).name)) | ~isempty(findstr('autosave',spkFiles(i).name)) + cluFiles(tempFiles==1)=[]; + tempFiles = zeros(length(resFiles),1); + for i = 1:length(resFiles) + if ~isempty(findstr('temp',resFiles(i).name)) | ~isempty(findstr('autosave',resFiles(i).name)) tempFiles(i) = 1; end end - spkFiles(tempFiles==1)=[]; -end + if any(getWaveforms) + resFiles(tempFiles==1)=[]; + tempFiles = zeros(length(spkFiles),1); + for i = 1:length(spkFiles) + if ~isempty(findstr('temp',spkFiles(i).name)) | ~isempty(findstr('autosave',spkFiles(i).name)) + tempFiles(i) = 1; + end + end + spkFiles(tempFiles==1)=[]; + end -if isempty(cluFiles) - disp('no clu files found...') - spikes = []; - return -end + if isempty(cluFiles) + disp('no clu files found...') + spikes = []; + return + end -% ensures we load in sequential order (forces compatibility with FMAT -% ordering) -for i = 1:length(cluFiles) - temp = strsplit(cluFiles(i).name,'.'); - shanks(i) = str2num(temp{length(temp)}); -end -[shanks ind] = sort(shanks); -cluFiles = cluFiles(ind); %Bug here if there are any files x.clu.x that are not your desired clus -resFiles = resFiles(ind); -if any(getWaveforms) - spkFiles = spkFiles(ind); -end + % ensures we load in sequential order (forces compatibility with FMAT + % ordering) + for i = 1:length(cluFiles) + temp = strsplit(cluFiles(i).name,'.'); + shanks(i) = str2num(temp{length(temp)}); + end + [shanks ind] = sort(shanks); + cluFiles = cluFiles(ind); %Bug here if there are any files x.clu.x that are not your desired clus + resFiles = resFiles(ind); + if any(getWaveforms) + spkFiles = spkFiles(ind); + end -% check if there are matching #'s of files -if length(cluFiles) ~= length(resFiles) & length(cluFiles) ~= length(spkFiles) - error('found an incorrect number of res/clu/spk files...') -end + % check if there are matching #'s of files + if length(cluFiles) ~= length(resFiles) & length(cluFiles) ~= length(spkFiles) + error('found an incorrect number of res/clu/spk files...') + end -% use the .clu files to get spike ID's and generate UID and spikeGroup -% use the .res files to get spike times -count = 1; + % use the .clu files to get spike ID's and generate UID and spikeGroup + % use the .res files to get spike times + count = 1; -if isempty(sessionInfo.spikeGroups.groups) - sessionInfo.spikeGroups = sessionInfo.AnatGrps; -end -for i=1:length(cluFiles) - disp(['working on ' cluFiles(i).name]) - - temp = strsplit(cluFiles(i).name,'.'); - shankID = str2num(temp{length(temp)}); %shankID is the spikegroup number - clu = load(fullfile(basepath,cluFiles(i).name)); - clu = clu(2:end); % toss the first sample to match res/spk files - res = load(fullfile(basepath,resFiles(i).name)); - spkGrpChans = sessionInfo.spikeGroups.groups{shankID}; % we'll eventually want to replace these two lines - - if any(getWaveforms) && sum(clu)>0 %bug fix if no clusters - nSamples = sessionInfo.spikeGroups.nSamples(shankID); - % load waveforms - chansPerSpikeGrp = length(sessionInfo.spikeGroups.groups{shankID}); - fid = fopen(fullfile(basepath,spkFiles(i).name),'r'); - wav = fread(fid,[1 inf],'int16=>int16'); - try %bug in some spk files... wrong number of samples? - wav = reshape(wav,chansPerSpikeGrp,nSamples,[]); - catch - if strcmp(getWaveforms,'force') - wav = nan(chansPerSpikeGrp,nSamples,length(clu)); - display([spkFiles(i).name,' error.']) - else - error(['something is wrong with ',spkFiles(i).name,... - ' Use ''getWaveforms'', false to skip waveforms or ',... - '''getWaveforms'', ''force'' to write nans on bad shanks.']) - end - end - wav = permute(wav,[3 1 2]); + if isempty(sessionInfo.spikeGroups.groups) + sessionInfo.spikeGroups = sessionInfo.AnatGrps; end - - cells = unique(clu); - % remove MUA and NOISE clusters... - cells(cells==0) = []; - cells(cells==1) = []; % consider adding MUA as another input argument...? - - for c = 1:length(cells) - spikes.UID(count) = count; % this only works if all shanks are loaded... how do we optimize this? - ind = find(clu == cells(c)); - spikes.times{count} = res(ind) ./ spikes.samplingRate; - spikes.shankID(count) = shankID; - spikes.cluID(count) = cells(c); - - %Waveforms - if any(getWaveforms) - wvforms = squeeze(mean(wav(ind,:,:)))-mean(mean(mean(wav(ind,:,:)))); % mean subtract to account for slower (theta) trends - if prod(size(wvforms))==length(wvforms)%in single-channel groups wvforms will squeeze too much and will have amplitude on D1 rather than D2 - wvforms = wvforms';%fix here - end - for t = 1:size(wvforms,1) - [a(t) b(t)] = max(abs(wvforms(t,:))); - end - [aa bb] = max(a,[],2); - spikes.rawWaveform{count} = wvforms(bb,:); - spikes.maxWaveformCh(count) = spkGrpChans(bb); - %Regions (needs waveform peak) - if isfield(sessionInfo,'region') %if there is regions field in your metadata - spikes.region{count} = sessionInfo.region{find(spkGrpChans(bb)==sessionInfo.channels)}; - elseif isfield(sessionInfo,'Units') %if no regions, but unit region from xml via Loadparamteres - %Find the xml Unit that matches group/cluster - unitnum = cellfun(@(X,Y) X==spikes.shankID(count) && Y==spikes.cluID(count),... - {sessionInfo.Units(:).spikegroup},{sessionInfo.Units(:).cluster}); - if sum(unitnum) == 0 - display(['xml Missing Unit - spikegroup: ',... - num2str(spikes.shankID(count)),' cluster: ',... - num2str(spikes.cluID(count))]) - spikes.region{count} = 'missingxml'; - else %possible future bug: two xml units with same group/clu... - spikes.region{count} = sessionInfo.Units(unitnum).structure; + for i=1:length(cluFiles) + disp(['working on ' cluFiles(i).name]) + + temp = strsplit(cluFiles(i).name,'.'); + shankID = str2num(temp{length(temp)}); %shankID is the spikegroup number + clu = load(fullfile(basepath,cluFiles(i).name)); + clu = clu(2:end); % toss the first sample to match res/spk files + res = load(fullfile(basepath,resFiles(i).name)); + spkGrpChans = sessionInfo.spikeGroups.groups{shankID}; % we'll eventually want to replace these two lines + + if any(getWaveforms) && sum(clu)>0 %bug fix if no clusters + nSamples = sessionInfo.spikeGroups.nSamples(shankID); + % load waveforms + chansPerSpikeGrp = length(sessionInfo.spikeGroups.groups{shankID}); + fid = fopen(fullfile(basepath,spkFiles(i).name),'r'); + wav = fread(fid,[1 inf],'int16=>int16'); + try %bug in some spk files... wrong number of samples? + wav = reshape(wav,chansPerSpikeGrp,nSamples,[]); + catch + if strcmp(getWaveforms,'force') + wav = nan(chansPerSpikeGrp,nSamples,length(clu)); + display([spkFiles(i).name,' error.']) + else + error(['something is wrong with ',spkFiles(i).name,... + ' Use ''getWaveforms'', false to skip waveforms or ',... + '''getWaveforms'', ''force'' to write nans on bad shanks.']) end + end + wav = permute(wav,[3 1 2]); + end + + cells = unique(clu); + % remove MUA and NOISE clusters... + cells(cells==0) = []; + cells(cells==1) = []; % consider adding MUA as another input argument...? + + for c = 1:length(cells) + spikes.UID(count) = count; % this only works if all shanks are loaded... how do we optimize this? + ind = find(clu == cells(c)); + spikes.times{count} = res(ind) ./ spikes.samplingRate; + spikes.shankID(count) = shankID; + spikes.cluID(count) = cells(c); + + %Waveforms + if any(getWaveforms) + wvforms = squeeze(mean(wav(ind,:,:)))-mean(mean(mean(wav(ind,:,:)))); % mean subtract to account for slower (theta) trends + if prod(size(wvforms))==length(wvforms)%in single-channel groups wvforms will squeeze too much and will have amplitude on D1 rather than D2 + wvforms = wvforms';%fix here + end + for t = 1:size(wvforms,1) + [a(t) b(t)] = max(abs(wvforms(t,:))); + end + [aa bb] = max(a,[],2); + spikes.rawWaveform{count} = wvforms(bb,:); + spikes.maxWaveformCh(count) = spkGrpChans(bb); + %Regions (needs waveform peak) + if isfield(sessionInfo,'region') %if there is regions field in your metadata + spikes.region{count} = sessionInfo.region{find(spkGrpChans(bb)==sessionInfo.channels)}; + elseif isfield(sessionInfo,'Units') %if no regions, but unit region from xml via Loadparamteres + %Find the xml Unit that matches group/cluster + unitnum = cellfun(@(X,Y) X==spikes.shankID(count) && Y==spikes.cluID(count),... + {sessionInfo.Units(:).spikegroup},{sessionInfo.Units(:).cluster}); + if sum(unitnum) == 0 + display(['xml Missing Unit - spikegroup: ',... + num2str(spikes.shankID(count)),' cluster: ',... + num2str(spikes.cluID(count))]) + spikes.region{count} = 'missingxml'; + else %possible future bug: two xml units with same group/clu... + spikes.region{count} = sessionInfo.Units(unitnum).structure; + end + end + clear a aa b bb end - clear a aa b bb - end - - count = count + 1; + + count = count + 1; + end end -end -if ~isempty(onlyLoad) - toRemove = true(size(spikes.UID)); - for cc = 1:size(onlyLoad,1) - whichUID = ismember(spikes.shankID,onlyLoad(cc,1)) & ismember(spikes.cluID,onlyLoad(cc,2)); - toRemove(whichUID) = false; - if ~any(whichUID) - display(['No unit with shankID:',num2str(onlyLoad(cc,1)),... - ' cluID:',num2str(onlyLoad(cc,2))]) + if ~isempty(onlyLoad) + toRemove = true(size(spikes.UID)); + for cc = 1:size(onlyLoad,1) + whichUID = ismember(spikes.shankID,onlyLoad(cc,1)) & ismember(spikes.cluID,onlyLoad(cc,2)); + toRemove(whichUID) = false; + if ~any(whichUID) + display(['No unit with shankID:',num2str(onlyLoad(cc,1)),... + ' cluID:',num2str(onlyLoad(cc,2))]) + end end + spikes = removeCells(toRemove,spikes,getWaveforms); end - spikes = removeCells(toRemove,spikes,getWaveforms); -end -spikes.sessionName = sessionInfo.FileName; -end + spikes.sessionName = sessionInfo.FileName; + end -%% save to buzcode format (before exclusions) -if saveMat - save([basepath filesep sessionInfo.FileName '.spikes.cellinfo.mat'],'spikes') -end + %% save to buzcode format (before exclusions) + if saveMat + save([basepath filesep sessionInfo.FileName '.spikes.cellinfo.mat'],'spikes') + end + + + %% EXCLUSIONS %% + %filter by spikeGroups input + if ~strcmp(spikeGroups,'all') + [toRemove] = ~ismember(spikes.shankID,spikeGroups); + spikes = removeCells(toRemove,spikes,getWaveforms); + end -%% EXCLUSIONS %% + %filter by region input + if ~isempty(region) + if ~isfield(spikes,'region') %if no region information in metadata + error(['You selected to load cells from region "',region,... + '", but there is no region information in your sessionInfo']) + end -%filter by spikeGroups input -if ~strcmp(spikeGroups,'all') - [toRemove] = ~ismember(spikes.shankID,spikeGroups); - spikes = removeCells(toRemove,spikes,getWaveforms); -end + toRemove = ~ismember(spikes.region,region); + if sum(toRemove)==length(spikes.UID) %if no cells from selected region + warning(['You selected to load cells from region "',region,... + '", but none of your cells are from that region']) + end -%filter by region input -if ~isempty(region) - if ~isfield(spikes,'region') %if no region information in metadata - error(['You selected to load cells from region "',region,... - '", but there is no region information in your sessionInfo']) - end - - toRemove = ~ismember(spikes.region,region); - if sum(toRemove)==length(spikes.UID) %if no cells from selected region - warning(['You selected to load cells from region "',region,... - '", but none of your cells are from that region']) + spikes = removeCells(toRemove,spikes,getWaveforms); end - - spikes = removeCells(toRemove,spikes,getWaveforms); -end - -%filter by UID input -if ~isempty(UID) - [toRemove] = ~ismember(spikes.UID,UID); - spikes = removeCells(toRemove,spikes,getWaveforms); -end -%% Generate spindices matrics -spikes.numcells = length(spikes.UID); -for cc = 1:spikes.numcells - groups{cc}=spikes.UID(cc).*ones(size(spikes.times{cc})); -end -if spikes.numcells>0 - alltimes = cat(1,spikes.times{:}); groups = cat(1,groups{:}); %from cell to array - [alltimes,sortidx] = sort(alltimes); groups = groups(sortidx); %sort both - spikes.spindices = [alltimes groups]; -end + %filter by UID input + if ~isempty(UID) + [toRemove] = ~ismember(spikes.UID,UID); + spikes = removeCells(toRemove,spikes,getWaveforms); + end -%% Check if any cells made it through selection -if isempty(spikes.times) | spikes.numcells == 0 - spikes = []; -end + %% Generate spindices matrics + spikes.numcells = length(spikes.UID); + for cc = 1:spikes.numcells + groups{cc}=spikes.UID(cc).*ones(size(spikes.times{cc})); + end + if spikes.numcells>0 + alltimes = cat(1,spikes.times{:}); groups = cat(1,groups{:}); %from cell to array + [alltimes,sortidx] = sort(alltimes); groups = groups(sortidx); %sort both + spikes.spindices = [alltimes groups]; + end + %% Check if any cells made it through selection + if isempty(spikes.times) | spikes.numcells == 0 + spikes = []; + end + + +else % Load from NWB file + + nNeurons = length(nwb2.units.id.data.load); + % Assign neuron to region based on the region that its Shank belongs to + % Get a single electrode that belongs in that shank + electrodes2Shank = nwb2.general_extracellular_ephys_electrodes.vectordata.get('group_name').data.load; + electrodes2Region = nwb2.general_extracellular_ephys_electrodes.vectordata.get('location').data.load; + + neurons2ShankNames = cell(1,nNeurons); + all_region = cell(1,nNeurons); + + for iNeuron = 1:length(nwb2.units.electrode_group.data) + path = nwb2.units.electrode_group.data(iNeuron).path; + inds = strfind(path, '/'); + neurons2ShankNames{iNeuron} = path(inds(end)+1:end); + + ElectrodesInThatShank = find(strcmp(strrep(electrodes2Shank,' ',''),neurons2ShankNames{iNeuron})); + all_region{iNeuron} = electrodes2Region{ElectrodesInThatShank(1)}; + + end + + uniqueShanks = unique(neurons2ShankNames,'legacy'); + shankID_of_selected_Neurons = zeros(1, nNeurons); + for iNeuron = 1:nNeurons + shankID_of_selected_Neurons(iNeuron) = find(strcmp(uniqueShanks, neurons2ShankNames{iNeuron})); + end + + %% Make the check here of what will be loaded + selected_neurons_UID = false(1,nNeurons); + selected_neurons_spikeGroups = false(1,nNeurons); + selected_neurons_region = false(1,nNeurons); + + % Check which neurons were selected + if ~isempty(UID) + selected_neurons_UID(UID) = true; + end + %Check which Shanks were selected + if ~isempty(spikeGroups) + for iShank = spikeGroups + selected_neurons_spikeGroups(find(shankID_of_selected_Neurons==iShank)) = true; + end + end + %Check which regions were selected + if ~isempty(region) + selected_neurons_region(find(contains(all_region,region))) = true; + end + + UID = find(selected_neurons_UID | selected_neurons_spikeGroups | selected_neurons_region); + + if isempty(UID) + warning(['Warning: The selection made didnt specify any subgroup of neurons. Loading all neurons!']) + UID = 1:nNeurons; + end + + + %% The first time that this function is loaded, make sure to save a .spikes.cellInfo.mat that contains information from all neurons + + if saveMat + UID_forSaveMAt = 1:nNeurons; + temp = get_the_spikes_from_selected_UIDs(UID_forSaveMAt, nwb2, sessionInfo, saveMat, all_region, neurons2ShankNames, the_path); clear temp + end + + spikes = get_the_spikes_from_selected_UIDs(UID, nwb2, sessionInfo, 0, all_region, neurons2ShankNames, the_path); + +end + end @@ -379,4 +478,60 @@ +% NWB function +function spikes = get_the_spikes_from_selected_UIDs(UID, nwb2, sessionInfo, saveMat, all_region, neurons2ShankNames, the_path) + + %% Get Spikes + spikes = struct; + spikes.samplingRate = sessionInfo.rates.wideband; + spikes.UID = UID; + + times = cell(1,length(UID)); + rawWaveform = cell(1,length(UID)); + spindices = []; + + entry = 0; + for iNeuron = UID + entry = entry+1; + if iNeuron == 1 + times_temp = nwb2.units.spike_times.data.load(1:sum(nwb2.units.spike_times_index.data.load(iNeuron))); + else + times_temp = nwb2.units.spike_times.data.load(sum(nwb2.units.spike_times_index.data.load(iNeuron-1))+1:sum(nwb2.units.spike_times_index.data.load(iNeuron))); + end + times{entry} = times_temp(times_temp~=0)'; + + % FIX + if ~isempty(nwb2.units.waveform_mean) + rawWaveform{entry} = nwb2.units.waveform_mean.data.load([iNeuron, 1], [iNeuron,Inf]); % The template waveform should be filled by: nwb2.units.waveform_mean + else + rawWaveform{entry} = []; + end + spindices = [spindices ; times{entry} ones(length(times{entry}),1)*iNeuron]; + end + % Spindices have to be sorted according to when each spike occured + [~,sortedIndices] = sort(spindices(:,1)); + spindices = spindices(sortedIndices,:); + + uniqueShanks = unique(neurons2ShankNames,'legacy'); + shankID_of_selected_Neurons = zeros(1, length(UID)); + for iNeuron = 1:length(UID) + shankID_of_selected_Neurons(iNeuron) = find(strcmp(uniqueShanks, neurons2ShankNames{UID(iNeuron)})); + end + + spikes.times = times; + spikes.shankID = shankID_of_selected_Neurons; + spikes.rawWaveform = rawWaveform; + spikes.sessionName = sessionInfo.FileName; + spikes.numcells = length(UID); + spikes.spindices = spindices; % This holds the timing of each spike, sorted, and the neuron it belongs to. + + % FIX + spikes.cluID = nwb2.units.vectordata.get('cluID').data.load'; + spikes.maxWaveformCh = nwb2.units.vectordata.get('maxWaveformCh').data.load'; + spikes.region = all_region(UID); + + if saveMat + save(fullfile(the_path, [sessionInfo.FileName '.spikes.cellinfo.mat']),'spikes') + end +end diff --git a/io/bz_LoadBehavior.m b/io/bz_LoadBehavior.m index ed5d62d2..95040679 100755 --- a/io/bz_LoadBehavior.m +++ b/io/bz_LoadBehavior.m @@ -4,6 +4,8 @@ % % %DLevenstein 2017 + +% Added NWB support: Konstantinos Nasiotis 2019 %% if ~exist('basePath','var') @@ -13,13 +15,20 @@ if ~exist('behaviorName','var') allBehaviorFiles = dir(fullfile(basePath,[baseName,'.','*','.behavior.mat'])); - [s,v] = listdlg('PromptString','Which behavior.mat would you like to load?',... - 'ListString',{allBehaviorFiles.name},'SelectionMode','single'); - if isempty(s) - behavior = []; filename = []; - return + + if ~isempty(allBehaviorFiles) + [s,v] = listdlg('PromptString','Which behavior.mat would you like to load?',... + 'ListString',{allBehaviorFiles.name},'SelectionMode','single'); + if isempty(s) + behavior = []; filename = []; + return + end + filename = fullfile(basePath,allBehaviorFiles(s).name); + else + check_for_NWB = 1; + filename = []; end - filename = fullfile(basePath,allBehaviorFiles(s).name); + else filename = fullfile(basePath,[baseName,'.',behaviorName,'.behavior.mat']); end @@ -27,12 +36,25 @@ if exist(filename,'file') behaviorstruct = load(filename); + check_for_NWB = 0; else - warning([filename,' does not exist...']) - behavior = []; - return + check_for_NWB = 1; end + +if check_for_NWB + if ~exist('behaviorName','var') + behaviorName = []; + end + behaviorstruct = loadBehaviorNWB(basePath, behaviorName); + if isempty(behaviorstruct) + behavior = []; filename = []; + return + end +end + + + varsInFile = fieldnames(behaviorstruct); if numel(varsInFile)==1 @@ -50,5 +72,192 @@ end + + + +%% Function that loads NWB fields to a Behavior structure +function behaviorstruct = loadBehaviorNWB(basePath, behaviorName) + + % Check for NWB file since no .behavior.mat files were found + d = dir(fullfile(basePath, '*.nwb')); + if length(d) > 1 % If more than one .nwb files exist in the directory, select which one to load from + warning('there is more than one .nwb file in this directory'); + [iNWBFile, ~] = listdlg('PromptString','Which NWB file would you like to load?',... + 'ListString',{d.name},'SelectionMode','single'); + d = d(iNWBFile); + elseif length(d) == 0 + d = dir('*nwb'); + if isempty(d) + disp('No .behavior.mat or NWB files present in this directory') + behaviorstruct = []; + return + end + end + nwb_file = fullfile(d.folder, d.name); + + + + if exist(nwb_file,'file') == 2 + nwb2 = nwbRead(nwb_file); + + % Check if the behavior field exists in the dataset + all_processing_keys = keys(nwb2.processing); + behavior_FIELD_exists = sum(ismember(all_processing_keys, 'behavior'))==1; + + if behavior_FIELD_exists + allBehaviorKeys = keys(nwb2.processing.get('behavior').nwbdatainterface); + + % Remove the trials_ keys since they are just the trials sub-field + % within the events field of the behavior.mat files + ii = []; + for iKey = 1:length(allBehaviorKeys) + if isempty(strfind(allBehaviorKeys{iKey}, 'trials_')) + ii = [ii,iKey]; + end + end + allBehaviorKeys = {allBehaviorKeys{ii}}; + + + behavior_exist_here = ~isempty(allBehaviorKeys); + if ~behavior_exist_here + disp('No entry is filled in the Behavior field') + behaviorstruct = []; + return + else + disp(' ') + disp('The following behavior types are present in this dataset') + disp('------------------------------------------------') + for iBehavior = 1:length(allBehaviorKeys) + disp(allBehaviorKeys{iBehavior}) + end + disp(' ') + end + + + % Check if a specific behavior was called to be loaded. If not, display a + % pop-up list with the available behaviors for selection + + % NOTE: THE PREFIX TRIALS_ WILL BE IGNORED SINCE IT SIMPLY + % IMPLIES THE TRIALS THAT ARE WITHIN THE BEHAVIOR.MAT FILE + + if isempty(behaviorName) + [iBehavior, ~] = listdlg('PromptString','Which behavior type would you like to load?',... + 'ListString',allBehaviorKeys,'SelectionMode','single'); + else + if sum(ismember(allBehaviorKeys, behaviorName))==1 % If the behavior exists, find its index + iBehavior = find(ismember(allBehaviorKeys, behaviorName)); + else + disp('The requested behavior doesnt exist in this NWB file') + behaviorstruct = []; + return + end + end + + %% Check if the Behavior fields are filled from a file that meets standard Buzcode standards (*.behavior.mat) + % This flag is added from the addBehavior function + + if isempty(strfind(nwb2.processing.get('behavior').description, 'behavior.mat file')) + error('The behavior in this nwb file was not added from a file that meets the standard Buzcode format') + end + + %% If the check for standard format passed, get the behavior + % fields + + % Fill the behavior + behaviorstruct = struct; % Initialize + + % Select the key that is within the Behavior selection + allBehaviorfields = keys(nwb2.processing.get('behavior').nwbdatainterface.get(allBehaviorKeys(iBehavior)).spatialseries); + + %% Fill the Behavior Structure - This is added to the NWB file from addBehavior + + for iField = 1:length(allBehaviorfields) + + if strcmp('behaviorinfo',allBehaviorfields{iField}) + selected_behavior = nwb2.processing.get('behavior').nwbdatainterface.get(allBehaviorKeys(iBehavior)).spatialseries.get(allBehaviorfields{iField}); % This is for easier reference through the code + behaviorstruct.behavior.samplingRate = 1000 ./ nanmean(diff(selected_behavior.timestamps.load))./1000; + behaviorstruct.behavior.units = selected_behavior.data_unit; + behaviorstruct.behavior.rotationType = ''; + + behaviorstruct.behavior.timestamps = selected_behavior.timestamps.load; + + behaviorstruct.behavior.behaviorinfo.description = selected_behavior.description; + behaviorstruct.behavior.behaviorinfo.substructnames = {allBehaviorfields{~ismember(allBehaviorfields, 'behaviorinfo')}}; + behaviorstruct.behavior.behaviorinfo.errorPerMarker = selected_behavior.data.load; + behaviorstruct.behavior.behaviorinfo.acquisitionsystem = selected_behavior.control_description.load; + + else + selected_behavior = nwb2.processing.get('behavior').nwbdatainterface.get(allBehaviorKeys(iBehavior)).spatialseries.get(allBehaviorfields{iField}); % This is for easier reference through the code + + axis = {'x','y','z','w'}; % I assume that this is the order that the axis are filled + + nAxis = selected_behavior.data.dims(1); % 4 x 263958 for orientation + + for iAxis = 1:nAxis + behaviorstruct.behavior.(allBehaviorfields{iField}).(axis{iAxis}) = selected_behavior.data.load([1,iAxis],[Inf,iAxis]); + end + end + end + + %% Fill the Events - Trials subStructure - This is added to the NWB file from addTrials + trials_object = nwb2.processing.get('behavior').nwbdatainterface.get(['trials_' allBehaviorKeys{iBehavior}]); % Time intervals object + + % SECOND CHECK THAT THE TRIALS FIELD HAS BEEN FILLED FROM A + % BEHAVIOR.MAT FILE + if isempty(strfind(trials_object.description, 'behavior.mat file')) + error('The trials in this nwb file was not added from a file that meets the standard Buzcode format') + end + + nTrials = trials_object.start_time.data.dims; + + for iTrial = 1:nTrials + for iColname = 1:length(trials_object.colnames) + + singleTrialData = trials_object.vectordata.get('trials').data.load([1,iTrial,iColname],[Inf,iTrial,iColname]); + singleTrialData = singleTrialData(~isnan(singleTrialData)); % I concatenated with NaNs in the addTrials function in order to use a matrix and not a cell array + + if isempty(strfind(trials_object.colnames{iColname},'orientation')) + behaviorstruct.behavior.events.trials{iTrial}.(trials_object.colnames{iColname}) = singleTrialData; + else + behaviorstruct.behavior.events.trials{iTrial}.orientation.(trials_object.colnames{iColname}(end)) = singleTrialData; + end + end + behaviorstruct.behavior.events.trials{iTrial}.direction = trials_object.vectordata.get('direction').data.load(iTrial); + behaviorstruct.behavior.events.trials{iTrial}.type = trials_object.vectordata.get('type').data.load(iTrial); + end + + + for iMap = 1: trials_object.vectordata.get('map').data.dims(2) % Number of conditions + axis = {'x','y','z'}; + for iAxis = 1:length(axis) + behaviorstruct.behavior.events.map{iMap}.(axis{iAxis}) = trials_object.vectordata.get('map').data.load([1,iMap,iAxis], [Inf,iMap,iAxis]); + end + end + + behaviorstruct.behavior.events.trialIntervals = [trials_object.start_time trials_object.stop_time]; + behaviorstruct.behavior.events.conditionType = trials_object.vectordata.get('conditionType').data.load; + + %% Optional + % reorder for easy comparison with the tutorial example + C = {'position','timestamps','samplingRate','units','orientation','rotationType','events','behaviorinfo'}; + behaviorstruct.behavior = orderfields(behaviorstruct.behavior,C); + + + else + disp('No behavior in this NWB file') + behaviorstruct = []; + return + end + + else + disp('No .behavior.mat or NWB files were found') + behaviorstruct = []; + return + end + + end + + +end \ No newline at end of file diff --git a/io/bz_LoadCellinfo.m b/io/bz_LoadCellinfo.m index cb604440..42c4c313 100755 --- a/io/bz_LoadCellinfo.m +++ b/io/bz_LoadCellinfo.m @@ -28,6 +28,8 @@ % filename filename loaded % %DLevenstein 2018 + +% Added support for NWB: Konstantinos Nasiotis 2019 %% p = inputParser; addParameter(p,'dataset',false) @@ -104,13 +106,19 @@ if ~exist('cellinfoName','var') || isempty(cellinfoName) allCellinfoFiles = dir(fullfile(basePath,[baseName,'.','*','.cellinfo.mat'])); - [s,v] = listdlg('PromptString','Which cellinfo.mat would you like to load?',... - 'ListString',{allCellinfoFiles.name},'SelectionMode','single'); - if isempty(s) - cellinfo = []; filename = []; - return + + if ~isempty(allCellinfoFiles) + [s,v] = listdlg('PromptString','Which cellinfo.mat would you like to load?',... + 'ListString',{allCellinfoFiles.name},'SelectionMode','single'); + if isempty(s) + cellinfo = []; filename = []; + return + end + filename = fullfile(basePath,allCellinfoFiles(s).name); + else + filename = []; end - filename = fullfile(basePath,allCellinfoFiles(s).name); + else filename = fullfile(basePath,[baseName,'.',cellinfoName,'.cellinfo.mat']); end @@ -118,10 +126,14 @@ if exist(filename,'file') cellinfostruct = load(filename); -else - warning([filename,' does not exist...']) - cellinfo = []; - return +else % Check for NWB file + + cellinfostruct = loadCellInfoNWB(basePath); + + if isempty(cellinfostruct) + cellinfo = []; + return + end end varsInFile = fieldnames(cellinfostruct); @@ -140,6 +152,80 @@ warning('Your cellinfo structure does not meet buzcode standards. Sad.') end + + + +%% Function that loads NWB fields to a Behavior structure +function cellinfostruct = loadCellInfoNWB(basePath) + + % Check for NWB file since no .behavior.mat files were found + d = dir(fullfile(basePath, '*.nwb')); + if length(d) > 1 % If more than one .nwb files exist in the directory, select which one to load from + warning('there is more than one .nwb file in this directory'); + [iNWBFile, ~] = listdlg('PromptString','Which NWB file would you like to load?',... + 'ListString',{d.name},'SelectionMode','single'); + d = d(iNWBFile); + elseif length(d) == 0 + d = dir('*nwb'); + if isempty(d) + disp('No .cellinfo.mat or NWB files present in this directory') + cellinfostruct = []; + return + end + end + + % If an NWB file exists, load the info according the Buzcode standard + nwb_file = fullfile(d.folder, d.name); + + nwb2 = nwbRead(nwb_file); + + if ~isempty(nwb2.units) + + nNeurons = length(nwb2.units.id.data.load); + + cellinfostruct.spikes.UID = double(nwb2.units.id.data.load'); + cellinfostruct.spikes.shankID = nwb2.units.vectordata.get('shankID').data.load'; + cellinfostruct.spikes.cluID = nwb2.units.vectordata.get('cluID').data.load'; + cellinfostruct.spikes.maxWaveformCh = nwb2.units.vectordata.get('maxWaveformCh').data.load'; + cellinfostruct.spikes.region = nwb2.units.vectordata.get('region').data.load'; + cellinfostruct.spikes.sessionName = nwb2.identifier; + + for iNeuron = 1:nNeurons + cellinfostruct.spikes.rawWaveform{1,iNeuron} = nwb2.units.waveform_mean.data.load([iNeuron, 1], [iNeuron,Inf]); + + if iNeuron == 1 + cellinfostruct.spikes.times{1,iNeuron} = nwb2.units.spike_times.data.load(1, double(nwb2.units.spike_times_index.data.load(iNeuron))); + else + cellinfostruct.spikes.times{1,iNeuron} = nwb2.units.spike_times.data.load(double(nwb2.units.spike_times_index.data.load(iNeuron-1))+1, double(nwb2.units.spike_times_index.data.load(iNeuron))); + end + end + + %% Optional + % reorder for easy comparison with the tutorial example + C = {'UID','times','shankID','cluID','rawWaveform','maxWaveformCh','region','sessionName'}; + cellinfostruct.spikes = orderfields(cellinfostruct.spikes,C); + + else + disp('No units information in the selected NWB file'); + cellinfostruct = []; + return + end + +end + + + + + + + + + + + + + + end diff --git a/io/bz_LoadEvents.m b/io/bz_LoadEvents.m index 32219a5d..bb1deffa 100755 --- a/io/bz_LoadEvents.m +++ b/io/bz_LoadEvents.m @@ -8,6 +8,8 @@ %Future update: 'all' (nonfunctional) to load all events.mat files for a given recording. % %DLevenstein 2017 + +% Added NWB support: Konstantinos Nasiotis 2019 %% if ~exist('basePath','var') basePath = pwd; @@ -16,13 +18,19 @@ if ~exist('eventsName','var') allEventsFiles = dir(fullfile(basePath,[baseName,'.','*','.events.mat'])); - [s,v] = listdlg('PromptString','Which events.mat would you like to load?',... - 'ListString',{allEventsFiles.name},'SelectionMode','single'); - if isempty(s) - events = []; filename = []; - return + + if ~isempty(allEventsFiles) + [s,v] = listdlg('PromptString','Which events.mat would you like to load?',... + 'ListString',{allEventsFiles.name},'SelectionMode','single'); + if isempty(s) + events = []; filename = []; + return + end + filename = fullfile(basePath,allEventsFiles(s).name); + else + check_for_NWB = 1; + filename = []; end - filename = fullfile(basePath,allEventsFiles(s).name); else filename = fullfile(basePath,[baseName,'.',eventsName,'.events.mat']); end @@ -30,12 +38,24 @@ if exist(filename,'file') eventstruct = load(filename); + check_for_NWB = 0; else - warning([filename,' does not exist...']) - events = []; - return + check_for_NWB = 1; +end + + +if check_for_NWB + if ~exist('eventsName','var') + eventsName = []; + end + eventstruct = loadEventsNWB(basePath, eventsName); + if isempty(eventstruct) + events = []; filename = []; + return + end end + varsInFile = fieldnames(eventstruct); if numel(varsInFile)==1 @@ -52,5 +72,93 @@ warning('Your events structure does not meet buzcode standards. Sad.') end + + +%% Function that loads NWB fields to an events structure +function eventstruct = loadEventsNWB(basePath, eventsName) + % Check for NWB file since no .events.mat files were found + d = dir(fullfile(basePath, '*.nwb')); + if length(d) > 1 % If more than one .nwb files exist in the directory, select which one to load from + warning('there is more than one .nwb file in this directory'); + [iNWBFile, ~] = listdlg('PromptString','Which NWB file would you like to load?',... + 'ListString',{d.name},'SelectionMode','single'); + d = d(iNWBFile); + elseif length(d) == 0 + d = dir('*nwb'); + if isempty(d) + warning('No .events.mat or NWB files present in this directory') + eventstruct = []; + return + end + end + nwb_file = fullfile(d.folder, d.name); + + if exist(nwb_file,'file') == 2 + nwb2 = nwbRead(nwb_file); + + % Check if the events field exists in the dataset + events_key_exists = sum(ismember(keys(nwb2.processing),'events')); + + if events_key_exists + all_event_keys = keys(nwb2.processing.get('events').nwbdatainterface); + else + disp('No events in this .nwb file') + eventstruct = []; + return + end + + if isempty(all_event_keys) + disp('No events in this .nwb file') + eventstruct = []; + return + + else + % Check if a specific event was called to be loaded. If not, display a + % pop-up list with the available events for selection + if isempty(eventsName) + [iEvent, ~] = listdlg('PromptString','Which event type would you like to load?',... + 'ListString',all_event_keys,'SelectionMode','single'); + else + if sum(ismember(all_event_keys, eventsName))==1 % If the event exists, find its index + iEvent = find(ismember(all_event_keys, eventsName)); + else + disp('The selected event doesnt exist in this NWB file') + eventstruct = []; + return + end + end + + %% Get the selected event from the NWB file + eventstruct = struct; + + event_module = nwb2.processing.get('events').nwbdatainterface.get(all_event_keys{iEvent}); % For convenience + % Check if the timestamps are start-stop or just start + % For information about the way this is coded, check: + % https://nwb-schema.readthedocs.io/en/latest/format.html#sec-intervalseries + first_Element = event_module.data.load(1); + + if ischar(first_Element) % This means I only have event start data + eventstruct.events.timestamps = event_module.timestamps.load;% neuroscope compatible matrix with 1 column - [starts] (in seconds) + else + pointers = double(event_module.data.load)'; % positive values imply start, negative imply stop + eventstruct.events.timestamps = [event_module.timestamps.load(pointers>0)' event_module.timestamps.load(pointers<0)'];% neuroscope compatible matrix with 2 columns - [starts stops] (in seconds) + end + + + eventstruct.events.detectorinfo.detectorname = 'N/A';% substructure with information about the detection method (fields below) + eventstruct.events.detectorinfo.detectionparms = 'N/A';% parameters used for detection + eventstruct.events.detectorinfo.detectiondate = 'N/A';% date of detection + eventstruct.events.detectorinfo.detectionintervals = 'N/A';% [start stop] pairs of intervals used for detection (optional) +% eventstruct.events.detectorinfo.detectionchannel = ;% channel used for detection (optional) + end + + else + disp('No .events.mat or NWB files were found') + eventstruct = []; + return + end + end + +end \ No newline at end of file