""" This script loads a complete session in the blackrock format and converts it to a single nix file """ import os import numpy as np import quantities as pq import neo from neo.test.tools import (assert_same_sub_schema, assert_same_annotations, assert_same_array_annotations) from elephant.signal_processing import butter from reachgraspio import reachgraspio # Choose which session you want to convert into a nix file session = "i140703-001" # session = "l101210-001" # Input data. i.e., original Blackrock files and odML dataset_dir = '../datasets_blackrock' session_path = f'{dataset_dir}/{session}' odml_dir = os.path.join('..', 'datasets_blackrock') # Output for the nix files nix_dataset_dir = '../datasets_nix' nix_session_path = f'{nix_dataset_dir}/{session}' ##### LOAD BLACKROCK FILES ############ session = reachgraspio.ReachGraspIO( filename=session_path, odml_directory=odml_dir, verbose=False) block = session.read_block(lazy=False, load_waveforms=True) # ============================================================================= # Create offline filtered LFP # # Here, we construct one offline filtered LFP from each ns5 (monkey L) or ns6 # (monkey N) raw recording trace. For monkey N, this filtered LFP can be # compared to the LFPs in the ns2 file (note that monkey L contains only # behavioral signals in the ns2 file). Also, we assign telling names to each # Neo AnalogSignal, which is used for plotting later on in this script. # ============================================================================= nsx_to_anasig_name = {2: 'LFP signal (online filtered)', 5: 'raw signal', 6: 'raw signal'} # this factor was experimentally determined as an approximate shift introduced # by the online filtering. Here, we integrate this shift such that the offline # filtered signal is aligned to the online filtered signal (if that were # available) time_shift_factor = -0.42*pq.ms filtered_anasig = None raw_anasig = None # identify neuronal signals and provide labels for plotting for anasig in block.segments[0].analogsignals: # skip non-neuronal signals if not anasig.annotations['neural_signal']: continue # identify nsx source of signals in this AnalogSignal object if 'nsx' in anasig.annotations: nsx = anasig.annotations['nsx'] else: nsx = np.unique(anasig.array_annotations['nsx']) assert len(nsx) == 1, 'Different nsx sources in AnalogSignal' nsx = nsx[0] if nsx == 2: # AnalogSignal is LFP from ns2 filtered_anasig = anasig elif nsx in [5, 6]: # AnalogSignal is raw signal from ns5 or ns6 raw_anasig = anasig if filtered_anasig is None: print("Filtering raw time series to obtain LFP") for f in range(raw_anasig.shape[1]): # filtering must be done channel by channel for memory reasons (requires approx. 32 GB RAM) print(f"Processing channel {f}") filtered_signal = butter( raw_anasig[:, f], highpass_freq=None, lowpass_freq=250 * pq.Hz, filter_function='sosfiltfilt', order=4) downsampled_signal=filtered_signal.downsample(30).time_shift(time_shift_factor) # first run? Create a new Analogsignal if f == 0: offline_filtered_anasig = neo.AnalogSignal( np.zeros((downsampled_signal.shape[0], raw_anasig.shape[1])) *\ downsampled_signal.units, t_start=downsampled_signal.t_start, sampling_rate=downsampled_signal.sampling_rate) offline_filtered_anasig[:, f] = downsampled_signal if 'nsx' in anasig.annotations: nsx = anasig.annotations['nsx'] else: nsx = anasig.array_annotations["nsx"][0] offline_filtered_anasig.name = f"NeuralTimeSeriesDownsampled" offline_filtered_anasig.description = "Downsampled continuous neuronal recordings, where the downsampling was " \ "performed off-line during post-processing" # Attach all offline filtered LFPs to the segment of data block.segments[0].analogsignals.append(offline_filtered_anasig) ##### SAVE NIX FILE ################### nix_filename = nix_session_path + '.nix' if os.path.exists(nix_filename): print('Nix file already exists and will not be overwritten.') else: with neo.NixIO(nix_filename) as io: print(f'Saving nix file at {nix_filename}') io.write_block(block) ##### VALIDATION OF FILE CONTENT ###### with neo.NixIO(nix_filename, mode='ro') as io: blocks = io.read_all_blocks() assert len(blocks) == 1 block_new = blocks[0] for seg_old, seg_new in zip(block.segments, block_new.segments): for anasig_old, anasig_new in zip(seg_old.analogsignals, seg_new.analogsignals): # ignoring differences in the file_origin attribute anasig_old.file_origin = anasig_new.file_origin assert_same_sub_schema(anasig_old, anasig_new) assert_same_annotations(anasig_old, anasig_new) assert_same_array_annotations(anasig_old, anasig_new) del anasig_old print(f'AnalogSignals are equivalent.') for st_old, st_new in zip(seg_old.spiketrains, seg_new.spiketrains): # ignoring differences in the file_origin attribute st_old.file_origin = st_new.file_origin assert_same_sub_schema(st_old, st_new) assert_same_annotations(st_old, st_new) assert_same_array_annotations(st_old, st_new) del st_old print(f'Spiketrains are equivalent.') for ev_old, ev_new in zip(seg_old.events, seg_new.events): # ignoring differences in the file_origin attribute ev_old.file_origin = ev_new.file_origin # ignore list-array type changes if 'color_codes' in ev_old.annotations: ev_old.annotations['color_codes'] = list(ev_old.annotations['color_codes']) assert_same_sub_schema(ev_old, ev_new) assert_same_annotations(ev_old, ev_new) assert_same_array_annotations(ev_old, ev_new) del ev_old print(f'Events are equivalent.') for ep_old, ep_new in zip(seg_old.epochs, seg_new.epochs): # ignoring differences in the file_origin attribute ep_old.file_origin = ep_new.file_origin assert_same_sub_schema(ep_old, ep_new) assert_same_annotations(ep_old, ep_new) assert_same_array_annotations(ep_old, ep_new) del ep_old print(f'Epochs are equivalent.')