Visualizing 2P Responses to Stimulus¶

Some of this content is adapted from the Allen SDK Documentation.

In 2-Photon Image NWB Files, the data stored measure fluorescence. From the raw 2P images, the The Allen Institute uses the Suite2P segmentation algorithm to identify regions of interest (ROIs). For more information on ROIs go to the Visualizing Raw 2-Photon Images notebook. For each ROI, two measures of fluorescence are stored; corrected fluorescence and DF/F. Corrected fluorescence (hereafter referred to as ‘fluorescence’) is the brightness of the region with unknown units. DF/F is a normalized metric that takes into account baseline fluorescence value in order to provide a measure that can be compared from cell to cell. These traces are stored in the Processing section. This notebook uses these fluorescence traces and the stimulus information from the tables in the Intervals section of the file to visualize neuron responses to stimulus.

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_download_open
except:
    !git clone https://github.com/AllenInstitute/openscope_databook.git
    %cd openscope_databook
    %pip install -e .

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np

%matplotlib inline

Downloading Ophys File¶

Change the values below to download the file you’re interested in. Set dandiset_id and dandi_filepath to correspond to the dandiset id and filepath of the file you want. If you’re accessing an embargoed dataset, set dandi_api_key to your DANDI API key. If you want to stream a file instead of downloading it, use dandi_stream_open instead. Checkout Streaming an NWB File with remfile for more details on this.

dandiset_id = "000488"
dandi_filepath = "sub-416366/sub-416366_ses-774663713_behavior+image+ophys.nwb"
download_loc = "."
dandi_api_key = None

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_download_open(dandiset_id, dandi_filepath, download_loc, dandi_api_key=dandi_api_key)
nwb = io.read()

File already exists
Opening file

Extracting 2P Data and Stimulus Data¶

Below, the fluorescence traces and timestamps are read from the file’s Processing section. Note that the exact format to access these traces can vary between newer and older NWB files, so some adjustments may be necessary. Additionally, the stimulus data is also read from the NWB file’s Intervals section. Stimulus information is stored as a series of tables depending on the type of stimulus shown in the session. One such table is displayed below, as well as the names of all the tables.

# for newer OpenScope multiplane NWBs, each plane is its own processing module;
# set plane_name to the plane you want (e.g. "VISp_0", "VISp_1", "VISl_4", etc.)
plane_name = "VISp_0"
dff = nwb.processing[plane_name]["dff_timeseries"]
dff_trace = dff.roi_response_series["dff_timeseries"].data
dff_timestamps = dff.roi_response_series["dff_timeseries"].timestamps

# for older OpenScope NWBs, dff is stored under "ophys"
# dff = nwb.processing["ophys"]["dff"]
# dff_trace = dff.roi_response_series["traces"].data
# dff_timestamps = dff.roi_response_series["traces"].timestamps

# for even older NWB files the series key may differ:
# dff_trace = dff.roi_response_series["RoiResponseSeries"].data
# dff_timestamps = dff.roi_response_series["RoiResponseSeries"].timestamps

print(dff_trace.shape)
print(dff_timestamps.shape)

(46303, 345)
(46303,)

stimulus_names = list(nwb.intervals.keys())
print(stimulus_names)

['Control block 1_presentations', 'Control block 2_presentations', 'Control block 3_presentations', 'Control block 4_presentations', 'RF mapping_presentations', 'Sensory-motor mismatch block_presentations', 'Zebra_presentations', 'spontaneous_presentations']

stim_name = "trials"
stim_name = "Sensory-motor mismatch block_presentations"
stim_table = nwb.intervals[stim_name]
print(stim_table.colnames)
stim_table[:10]

('start_time', 'stop_time', 'stim_name', 'stim_type', 'stim_block', 'Orientation', 'SpatialFrequency', 'TemporalFrequency', 'contrast', 'phase', 'DiameterX', 'DiameterY', 'X', 'Y', 'Duration', 'Delay', 'BlockNumber', 'BlockLabel', 'TrialNumber', 'SequenceNumber', 'TrialInSequence', 'TrialType', 'BlockType', 'stim_index', 'timeseries')

Getting Stimulus Epochs¶

Here, epochs are extracted from the stimulus tables. In this case, an epoch is a continuous period of time during a session where a particular type of stimulus is shown. The output here is a list of epochs, where an epoch is a tuple of four values; the stimulus name, the stimulus block, the starting time and the ending time. Since stimulus information can vary significantly between experiments and NWB files, you may need to tailor the code below to extract epochs for the file you’re interested in. For normally, formatted stimulus tables, the code below should work, provided you set epoch_criterion to be the name of a valid column in the stimulus tables above. In this example, we define an epoch as ending when the “stimulus_type” value changes for a row.

epoch_criterion = "stimulus_type" # also try "stimulus_key"
# epoch_criterion = "stim_block"

### extract epoch times from stim table where stimulus rows have a different 'block' than following row
### returns list of epochs, where an epoch is of the form (stimulus name, stimulus block, start time, stop time)
def extract_epochs(stim_name, stim_table, epochs):    
    # set first epoch bounds to first start and stop times in stim table
    epoch_start = stim_table.start_time[0]
    epoch_stop = stim_table.stop_time[0]

    # for each row, try to extend current epoch stop_time
    for i in range(len(stim_table)):
        this_block = stim_table[epoch_criterion][i]
        # if end of table, end the current epoch
        if i+1 >= len(stim_table):
            epochs.append((stim_name, this_block, epoch_start, epoch_stop))
            break
            
        next_block = stim_table[epoch_criterion][i+1]
        # if next row is the same stim block, push back epoch_stop time
        if next_block == this_block:
            epoch_stop = stim_table.stop_time[i+1]
        # otherwise, end the current epoch, start new epoch
        else:
            epochs.append((stim_name, this_block, epoch_start, epoch_stop))
            epoch_start = stim_table.start_time[i+1]
            epoch_stop = stim_table.stop_time[i+1]
    
    return epochs

### extract epochs from all valid stimulus tables
epochs = []
for stim_name in stimulus_names:
    stim_table = nwb.intervals[stim_name]
    try:
        epochs = extract_epochs(stim_name, stim_table, epochs)
    except:
        continue

# epochs take the form (stimulus name, stimulus block, start time, stop time)
print(len(epochs))
# sorts epochs based on epoch start time 
epochs.sort(key=lambda x: x[2])
for epoch in epochs:
    print(epoch)

Fluorescence Activity Throughout Epochs¶

Below is a view of the fluorescence activity of ROIs throughout a session where epochs are shown as colored sections. First, set trace_to_display to choose “fluorescence” or “dff”. You can view the activity of individual ROIs, all ROIs, or an average of all ROIs. Set roi_num to be the ID of the ROI you want to view, or to None to view all ROIs. Alternatively, set display_average to True to see the average of all ROIs. Set time_start to the starting bound in seconds of the session, you’d like to see, and time_end to the ending bound. You may want to use the output above to inform this choice. As mentioned above, if your file’s stimulus information differs significantly, the plotting code below may need to be modified to appropriately display the epochs.

roi_num = 13 # chosen from dff_trace.shape[1]
display_average_trace = False
time_start = 0
time_end = 4000

### choose data and plot properties based on user input

avg_dff_trace = np.average(dff_trace, axis=1)

if display_average_trace:
    plot_title = "Average DF/F Throughout Epochs"
    ylabel = "Average DF/F (%)"
    display_trace = avg_dff_trace*100
else:
    plot_title = "DF/F Throughout Epochs"
    ylabel = "DF/F (%)"
    if roi_num:
        display_trace = dff_trace[:,roi_num]*100
    else:
        display_trace = dff_trace*100
        # ^^^ * 100 to yield percentage for DF/F

### make plot of chosen fluorescence trace over time with colored epoch sections

fig, ax = plt.subplots(figsize=(15,5))

# filter epochs which aren't at least partially in the time window
bounded_epochs = {epoch for epoch in epochs if epoch[2] < time_end and epoch[3] > time_start}

# assign unique color to each stimulus name
stim_names = list({(epoch[0], epoch[1])for epoch in bounded_epochs})
colors = plt.cm.rainbow(np.linspace(0,1,len(stim_names)))
stim_color_map = {stim_names[i]:colors[i] for i in range(len(stim_names))}

epoch_key = {}
y_hi = np.amax(display_trace) # change these to manually set height of the plot
y_lo = np.amin(display_trace)
# draw colored rectangles for each epoch
for epoch in bounded_epochs:
    stim_name, stim_block, epoch_start, epoch_end = epoch
    color = stim_color_map[(stim_name, stim_block)]
    rec = ax.add_patch(mpl.patches.Rectangle((epoch_start, y_lo), epoch_end-epoch_start, y_hi, alpha=0.3, facecolor=color))
    epoch_key[(stim_name, stim_block)] = rec

ax.set_xlim(time_start, time_end)
ax.set_ylim(y_lo, y_hi)
ax.set_xlabel("Time (s)")
ax.set_ylabel(ylabel)
ax.set_title(plot_title)

fig.legend(epoch_key.values(), epoch_key.keys(), loc="lower right", bbox_to_anchor=(1.2, 0.25))
ax.plot(dff_timestamps[:], display_trace)

print(np.amax(display_trace))
plt.tight_layout()
plt.show()

418.55164

Region-wise Fluorescence Over Time¶

Here a 2D image is made showing ROI fluorescence values over time. Without organizing the ROIs, there are limits to how useful it can be, but it can serve to show responses to stimulus times. Set trace_to_display to “fluorescence” or “dff” to select the trace you want to see. Set interval_start and interval_end to the bounds of time, in seconds, you’d like to see displayed. These bounds get translated into the indices of the trace that match most closely to the time that was provided. If the bounds are outside the measurement period, an error will be thrown. You may also set start_roi and end_roi to narrow the number of ROIs that will be displayed.

interval_start = 500
interval_end = 800
start_roi = 0
end_roi = dff_trace.shape[1] # show all ROIs

# translate interval time bounds into approximate data indices
interval_start_idx, interval_end_idx = np.searchsorted(dff_timestamps, (interval_start, interval_end))

if interval_start_idx >= interval_end_idx:
    raise IndexError(f"Interval end is not greater than interval start; ({interval_start_idx,interval_end_idx}); Interval has non-positive length")
    
print(interval_start_idx, interval_end_idx)

4944 8155

### display array of ROI fluorescence traces as 2D image with color
fig, ax = plt.subplots(figsize=(10,10))

vmin, vmax = 0, 5 # adjust these to change contrast of fluorescence values in the plot
img = ax.imshow(np.transpose(dff_trace[interval_start_idx:interval_end_idx]), extent=[interval_start,interval_end,start_roi,end_roi], aspect="auto", vmax=vmin, vmin=vmax)
cbar = plt.colorbar(img, shrink=0.5)
cbar.set_label("DF/F (%)")

ax.yaxis.set_major_locator(plt.NullLocator())
ax.set_ylabel("ROIs", fontsize=16)
ax.set_xlabel("Time (s)", fontsize=16)
ax.set_title("ROIs DF/F Over Time", fontsize=20)