Note
Go to the end to download the full example code.
Reading files with neo.rawio#
compare with read_files_neo_io.py
First we import a RawIO from neo.rawio For this example we will use PlexonRawIO
import urllib
from neo.rawio import PlexonRawIO
Get Plexon files
# We will be pulling these files down, but if you have a local file
# then all you need to do is specify the file location on your
# computer. NeuralEnsemble keeps a wide variety of freely accesible, small
# test files that can be used. So for this example we will take advantage of that
# fact.
url_repo = "https://web.gin.g-node.org/NeuralEnsemble/ephy_testing_data/raw/master/"
distantfile = url_repo + "plexon/File_plexon_3.plx"
localfile = "File_plexon_3.plx"
urllib.request.urlretrieve(distantfile, localfile)
('File_plexon_3.plx', <http.client.HTTPMessage object at 0x7f5b01fb58d0>)
Create a reader
# All it takes to create a reader is giving the filename (or dirname)
# Then we need to do the slow step of parsing the header with the
# `parse_header` function. This collects metadata as well as
# make all future steps much faster for us
reader = PlexonRawIO(filename="File_plexon_3.plx")
reader.parse_header()
print(reader)
# we can view metadata in the header
print(reader.header)
PlexonRawIO: File_plexon_3.plx
nb_block: 1
nb_segment: [1]
signal_streams: [Signals 0 (chans: 1)]
signal_channels: [V1]
spike_channels: [Wspk1u, Wspk2u, Wspk4u, Wspk5u ... Wspk29u , Wspk30u , Wspk31u , Wspk32u]
event_channels: []
{'nb_block': 1, 'nb_segment': [1], 'signal_streams': array([('Signals 0', '0')], dtype=[('name', '<U64'), ('id', '<U64')]), 'signal_channels': array([('V1', '0', 1000., 'int16', '', 2.44140625, 0., '0')],
dtype=[('name', '<U64'), ('id', '<U64'), ('sampling_rate', '<f8'), ('dtype', '<U16'), ('units', '<U64'), ('gain', '<f8'), ('offset', '<f8'), ('stream_id', '<U64')]), 'spike_channels': array([('Wspk1u', 'ch1#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk2u', 'ch2#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk4u', 'ch3#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk5u', 'ch4#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk6u', 'ch5#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk7u', 'ch6#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk8u', 'ch7#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk9u', 'ch8#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk10u', 'ch9#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk11u', 'ch10#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk12u', 'ch11#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk13u', 'ch12#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk14u', 'ch13#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk15u', 'ch14#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk16u', 'ch15#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk17u', 'ch16#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk19u', 'ch18#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk20u', 'ch19#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk21u', 'ch20#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk22u', 'ch21#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk23u', 'ch22#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk24u', 'ch23#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk25u', 'ch24#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk26u', 'ch25#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk27u', 'ch26#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk28u', 'ch27#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk29u', 'ch28#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk30u', 'ch29#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk31u', 'ch30#0', '', 7.32421875e-05, 0., -1, 30000.),
('Wspk32u', 'ch31#0', '', 7.32421875e-05, 0., -1, 30000.)],
dtype=[('name', '<U64'), ('id', '<U64'), ('wf_units', '<U64'), ('wf_gain', '<f8'), ('wf_offset', '<f8'), ('wf_left_sweep', '<i8'), ('wf_sampling_rate', '<f8')]), 'event_channels': array([], dtype=[('name', '<U64'), ('id', '<U64'), ('type', 'S5')])}
Read signal chunks This is how we read raw data. We choose indices that we want or we can use None to mean look at all channels. We also need to specify the block of data (block_index) as well as the segment (seg_index). Then we give the index start and stop. Since we often think in time: to go from time to index would just require the sample rate (so index = time / sampling_rate)
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
)
# raw_sigs are not voltages so to convert to voltages we do the follwing
float_sigs = reader.rescale_signal_raw_to_float(raw_sigs, dtype="float64")
# We can see that the shapes are the same, but that the datatypes
# are different once we've rescaled our data
print("Raw Data: ", raw_sigs.shape, raw_sigs.dtype)
print("Scaled Data: ", float_sigs.shape, float_sigs.dtype)
Raw Data: (1024, 1) int16
Scaled Data: (1024, 1) float64
Each rawio gives you access to lots of information about your data some of this information comes from functions other information is stored as metadata in the reader.header
sampling_rate = reader.get_signal_sampling_rate()
# Like above we need to indicate the block and segment
t_start = reader.get_signal_t_start(block_index=0, seg_index=0)
units = reader.header["signal_channels"][0]["units"]
# and we can display all of this information
print(f"{sampling_rate=}, {t_start=}, {units=}")
sampling_rate=1000.0, t_start=0.0, units=''
Some rawio’s and file formats also provide information about spikes If a rawio can’t read this data it will raise an error, but lucky for us PlexonRawIO does have spikes data!!
# Count units and spikes per unit
nb_unit = reader.spike_channels_count()
print(f"nb_unit: {nb_unit}\n") # nb_unit stands for number of units
print("spike_channel_index nb_spike")
for spike_channel_index in range(nb_unit):
nb_spike = reader.spike_count(block_index=0, seg_index=0, spike_channel_index=spike_channel_index)
print(f"{spike_channel_index}: {nb_spike}\n")
# Read spike times and rescale (just like analogsignal above)
spike_timestamps = reader.get_spike_timestamps(
block_index=0, seg_index=0, spike_channel_index=0, t_start=0.0, t_stop=10.0
)
print(f"{spike_timestamps.shape=}\n{spike_timestamps.dtype=}\n{spike_timestamps[:5]=}\n")
spike_times = reader.rescale_spike_timestamp(spike_timestamps, dtype="float64")
print(f"{spike_times.shape=}\n{spike_times.dtype=}\n{spike_times[:5]}\n")
nb_unit: 30
spike_channel_index nb_spike
0: 701
1: 716
2: 69
3: 12
4: 95
5: 37
6: 25
7: 15
8: 33
9: 392
10: 337
11: 333
12: 85
13: 13
14: 20
15: 750
16: 3
17: 1345
18: 826
19: 818
20: 744
21: 753
22: 34
23: 36
24: 30
25: 19
26: 2
27: 42
28: 15
29: 51
spike_timestamps.shape=(424,)
spike_timestamps.dtype=dtype('int64')
spike_timestamps[:5]=array([ 90, 420, 708, 1020, 1310])
spike_times.shape=(424,)
spike_times.dtype=dtype('float64')
[0.003 0.014 0.0236 0.034 0.04366667]
Some file formats can also give waveform information. We are lucky again our file has waveform data!! We forms are a 3d dataset of (nb_spike, nb_channel, nb_sample)
# Read spike waveforms
raw_waveforms = reader.get_spike_raw_waveforms(
block_index=0, seg_index=0, spike_channel_index=0, t_start=0.0, t_stop=10.0
)
print(f"{raw_waveforms.shape=}\n{raw_waveforms.dtype=}\n{raw_waveforms[0, 0, :4]=}\n")
float_waveforms = reader.rescale_waveforms_to_float(raw_waveforms, dtype="float32", spike_channel_index=0)
print(f"{float_waveforms.shape=}\n{float_waveforms.dtype=}{float_waveforms[0, 0, :4]=}\n")
raw_waveforms.shape=(424, 1, 64)
raw_waveforms.dtype=dtype('int16')
raw_waveforms[0, 0, :4]=array([-449, -206, 34, 40], dtype=int16)
float_waveforms.shape=(424, 1, 64)
float_waveforms.dtype=dtype('float32')float_waveforms[0, 0, :4]=array([-0.03288575, -0.01508789, 0.00249023, 0.00292969], dtype=float32)
RawIOs can also read event timestamps. But looks like our luck ran out let’s grab a new file to see this feature
# Read event timestamps and times (take another file)
distantfile = url_repo + "plexon/File_plexon_2.plx"
localfile = "File_plexon_2.plx"
urllib.request.urlretrieve(distantfile, localfile)
('File_plexon_2.plx', <http.client.HTTPMessage object at 0x7f5b01f91710>)
Since this is a new file we need to read initialize our reader as well as parse the header
reader = PlexonRawIO(filename="File_plexon_2.plx")
reader.parse_header()
# if we look at this header we see it is different than the header above
print(reader.header)
{'nb_block': 1, 'nb_segment': [1], 'signal_streams': array([], dtype=[('name', '<U64'), ('id', '<U64')]), 'signal_channels': array([],
dtype=[('name', '<U64'), ('id', '<U64'), ('sampling_rate', '<f8'), ('dtype', '<U16'), ('units', '<U64'), ('gain', '<f8'), ('offset', '<f8'), ('stream_id', '<U64')]), 'spike_channels': array([('sig001', 'ch1#1', '', 0.00073242, 0., -1, 0.),
('sig001', 'ch1#2', '', 0.00073242, 0., -1, 0.),
('sig001', 'ch1#3', '', 0.00073242, 0., -1, 0.),
('sig002', 'ch2#1', '', 0.00073242, 0., -1, 0.),
('sig002', 'ch2#2', '', 0.00073242, 0., -1, 0.),
('sig002', 'ch2#3', '', 0.00073242, 0., -1, 0.)],
dtype=[('name', '<U64'), ('id', '<U64'), ('wf_units', '<U64'), ('wf_gain', '<f8'), ('wf_offset', '<f8'), ('wf_left_sweep', '<i8'), ('wf_sampling_rate', '<f8')]), 'event_channels': array([('Event001', '1', b'event'), ('Event002', '2', b'event'),
('Event003', '3', b'event'), ('Event004', '4', b'event'),
('Event005', '5', b'event'), ('Event006', '6', b'event'),
('Event007', '7', b'event'), ('Event008', '8', b'event'),
('Event009', '9', b'event'), ('Event010', '10', b'event'),
('Event011', '11', b'event'), ('Event012', '12', b'event'),
('Event013', '13', b'event'), ('Event014', '14', b'event'),
('Event015', '15', b'event'), ('Event016', '16', b'event'),
('Strobed', '257', b'event'), ('Start\x0018', '258', b'event'),
('Stop\x00019', '259', b'event'), ('Keyboard1', '101', b'event'),
('Keyboard2', '102', b'event'), ('Keyboard3', '103', b'event'),
('Keyboard4', '104', b'event'), ('Keyboard5', '105', b'event'),
('Keyboard6', '106', b'event'), ('Keyboard7', '107', b'event'),
('Keyboard8', '108', b'event'), ('Keyboard9', '109', b'event')],
dtype=[('name', '<U64'), ('id', '<U64'), ('type', 'S5')])}
Now let’s look at some event data. This could be things like stimuli applied during the course of an ephys recording
nb_event_channel = reader.event_channels_count()
print(f"nb_event_channel: {nb_event_channel}")
# now iterate through the channels
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(f"chan_index: {chan_index} nb_event: {nb_event}\n")
nb_event_channel: 28
chan_index: 0 nb_event: 1
chan_index: 1 nb_event: 0
chan_index: 2 nb_event: 0
chan_index: 3 nb_event: 0
chan_index: 4 nb_event: 0
chan_index: 5 nb_event: 0
chan_index: 6 nb_event: 0
chan_index: 7 nb_event: 0
chan_index: 8 nb_event: 0
chan_index: 9 nb_event: 0
chan_index: 10 nb_event: 0
chan_index: 11 nb_event: 0
chan_index: 12 nb_event: 0
chan_index: 13 nb_event: 0
chan_index: 14 nb_event: 0
chan_index: 15 nb_event: 0
chan_index: 16 nb_event: 0
chan_index: 17 nb_event: 1
chan_index: 18 nb_event: 0
chan_index: 19 nb_event: 0
chan_index: 20 nb_event: 0
chan_index: 21 nb_event: 0
chan_index: 22 nb_event: 0
chan_index: 23 nb_event: 0
chan_index: 24 nb_event: 0
chan_index: 25 nb_event: 0
chan_index: 26 nb_event: 0
chan_index: 27 nb_event: 0
Finally we can get our actual event timestamps. Some file formats provide the real timestamps (timestamps in s/ms) others have raw timestamps (in samples) so we can do the same style of functions as above. Get the raw timestamps and convert to real times with a rescale function.
ev_timestamps, ev_durations, ev_labels = reader.get_event_timestamps(
block_index=0, seg_index=0, event_channel_index=0, t_start=None, t_stop=None
)
print(f"{ev_timestamps=}\n{ev_durations=}\n{ev_labels=}\n")
ev_times = reader.rescale_event_timestamp(ev_timestamps, dtype="float64")
print(f"{ev_times=}\n")
ev_timestamps=array([1268])
ev_durations=None
ev_labels=array(['0'], dtype='<U5')
ev_times=array([0.0317])
Total running time of the script: (0 minutes 46.768 seconds)