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⚠️

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 .

c:\Users\carter.peene\AppData\Local\Programs\Python\Python39\lib\site-packages\numpy\_distributor_init.py:30: UserWarning: loaded more than 1 DLL from .libs:
c:\Users\carter.peene\AppData\Local\Programs\Python\Python39\lib\site-packages\numpy\.libs\libopenblas.FB5AE2TYXYH2IJRDKGDGQ3XBKLKTF43H.gfortran-win_amd64.dll
c:\Users\carter.peene\AppData\Local\Programs\Python\Python39\lib\site-packages\numpy\.libs\libopenblas.XWYDX2IKJW2NMTWSFYNGFUWKQU3LYTCZ.gfortran-win_amd64.dll
  warnings.warn("loaded more than 1 DLL from .libs:"

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

from scipy import interpolate

%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 fsspec for more details on this.

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

# 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()

A newer version (0.59.0) of dandi/dandi-cli is available. You are using 0.55.1

File already exists
Opening file

c:\Users\carter.peene\AppData\Local\Programs\Python\Python39\lib\site-packages\hdmf\spec\namespace.py:531: UserWarning: Ignoring cached namespace 'hdmf-common' version 1.5.1 because version 1.8.0 is already loaded.
  warn("Ignoring cached namespace '%s' version %s because version %s is already loaded."
c:\Users\carter.peene\AppData\Local\Programs\Python\Python39\lib\site-packages\hdmf\spec\namespace.py:531: UserWarning: Ignoring cached namespace 'core' version 2.6.0-alpha because version 2.5.0 is already loaded.
  warn("Ignoring cached namespace '%s' version %s because version %s is already loaded."
c:\Users\carter.peene\AppData\Local\Programs\Python\Python39\lib\site-packages\hdmf\spec\namespace.py:531: UserWarning: Ignoring cached namespace 'hdmf-experimental' version 0.2.0 because version 0.5.0 is already loaded.
  warn("Ignoring cached namespace '%s' version %s because version %s is already loaded."

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.

dff = nwb.processing["ophys"]["dff"]
dff_trace = dff.roi_response_series["traces"].data
dff_timestamps = dff.roi_response_series["traces"].timestamps

# accessing the above data may look different for older nwb files, like the following
# dff_trace = dff.roi_response_series["RoiResponseSeries"].data
# dff_timestamps = dff.roi_response_series["RoiResponseSeries"].timestamps

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

(115870, 178)
(115870,)

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

['trials']

stim_table = nwb.intervals["trials"]
print(stim_table.colnames)
stim_table[:10]

('start_time', 'stop_time', 'stimulus_type', 'image_id', 'repeat', 'stimulus_key', 'fraction_occlusion', 'duration', 'sweep')

	start_time	stop_time	stimulus_type	image_id	repeat	stimulus_key	fraction_occlusion	duration	sweep
id
0	102.580917	102.814461	randomized_control_pre	103.0	NaN	ophys_pilot_randomized_control_A	NaN	0.25	0
1	102.831096	103.064622	randomized_control_pre	68.0	NaN	ophys_pilot_randomized_control_A	NaN	0.25	1
2	103.081333	103.314856	randomized_control_pre	26.0	NaN	ophys_pilot_randomized_control_A	NaN	0.25	2
3	103.331516	103.565072	randomized_control_pre	71.0	NaN	ophys_pilot_randomized_control_A	NaN	0.25	3
4	103.581687	103.815235	randomized_control_pre	110.0	NaN	ophys_pilot_randomized_control_A	NaN	0.25	4
5	103.831911	104.065457	randomized_control_pre	6.0	NaN	ophys_pilot_randomized_control_A	NaN	0.25	5
6	104.082153	104.315617	randomized_control_pre	112.0	NaN	ophys_pilot_randomized_control_A	NaN	0.25	6
7	104.332359	104.565821	randomized_control_pre	17.0	NaN	ophys_pilot_randomized_control_A	NaN	0.25	7
8	104.582515	104.816081	randomized_control_pre	13.0	NaN	ophys_pilot_randomized_control_A	NaN	0.25	8
9	104.832731	105.066273	randomized_control_pre	111.0	NaN	ophys_pilot_randomized_control_A	NaN	0.25	9

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"

### 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)

6
('trials', 'randomized_control_pre', 102.58091661959725, 207.73271267031623)
('trials', 'oddball', 267.7980219847227, 2269.8206912029705)
('trials', 'transition_control', 2329.885995836317, 2690.2113244542556)
('trials', 'occlusion', 2750.777092659118, 3350.829824249638)
('trials', 'natural_movie_one', 3410.895166331465, 3711.1216876550166)
('trials', 'randomized_control_post', 3771.2371179895413, 3876.3055732709176)

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()

174.87656701939662

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)

13924 22968

### 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)

Text(0.5, 1.0, 'ROIs DF/F Over Time')