Visualizing Behavioral Data - OpenScope Databook

A number of different behavioral measurements during experiments and stored in NWB Files. These include various kinds of eye tracking. Eye tracking provides a proxy into the mouse visual attention and its overall brain state. Also included is data of the running wheel that the mouse is on during recording. Oftentimes other measurements are taken like mouse lick times, but this is not included here. This notebook just has basic code to take these behavioral data from an ecephys NWB file and plot them.

Environment Setup¶

⚠️Note: If running on a new environment, run this cell once and then restart the kernel⚠️

import warnings
warnings.filterwarnings('ignore')

try:
    from databook_utils.dandi_utils import dandi_stream_open
except:
    !git clone https://github.com/AllenInstitute/openscope_databook.git
    %cd openscope_databook
    %pip install -e .

import os

import matplotlib.pyplot as plt
import numpy as np

%matplotlib inline

Streaming NWB File¶

Streaming a file from DANDI requires information about the file of interest. The current information below is for data that is private to the Allen Institute. Set dandiset_id to be the ID of the dandiset you want, and set dandi_filepath to be the path of the file within the dandiset. The filepath can be found if you press on the i icon of a file and copy the path field that shows up in the resulting JSON. If you are accessing embargoed data, you will need to set dandi_api_key to your DANDI API key.

dandiset_id = "000336"
dandi_filepath = "sub-621602/sub-621602_ses-1193555033-acq-1193675745_ophys.nwb"

# use for OpenScope Predictive Processing Community Project NWBs
# dandiset_id = "001768"
# dandi_filepath = "sub-839909/sub-839909_ses-multiplane-ophys-839909-2026-02-20-12-53-27_ophys.nwb"

# This can sometimes take a while depending on the size of the file
io = dandi_stream_open(dandiset_id, dandi_filepath)
nwb = io.read()

Extracting Eye Tracking Data¶

Our datasets include eye data with eye tracking, corneal reflection tracking, and pupil tracking. These are different visual components identified in footage of the mouse during the experiment. Any of these could be usable for the following analyses. Each eye tracking type includes measurements of the angle of the eye and the height, width, and area (in terms of pixels on the camera), and the coordinates of the center of the component.

The cell below auto-detects which NWB format is present and normalizes the data into plain numpy arrays used by all cells below:

Older NWBs (standard ecephys/ophys): eye tracking in nwb.acquisition["EyeTracking"] as SpatialSeries. Area is available.
Newer NWBs (multiplane ophys): eye tracking in nwb.processing["eye_tracking"] as DynamicTables. Area may not be available.

if "EyeTracking" in nwb.acquisition:
    # Older NWBs: eye tracking in nwb.acquisition["EyeTracking"] as SpatialSeries
    et = nwb.acquisition["EyeTracking"]
    eye_tracking = et.eye_tracking
    timestamps  = np.array(eye_tracking.timestamps)
    blink_times = np.array(et.likely_blink.data)
    eye_width   = np.array(eye_tracking.width)
    eye_height  = np.array(eye_tracking.height)
    eye_angle   = np.array(eye_tracking.angle)
    eye_xs      = np.array([pt[0] for pt in eye_tracking.data])
    eye_ys      = np.array([pt[1] for pt in eye_tracking.data])
    has_eye_area = hasattr(eye_tracking, "area") and eye_tracking.area is not None
    eye_area     = np.array(eye_tracking.area) if has_eye_area else None
else:
    # Newer NWBs: eye tracking in nwb.processing["eye_tracking"] as DynamicTables
    et_table    = nwb.processing["eye_tracking"]["ellipse"]
    timestamps  = np.array(et_table["timestamps"][:])
    blink_times = np.array(nwb.processing["eye_tracking"]["likely_blink_times"].data)
    eye_width   = np.array(et_table["width"][:])
    eye_height  = np.array(et_table["height"][:])
    eye_angle   = np.array(et_table["angle"][:])
    eye_xs      = np.array(et_table["data_x"][:])
    eye_ys      = np.array(et_table["data_y"][:])
    has_eye_area = "area" in et_table.colnames
    eye_area     = np.array(et_table["area"][:]) if has_eye_area else None

print(timestamps.shape)

(241673,)

Selecting a Period¶

The data can be large or messy. In order to visualize the data more cleanly and efficiently, you can just select a period of time within the data to plot. To do this, specify the start_time and end_time you’d like in terms of seconds. Below, the first and last timestamps from the data are printed to inform this choice. Once you input your start time and end time in seconds, those times will be translated into indices so the program can identify what slice of data to select.

print(timestamps[0])
print(timestamps[-2]) # not [-1] because that's NaN

0.20801
4027.98269

start_time = 1000
end_time = 1200

### translate times to data indices using timestamps

start_idx, end_idx = None, None
for i, ts in enumerate(timestamps):
    if not start_idx and ts >= start_time:
        start_idx = i
    if start_idx and ts >= end_time:
        end_idx = i
        break

if start_idx == None or end_idx == None:
    raise ValueError("Time bounds not found within eye tracking timestamps")

# make time axis
time_axis = np.arange(start_idx, end_idx)

Visualizing the Data¶

Below, several types of measurements are visualized from the eye tracking data that was selected.

Blink Times¶

fig, ax = plt.subplots()
ax.plot(time_axis, blink_times[start_idx:end_idx], linewidth=0.2)
ax.set_xlabel("time (s)")
ax.set_title("Blink Times")

Area¶

Below, eye height and width are plotted together, and area, the product of height and width, is also plotted. The units of height and width are defined in terms of pixels on the eye tracking camera.

fig, ax1 = plt.subplots()
ax1.set_xlabel('time (s)')
ax1.plot(time_axis, eye_width[start_idx:end_idx], color='b')
ax1.set_ylabel('width (pixels)', color='b')
ax2 = ax1.twinx()
ax2.plot(time_axis, eye_height[start_idx:end_idx], color='r')
ax2.set_ylabel('height (pixels)', color='r')
ax1.set_title("Eye Height and Width Over Time")

if has_eye_area:
    fig, ax1 = plt.subplots()
    ax1.set_xlabel('time (s)')
    ax1.set_ylabel('area (pixels)')
    ax1.set_title("Area Over Time")
    ax1.plot(time_axis, eye_area[start_idx:end_idx])
else:
    print("Area data not available for this NWB file.")

Angle¶

fig, ax = plt.subplots()
ax.set_xlabel('time (s)')
ax.set_ylabel('angle (degrees)')
ax.set_title("Eye Angle Over Time")
ax.plot(time_axis, eye_angle[start_idx:end_idx])

Eye Trace¶

With marker color representing time, the x and y coordinates of the eye’s view are traced below.

fig, ax = plt.subplots()
colors = plt.cm.viridis(np.linspace(0, 1, end_idx-start_idx))
ax.plot(eye_xs[start_idx:end_idx], eye_ys[start_idx:end_idx], zorder=0, linewidth=0.25)
ax.scatter(eye_xs[start_idx:end_idx], eye_ys[start_idx:end_idx], s=5, c=colors, zorder=1)

# change these to set the plot limits
# ax.set_xlim(310,360)
# ax.set_ylim(270,320)

ax.set_xlabel("x pixel")
ax.set_ylabel("y pixel")
ax.set_title("Eye Trace Through Time")
plt.show()

Running Data¶

Apart from the Eye Tracking Data, the mouse’s running data is tracked in the "running" field of the processing section of the NWB File. Speed, in cm/s is tracked in the "speed" field. Another field, "dx" is also recorded which represents the number of degrees the wheel has turned between timestamp. These are plotted below. It is important to note that, at the moment, running data is stored in NWB differently in our Ecephys files and our Ophys files. Change the code in the cell below if you’re accessing a different type of file.

running = nwb.processing["running"]
running_speed  = running["speed"] if "speed" in running.keys() else running["running_speed"]
wheel_rotation = running["dx"] if "dx" in running.keys() else running["running_wheel_rotation"]

Running Speed¶

Below are the running speed over time and wheel rotation over time. They should probably look very similar. Additionally, we show the distribution of the timestamp deltas. If the data was processed properly, there shouldn’t be many deviations. However, sometimes, individual frames might be dropped during processing from the experimental rigs which could show some diffs that are double what the expected value is.

speed_data = np.array(running_speed.data)
speed_timestamps = np.array(running_speed.timestamps)
print(speed_data.shape)
print(speed_timestamps.shape)

fig, ax = plt.subplots()
ax.set_xlabel("time (s)")
ax.set_ylabel("speed (cm/s)")
ax.set_title("Running Speed Over Time")
ax.plot(speed_timestamps, speed_data)

(236232,)
(236232,)

fig, ax = plt.subplots()
ax.hist(np.diff(speed_timestamps), 100)
ax.set_xlabel("Delta (s)")
ax.set_ylabel("# Deltas")
ax.set_title("Running Timestamp Deltas Histogram")
plt.show()

rotation_data = np.array(wheel_rotation.data[1:]) # cut out first value since it is high outlier
rotation_timestamps = np.array(wheel_rotation.timestamps[1:])
print(rotation_data.shape)
print(rotation_timestamps.shape)

fig, ax = plt.subplots()
ax.set_xlabel("time (s)")
ax.set_ylabel("rotation (degrees)")
ax.set_title("Wheel Rotation Over Time")
ax.plot(rotation_timestamps, rotation_data)

(236231,)
(236231,)