Working with NWB Files using h5py and pynwb

Import the libraries we’ll need. h5py is the low-level HDF5 reader; pynwb is the NWB-aware library; numpy, pandas, and matplotlib are for working with the data once we’ve read it.

import h5py
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from uuid import uuid4

from pynwb import NWBHDF5IO

The synthetic NWB file

In the exercises in this notebook you will read from a small synthetic file, synthetic_session.nwb. Run the cell below to generate the data and construct the file. It imitates data from a single mouse subject with data from an 8-channel probe split across V1 and hippocampus.

!python generate_synthetic_nwb_data.py

Wrote synthetic_session.nwb

Background

An NWB file is a single HDF5 file with. HDF5 organises data hierarchically: groups are like folders, datasets are like files holding numeric or string arrays, and both groups and datasets can carry attributes — short key/value pairs that store metadata (units, descriptions, neurodata types). When you open an NWB file with h5py, you see this raw layout, not the NWB-typed objects. Common top-level groups are:

• /acquisition — raw recorded signals (e.g. extracellular voltage)
• /processing — derived data organised into modules (e.g. LFP, behaviour)
• /general — subject, devices, electrode metadata
• /intervals/trials — trial table
• /units — sorted-unit spike times

Opening an NWB file with h5py is a quick way to confirm what’s inside.

Exercises

Example: open the file and print a list of the top-level groups.

with h5py.File("synthetic_session.nwb", "r") as f:
    print(list(f.keys()))

['acquisition', 'analysis', 'file_create_date', 'general', 'identifier', 'intervals', 'processing', 'session_description', 'session_start_time', 'specifications', 'stimulus', 'timestamps_reference_time', 'units']

Alternative solution:

f = h5py.File("synthetic_session.nwb", "r")
print(list(f.keys()))
f.close()

['acquisition', 'analysis', 'file_create_date', 'general', 'identifier', 'intervals', 'processing', 'session_description', 'session_start_time', 'specifications', 'stimulus', 'timestamps_reference_time', 'units']

Exercise: print the list of entries inside the acquisition group.

with h5py.File("synthetic_session.nwb", "r") as f:
    print(list(f["acquisition"].keys()))

['ElectricalSeries']

Exercise: print the list of entries inside the general group.

with h5py.File("synthetic_session.nwb", "r") as f:
    print(list(f["general"].keys()))

['devices', 'experimenter', 'extracellular_ephys', 'institution', 'lab', 'session_id', 'subject']

Exercise: read and print the shape and dtype of the raw signal dataset at acquisition/ElectricalSeries/data.

with h5py.File("synthetic_session.nwb", "r") as f:
    dset = f["acquisition/ElectricalSeries/data"]
    print(dset.shape, dset.dtype)

(5000, 8) float32

Exercise: read the attributes attached to the same dataset. Which attribute tells you the physical unit of the recorded values?

with h5py.File("synthetic_session.nwb", "r") as f:
    attrs = dict(f["acquisition/ElectricalSeries/data"].attrs)
print(attrs)
# The "unit" attribute holds the physical unit (here, 'volts').

{'conversion': np.float64(1.0), 'offset': np.float64(0.0), 'resolution': np.float64(-1.0), 'unit': 'volts'}

Exercise: walk the entire file with visit and print every group and dataset path. This gives you a one-shot map of what’s stored.

with h5py.File("synthetic_session.nwb", "r") as f:
    f.visit(lambda name: print(name))

acquisition
acquisition/ElectricalSeries
acquisition/ElectricalSeries/data
acquisition/ElectricalSeries/electrodes
acquisition/ElectricalSeries/starting_time
analysis
file_create_date
general
general/devices
general/devices/probe-A
general/experimenter
general/extracellular_ephys
general/extracellular_ephys/electrodes
general/extracellular_ephys/electrodes/filtering
general/extracellular_ephys/electrodes/group
general/extracellular_ephys/electrodes/group_name
general/extracellular_ephys/electrodes/id
general/extracellular_ephys/electrodes/location
general/extracellular_ephys/electrodes/x
general/extracellular_ephys/electrodes/y
general/extracellular_ephys/electrodes/z
general/extracellular_ephys/shank-HPC
general/extracellular_ephys/shank-V1
general/institution
general/lab
general/session_id
general/subject
general/subject/age
general/subject/description
general/subject/sex
general/subject/species
general/subject/subject_id
identifier
intervals
intervals/trials
intervals/trials/id
intervals/trials/start_time
intervals/trials/stim_type
intervals/trials/stop_time
processing
processing/behavior
processing/behavior/BehavioralTimeSeries
processing/behavior/BehavioralTimeSeries/pupil_diameter
processing/behavior/BehavioralTimeSeries/pupil_diameter/data
processing/behavior/BehavioralTimeSeries/pupil_diameter/starting_time
processing/behavior/BehavioralTimeSeries/running_speed
processing/behavior/BehavioralTimeSeries/running_speed/data
processing/behavior/BehavioralTimeSeries/running_speed/starting_time
processing/ecephys
processing/ecephys/LFP
processing/ecephys/LFP/lfp
processing/ecephys/LFP/lfp/data
processing/ecephys/LFP/lfp/starting_time
session_description
session_start_time
specifications
specifications/core
specifications/core/2.9.0
specifications/core/2.9.0/namespace
specifications/core/2.9.0/nwb.base
specifications/core/2.9.0/nwb.behavior
specifications/core/2.9.0/nwb.device
specifications/core/2.9.0/nwb.ecephys
specifications/core/2.9.0/nwb.epoch
specifications/core/2.9.0/nwb.file
specifications/core/2.9.0/nwb.icephys
specifications/core/2.9.0/nwb.image
specifications/core/2.9.0/nwb.misc
specifications/core/2.9.0/nwb.ogen
specifications/core/2.9.0/nwb.ophys
specifications/core/2.9.0/nwb.retinotopy
specifications/hdmf-common
specifications/hdmf-common/1.8.0
specifications/hdmf-common/1.8.0/base
specifications/hdmf-common/1.8.0/namespace
specifications/hdmf-common/1.8.0/sparse
specifications/hdmf-common/1.8.0/table
specifications/hdmf-experimental
specifications/hdmf-experimental/0.5.0
specifications/hdmf-experimental/0.5.0/experimental
specifications/hdmf-experimental/0.5.0/namespace
specifications/hdmf-experimental/0.5.0/resources
stimulus
stimulus/presentation
stimulus/templates
timestamps_reference_time
units
units/id
units/spike_times
units/spike_times_index

Exercise: read the subject’s species from general/subject/species, assign it to a variable named species, and print it.

with h5py.File("synthetic_session.nwb", "r") as f:
    species = f["general/subject/species"][()]
print(species)

b'Mus musculus'

Below is a simplified file tree to illustrate the structure of the different contents of an NWB file. It can be handy to look back on if you are wondering about the relationship between different parts of the NWB file throughout the next sections.

NWBFile  ──────────────────── the whole session
│
├── acquisition ───────────── raw signals, by name
│     └── ElectricalSeries        ← TimeSeries (data + time base + unit)
│
├── processing ────────────── ProcessingModules (derived data)
│     └── "ecephys" ── ProcessingModule
│            └── LFP                  ──────────── a container / data interface
│                 └── ElectricalSeries "lfp"
|     └── "behavior" ── ProcessingModule
│            └── BehavioralTimeSeries ──────────── a container / data interface
│                 └── TimeSeries "running_speed"
│
├── intervals/trials ──────── DynamicTable  (one row per trial)
├── units ──────────────────  DynamicTable  (one row per neuron)
└── general/ ───────────────  metadata
      ├── subject                  ← Subject
      ├── devices/ …
      └── extracellular_ephys/electrodes ← DynamicTable (one row per electrode)

Background

h5py is great for inspection, but for analysis of NWB-files you usually want typed objects (objects with information about what the data is). For example, instead of just a raw array with no information about sampling rates and what kind of data it is, you get an ElectricalSeries object that, together with the data array, contains a table of electrodes and metadata you can access by attribute.

pynwb.NWBHDF5IO opens an NWB file and io.read() reads it and returns an NWBFile object, which is a a typed view of the same HDF5 data you just explored.

Important: keep the file open while you access lazily-loaded data. Once the NWBHDF5IO with-block exits, slicing .data[...] may fail. You can either read inside the with block to automatically close the file after you have read what you wanted to read, or you can explictly call io.close() to close it.

Exercises

Example: open the file with pynwb, read it into an NWBFile, and look at the ElectricalSeries in acquisition by printing it.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    print(nwb.acquisition["ElectricalSeries"])

ElectricalSeries pynwb.ecephys.ElectricalSeries at 0x124534297604848
Fields:
  comments: no comments
  conversion: 1.0
  data: <HDF5 dataset "data": shape (5000, 8), type "<f4">
  description: no description
  electrodes: electrodes <class 'hdmf.common.table.DynamicTableRegion'>
  offset: 0.0
  rate: 1000.0
  resolution: -1.0
  starting_time: 0.0
  starting_time_unit: seconds
  unit: volts

Alternative solution:

io = NWBHDF5IO("synthetic_session.nwb", "r")
nwb = io.read()
print(nwb.acquisition["ElectricalSeries"])
io.close()

ElectricalSeries pynwb.ecephys.ElectricalSeries at 0x124534297799856
Fields:
  comments: no comments
  conversion: 1.0
  data: <HDF5 dataset "data": shape (5000, 8), type "<f4">
  description: no description
  electrodes: electrodes <class 'hdmf.common.table.DynamicTableRegion'>
  offset: 0.0
  rate: 1000.0
  resolution: -1.0
  starting_time: 0.0
  starting_time_unit: seconds
  unit: volts

Exercise: open the file with pynwb, read it into an NWBFile, and print nwb.session_description.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    print(nwb.session_description)

Synthetic ephys session for course exercises

Exercise: open the file and print nwb.subject. Which species was recorded?

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    print(nwb.subject)

subject pynwb.file.Subject at 0x124534296111344
Fields:
  age: P60D
  age__reference: birth
  description: Synthetic subject for the course
  sex: M
  species: Mus musculus
  subject_id: mouse-01

Exercise: assign the raw signal in acquisition/ElectricalSeries to a variable es. Print the sampling rate, starting time, and unit of the raw signal in es.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    es = nwb.acquisition["ElectricalSeries"]
    print(f"rate: {es.rate} Hz; starting_time: {es.starting_time}, unit: {es.unit}")

rate: 1000.0 Hz; starting_time: 0.0, unit: volts

Exercise: Print the shape and data type of the raw signal in acquisition/ElectricalSeries. Note that for shape and dtype you need to access the data attribute of the ElectricalSeries object, i.e. es.data.shape and es.data.dtype.

print(es.data.shape, es.data.dtype)

(5000, 8) float32

Exercise: get the first 200 samples of channel 0 of the raw signal in acquisition/ElectricalSeries and assign to a variable named y. Hint: es.data to get the data. Then, run the cell below to plot the signal.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    es = nwb.acquisition["ElectricalSeries"]
    y = es.data[:200, 0]

fig, ax = plt.subplots()
ax.plot(y)
ax.set_ylabel(f"voltage ({es.unit})")
ax.set_title("Channel 0, first 200 samples")
plt.show()

Exercise: convert the electrode table to a pandas DataFrame, assign it to a variable named edf, and display the first few rows.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    edf = nwb.electrodes.to_dataframe()
edf.head()

Exercise: from the same electrode DataFrame, count how many electrodes are in V1 and how many are in HPC. Hint: See value_counts() in code reference table.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    edf = nwb.electrodes.to_dataframe()
print(edf["location"].value_counts())

location
V1     4
HPC    4
Name: count, dtype: int64

Background

A DynamicTable is a flexible row-and-column data structure used to store, among other things, trial data and units data (activity of neurons). It has certain required columns (e.g. start_time, stop_time for trials; spike_times for units) as well as optional user-defined columns. The easiest way to work with DynamicTables is to call .to_dataframe() and treat it like any other pandas DataFrame.

In our synthetic file:

• The table nwb.trials has 10 trials, each with a start_time, stop_time, and a stim_type of either 'A' or 'B'.
• The table nwb.units has 3 sorted units, each with an array of spike times in seconds.

Exercises

Example: load the trial table as a pandas DataFrame, and display it.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    trials_df = nwb.trials.to_dataframe()
trials_df

Exercise: print the trial table’s column names.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    trial_col_names = nwb.trials.colnames
trial_col_names

('start_time', 'stop_time', 'stim_type')

Exercise: load the units table as a pandas DataFrame, and display it.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    units_df = nwb.units.to_dataframe()
units_df

Example: get the spike times of unit 0 and print the first five timestamps.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    st0 = nwb.units["spike_times"][0]
print(st0[:5])

[0.41953509 0.49117048 0.74386423 1.40682199 1.53924686]

Exercise: get the spike times of unit 2 and print the first five timestamps.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    st2 = nwb.units["spike_times"][2]
print(st2[:5])

[1.17183329 1.49858341 1.64159373 2.0096868  2.10774683]

Exercise: get the spike times of all units and assign them to a variable named sts. Print the first five timestamps of unit 2 to check that you get the same result as in the previous exercise.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    sts = nwb.units["spike_times"][:]

print(sts[0][:5])

[0.41953509 0.49117048 0.74386423 1.40682199 1.53924686]

Demo: Run the cell below to read the spikes from all units and make a spike raster plot.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    spikes = [nwb.units["spike_times"][i] for i in range(len(nwb.units))]

fig, ax = plt.subplots()
ax.eventplot(spikes, lineoffsets=np.arange(len(spikes)), linelengths=0.8)
ax.set_xlabel("time (s)")
ax.set_ylabel("unit")
ax.set_yticks(np.arange(len(spikes)))
ax.set_title("Spike raster")
plt.show()

Background

Data processed from the raw signal, like filtered LFP, behavioural variables, and spike-sorting output, is stored in nwb.processing.

Processed data is grouped by topic into ProcessingModules, one per kind of data (for example one for electrophysiology-derived signals, one for behaviour). Within a module, the data is put in typed containers (such as LFP or BehavioralTimeSeries) that wrap the time series. Reading processed data is therefore a matter of navigating from nwb.processing to the right module, then to the container, and finally to the series inside it.

In our synthetic file, there are two processing modules:

• nwb.processing['ecephys'] contains an LFP container whose electrical_series['lfp'] is the raw signal downsampled to 250 Hz.
• nwb.processing['behavior'] contains a BehavioralTimeSeries container with two behavioural signals sampled at 50 Hz: running_speed (in cm/s) and pupil_diameter (in mm).

Exercises

Example: list the processing modules in the file.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    print(list(nwb.processing.keys()))

['behavior', 'ecephys']

Exercise: list the data interfaces inside the ecephys module.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    print(list(nwb.processing["ecephys"].data_interfaces.keys()))

['LFP']

Exercise: get the LFP time series ElectricalSeries and print its sampling rate. How does it compare to the raw signal’s 1000 Hz?

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    es_lfp = nwb.processing["ecephys"]["LFP"].electrical_series["lfp"]
    print("LFP rate:", es_lfp.rate, "Hz")

LFP rate: 250.0 Hz

Exercise: plot LFP channel 0 over time.

Solution

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    es_lfp = nwb.processing["ecephys"]["LFP"].electrical_series["lfp"]
    y = es_lfp.data[:, 0]
    t = es_lfp.starting_time + np.arange(y.shape[0]) / es_lfp.rate

fig, ax = plt.subplots()
ax.plot(t, y)
ax.set_xlabel("time (s)")
ax.set_ylabel(f"LFP ({es_lfp.unit})")
ax.set_title("LFP channel 0")
plt.show()

Now that you have read the LFP from the ecephys module, turn to the behavioural variables stored in the behavior module.

Example: get the behavior processing module and list the data interfaces (containers) inside it.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    behavior = nwb.processing["behavior"]
    print(list(behavior.data_interfaces.keys()))

['BehavioralTimeSeries']

Exercise: get the BehavioralTimeSeries container from the behavior module and list the behavioural series inside it. Hint: a BehavioralTimeSeries exposes its series through .time_series, just like the LFP container exposes .electrical_series.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    bts = nwb.processing["behavior"]["BehavioralTimeSeries"]
    print(list(bts.time_series.keys()))

['pupil_diameter', 'running_speed']

Exercise: get the running_speed series from the container and print its unit.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    bts = nwb.processing["behavior"]["BehavioralTimeSeries"]
    rs = bts["running_speed"]
    print(f"Running speed unit: {rs.unit}")

Running speed unit: cm/s

Exercise: plot running speed over time. Build the time axis from starting_time and rate the same way you did for the LFP.

Solution

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    rs = nwb.processing["behavior"]["BehavioralTimeSeries"]["running_speed"]
    y = rs.data[:]
    t = rs.starting_time + np.arange(y.shape[0]) / rs.rate

fig, ax = plt.subplots()
ax.plot(t, y)
ax.set_xlabel("time (s)")
ax.set_ylabel(f"running speed ({rs.unit})")
ax.set_title("Running speed")
plt.show()

Exercise: plot the pupil_diameter over time.

Solution

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    pupil = nwb.processing["behavior"]["BehavioralTimeSeries"]["pupil_diameter"]
    y = pupil.data[:]
    t = pupil.starting_time + np.arange(y.shape[0]) / pupil.rate

fig, ax = plt.subplots()
ax.plot(t, y)
ax.set_xlabel("time (s)")
ax.set_ylabel(f"pupil diameter ({pupil.unit})")
ax.set_title("Pupil diameter")
plt.show()

(Extra Challenge) Exercise: compute the mean LFP amplitude per electrode and group the result by brain region (V1 vs HPC) using the electrode table.

with NWBHDF5IO("synthetic_session.nwb", "r") as io:
    nwb = io.read()
    es_lfp = nwb.processing["ecephys"]["LFP"].electrical_series["lfp"]
    data = es_lfp.data[:]
    edf = nwb.electrodes.to_dataframe()

per_channel_mean = data.mean(axis=0)
edf = edf.copy()
edf["lfp_mean"] = per_channel_mean
print(edf.groupby("location")["lfp_mean"].mean())

location
HPC   -5.126826e-07
V1    -1.076278e-06
Name: lfp_mean, dtype: float32