Neo RawIO

For performance and memory consumption reasons a new layer has been added to Neo.

In brief:
  • neo.io is the user-oriented read/write layer. Reading consists of getting a tree of Neo objects from a data source (file, url, or directory). When reading, all Neo objects are correctly scaled to the correct units. Writing consists of making a set of Neo objects persistent in a file format.
  • neo.rawio is a low-level layer for reading data only. Reading consists of getting NumPy buffers (often int16/int64) of signals/spikes/events. Scaling to real values (microV, times, …) is done in a second step. Here the underlying objects must be consistent across Blocks and Segments for a given data source.

The neo.rawio API has been added for developers. The neo.rawio is close to what could be a C API for reading data but in Python/NumPy.

Not all IOs are implemented in neo.rawio but all classes implemented in neo.rawio are also available in neo.io.

Possible uses of the neo.rawio API are:
  • fast reading chunks of signals in int16 and do the scaling of units (uV) on a GPU while scaling the zoom. This should improve bandwith HD to RAM and RAM to GPU memory.
  • load only some small chunk of data for heavy computations. For instance the spike sorting module tridesclous does this.
The neo.rawio API is less flexible than neo.io and has some limitations:
  • read-only
  • AnalogSignals must have the same characteristcs across all Blocks and Segments: sampling_rate, shape[1], dtype
  • AnalogSignals should all have the same value of sampling_rate, otherwise they won’t be read at the same time.
  • Units must have SpikeTrain event if empty across all Block and Segment
  • Epoch and Event are processed the same way (with durations=None for Event).
For an intuitive comparison of neo.io and neo.rawio see:
  • example/read_file_neo_io.py
  • example/read_file_neo_rawio.py

One speculative benefit of the neo.rawio API should be that a developer should be able to code a new RawIO class with little knowledge of the Neo tree of objects or of the quantities package.

Basic usage

First create a reader from a class:

>>> from neo.rawio import PlexonRawIO
>>> reader = PlexonRawIO(filename='File_plexon_3.plx')

Then browse the internal header and display information:

>>> reader.parse_header()
>>> print(reader)
PlexonRawIO: File_plexon_3.plx
nb_block: 1
nb_segment:  [1]
signal_channels: [V1]
unit_channels: [Wspk1u, Wspk2u, Wspk4u, Wspk5u ... Wspk29u Wspk30u Wspk31u Wspk32u]
event_channels: []

You get the number of blocks and segments per block. You have information about channels: signal_channels, unit_channels, event_channels.

All this information is internally available in the header dict:

>>> for k, v in reader.header.items():
...    print(k, v)
signal_channels [('V1', 0,  1000., 'int16', '',  2.44140625,  0., 0)]
event_channels []
nb_segment [1]
nb_block 1
unit_channels [('Wspk1u', 'ch1#0', '',  0.00146484,  0., 0,  30000.)
('Wspk2u', 'ch2#0', '',  0.00146484,  0., 0,  30000.)
...

Read signal chunks of data and scale them:

>>> channel_indexes = None  #could be channel_indexes = [0]
>>> raw_sigs = reader.get_analogsignal_chunk(block_index=0, seg_index=0,
                    i_start=1024, i_stop=2048, channel_indexes=channel_indexes)
>>> float_sigs = reader.rescale_signal_raw_to_float(raw_sigs, dtype='float64')
>>> sampling_rate = reader.get_signal_sampling_rate()
>>> t_start = reader.get_signal_t_start(block_index=0, seg_index=0)
>>> units =reader.header['signal_channels'][0]['units']
>>> print(raw_sigs.shape, raw_sigs.dtype)
>>> print(float_sigs.shape, float_sigs.dtype)
>>> print(sampling_rate, t_start, units)
(1024, 1) int16
(1024, 1) float64
1000.0 0.0 V

There are 3 ways to select a subset of channels: by index (0 based), by id or by name. By index is not ambiguous 0 to n-1 (included), for some IOs channel_names (and sometimes channel_ids) have no guarantees to be unique, in such cases it would raise an error.

Example with BlackrockRawIO for the file FileSpec2.3001:

>>> raw_sigs = reader.get_analogsignal_chunk(channel_indexes=None) #Take all channels
>>> raw_sigs1 = reader.get_analogsignal_chunk(channel_indexes=[0,  2, 4])) #Take 0 2 and 4
>>> raw_sigs2 = reader.get_analogsignal_chunk(channel_ids=[1, 3, 5]) # Same but with there id (1 based)
>>> raw_sigs3 = reader.get_analogsignal_chunk(channel_names=['chan1', 'chan3', 'chan5'])) # Same but with there name
print(raw_sigs1.shape[1], raw_sigs2.shape[1], raw_sigs3.shape[1])
3, 3, 3

Inspect units channel. Each channel gives a SpikeTrain for each Segment. Note that for many formats a physical channel can have several units after spike sorting. So the nb_unit could be more than physical channel or signal channels.

>>> nb_unit = reader.unit_channels_count()
>>> print('nb_unit', nb_unit)
nb_unit 30
>>> for unit_index in range(nb_unit):
...     nb_spike = reader.spike_count(block_index=0, seg_index=0, unit_index=unit_index)
...     print('unit_index', unit_index, 'nb_spike', nb_spike)
unit_index 0 nb_spike 701
unit_index 1 nb_spike 716
unit_index 2 nb_spike 69
unit_index 3 nb_spike 12
unit_index 4 nb_spike 95
unit_index 5 nb_spike 37
unit_index 6 nb_spike 25
unit_index 7 nb_spike 15
unit_index 8 nb_spike 33
...

Get spike timestamps only between 0 and 10 seconds and convert them to spike times:

>>> spike_timestamps = reader.spike_timestamps(block_index=0, seg_index=0, unit_index=0,
                    t_start=0., t_stop=10.)
>>> print(spike_timestamps.shape, spike_timestamps.dtype, spike_timestamps[:5])
(424,) int64 [  90  420  708 1020 1310]
>>> spike_times =  reader.rescale_spike_timestamp( spike_timestamps, dtype='float64')
>>> print(spike_times.shape, spike_times.dtype, spike_times[:5])
(424,) float64 [ 0.003       0.014       0.0236      0.034       0.04366667]

Get spike waveforms between 0 and 10 s:

>>> raw_waveforms = reader.spike_raw_waveforms(  block_index=0, seg_index=0, unit_index=0,
                    t_start=0., t_stop=10.)
>>> print(raw_waveforms.shape, raw_waveforms.dtype, raw_waveforms[0,0,:4])
(424, 1, 64) int16 [-449 -206   34   40]
>>> float_waveforms = reader.rescale_waveforms_to_float(raw_waveforms, dtype='float32', unit_index=0)
>>> print(float_waveforms.shape, float_waveforms.dtype, float_waveforms[0,0,:4])
(424, 1, 64) float32 [-0.65771484 -0.30175781  0.04980469  0.05859375]

Count events per channel:

>>> reader = PlexonRawIO(filename='File_plexon_2.plx')
>>> reader.parse_header()
>>> nb_event_channel = reader.event_channels_count()
nb_event_channel 28
>>> print('nb_event_channel', nb_event_channel)
>>> for chan_index in range(nb_event_channel):
...     nb_event = reader.event_count(block_index=0, seg_index=0, event_channel_index=chan_index)
...     print('chan_index',chan_index, 'nb_event', nb_event)
chan_index 0 nb_event 1
chan_index 1 nb_event 0
chan_index 2 nb_event 0
chan_index 3 nb_event 0
...

Read event timestamps and times for chanindex=0 and with time limits (t_start=None, t_stop=None):

>>> ev_timestamps, ev_durations, ev_labels = reader.event_timestamps(block_index=0, seg_index=0, event_channel_index=0,
                    t_start=None, t_stop=None)
>>> print(ev_timestamps, ev_durations, ev_labels)
[1268] None ['0']
>>> ev_times = reader.rescale_event_timestamp(ev_timestamps, dtype='float64')
>>> print(ev_times)
[ 0.0317]