A short introduction to the Visual Behavior Neuropixels data and SDK. This notebook focuses on aligning neural data to visual and optotagging stimuli. For more information about task and behavioral data, check out the 'Aligning behavioral data to task events with the stimulus and trials tables' tutorial.
Also note that this project shares many features with the Visual Coding Neuropixels and Visual Behavior 2-Photon datasets. Users are encouraged to check out the documentation for these projects for additional information and context.
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from allensdk.brain_observatory.behavior.behavior_project_cache.\
behavior_neuropixels_project_cache \
import VisualBehaviorNeuropixelsProjectCache
C:\Users\svc_ccg\AppData\Local\Continuum\anaconda3\lib\site-packages\requests\__init__.py:91: RequestsDependencyWarning: urllib3 (1.26.9) or chardet (3.0.4) doesn't match a supported version! RequestsDependencyWarning)
The VisualBehaviorNeuropixelsProjectCache
is the main entry point to the Visual Behavior Neuropixels dataset. It allows you to download data for individual recording sessions and view cross-session summary information.
# this path determines where downloaded data will be stored
cache_dir = r"C:\Users\svc_ccg\Desktop\Data\vbn_cache"
cache = VisualBehaviorNeuropixelsProjectCache.from_s3_cache(
cache_dir=cache_dir)
# get the metadata tables
units_table = cache.get_unit_table()
channels_table = cache.get_channel_table()
probes_table = cache.get_probe_table()
behavior_sessions_table = cache.get_behavior_session_table()
ecephys_sessions_table = cache.get_ecephys_session_table()
C:\Users\svc_ccg\Documents\GitHub\AllenSDK\allensdk\api\cloud_cache\cloud_cache.py:521: OutdatedManifestWarning: You are loading visual-behavior-neuropixels-2022_project_manifest_v0.1.0.json. A more up to date version of the dataset -- visual-behavior-neuropixels-2022_project_manifest_v0.2.0.json -- exists online. To see the changes between the two versions of the dataset, run VisualBehaviorNeuropixelsProjectCache.compare_manifests('visual-behavior-neuropixels-2022_project_manifest_v0.1.0.json', 'visual-behavior-neuropixels-2022_project_manifest_v0.2.0.json') To load another version of the dataset, run VisualBehaviorNeuropixelsProjectCache.load_manifest('visual-behavior-neuropixels-2022_project_manifest_v0.2.0.json') warnings.warn(msg, OutdatedManifestWarning)
This dataset contains ephys recording sessions from 3 genotypes (C57BL6J, VIP-IRES-CrexAi32 and SST-IRES-CrexAi32). For each mouse, two recordings were made on consecutive days. One of these sessions used the image set that was familiar to the mouse from training. The other session used a new image set for which 6 of the images were novel two were familiar. As an example, let's grab a session from an SST mouse during a novel session.
ecephys_sessions_table = ecephys_sessions_table[0]
sst_novel_sessions = ecephys_sessions_table.loc[(ecephys_sessions_table['genotype'].str.contains('Sst')) &
(ecephys_sessions_table['experience_level']=='Novel')]
sst_novel_sessions.head()
behavior_session_id | date_of_acquisition | equipment_name | session_type | mouse_id | genotype | sex | project_code | age_in_days | unit_count | ... | channel_count | structure_acronyms | image_set | prior_exposures_to_image_set | session_number | experience_level | prior_exposures_to_omissions | file_id | abnormal_activity | abnormal_histology | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
ecephys_session_id | |||||||||||||||||||||
1052530003 | 1052572357 | 2020-09-24 14:34:32.031 | NP.0 | EPHYS_1_images_H_3uL_reward | 524926 | Sst-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt | F | NeuropixelVisualBehavior | 187 | 2005.0 | ... | 2304.0 | ['APN', 'CA1', 'CA3', 'DG-mo', 'DG-po', 'DG-sg... | H | 0.0 | 2 | Novel | 1.0 | 3 | [132] | ['Thalamus'] |
1053941483 | 1053960987 | 2020-10-01 17:03:58.362 | NP.1 | EPHYS_1_images_H_3uL_reward | 527749 | Sst-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt | M | NeuropixelVisualBehavior | 180 | 1543.0 | ... | 2304.0 | ['APN', 'CA1', 'CA3', 'DG-mo', 'DG-po', 'DG-sg... | H | 0.0 | 2 | Novel | 1.0 | 6 | False | NaN |
1046581736 | 1046626452 | 2020-08-27 14:39:26.423 | NP.0 | EPHYS_1_images_H_5uL_reward | 527294 | Sst-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt | M | NeuropixelVisualBehavior | 146 | 2193.0 | ... | 2304.0 | ['APN', 'CA1', 'CA3', 'DG-mo', 'DG-po', 'DG-sg... | H | 0.0 | 2 | Novel | 1.0 | 23 | [111] | ['Hippocampus'] |
1064644573 | 1064666428 | 2020-11-19 15:18:01.372 | NP.1 | EPHYS_1_images_H_3uL_reward | 544456 | Sst-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt | M | NeuropixelVisualBehavior | 120 | 2254.0 | ... | 2304.0 | ['APN', 'CA1', 'CA3', 'DG-mo', 'DG-po', 'DG-sg... | H | 0.0 | 2 | Novel | 1.0 | 27 | False | NaN |
1048189115 | 1048221709 | 2020-09-03 14:16:57.913 | NP.1 | EPHYS_1_images_H_3uL_reward | 509808 | Sst-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt | M | NeuropixelVisualBehavior | 264 | 1925.0 | ... | 2304.0 | ['APN', 'CA1', 'CA3', 'DG-mo', 'DG-po', 'DG-sg... | H | 0.0 | 2 | Novel | 1.0 | 37 | False | NaN |
5 rows × 21 columns
Now let's choose one of these sessions to look at more closely.
session_id = sst_novel_sessions.index.values[0]
session = cache.get_ecephys_session(
ecephys_session_id=session_id)
C:\Users\svc_ccg\AppData\Local\Continuum\anaconda3\lib\site-packages\hdmf\spec\namespace.py:533: UserWarning: Ignoring cached namespace 'hdmf-common' version 1.5.1 because version 1.5.0 is already loaded. % (ns['name'], ns['version'], self.__namespaces.get(ns['name'])['version'])) C:\Users\svc_ccg\AppData\Local\Continuum\anaconda3\lib\site-packages\hdmf\spec\namespace.py:533: UserWarning: Ignoring cached namespace 'core' version 2.4.0 because version 2.3.0 is already loaded. % (ns['name'], ns['version'], self.__namespaces.get(ns['name'])['version'])) C:\Users\svc_ccg\AppData\Local\Continuum\anaconda3\lib\site-packages\hdmf\spec\namespace.py:533: UserWarning: Ignoring cached namespace 'hdmf-experimental' version 0.2.0 because version 0.1.0 is already loaded. % (ns['name'], ns['version'], self.__namespaces.get(ns['name'])['version']))
We can get a high-level summary of this session by accessing its metadata
attribute:
session.metadata
{'equipment_name': 'NP.0', 'sex': 'F', 'age_in_days': 187, 'stimulus_frame_rate': 60.0, 'session_type': 'EPHYS_1_images_H_3uL_reward', 'date_of_acquisition': datetime.datetime(2020, 9, 24, 21, 34, 32, tzinfo=tzutc()), 'reporter_line': 'Ai32(RCL-ChR2(H134R)_EYFP)', 'cre_line': 'Sst-IRES-Cre', 'behavior_session_uuid': None, 'driver_line': ['Sst-IRES-Cre'], 'mouse_id': 524926, 'full_genotype': 'Sst-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt', 'behavior_session_id': 1052572357, 'ecephys_session_id': 1052530003}
Now we can get the unit and channel data for this session. Merging these dataframes will give us CCF coordinates for each unit
units = session.get_units()
channels = session.get_channels()
unit_channels = units.merge(channels, left_on='peak_channel_id', right_index=True)
Let's look at which brain structures were recorded during this session
unit_channels.groupby('structure_acronym')['cluster_id'].count()
structure_acronym APN 214 CA1 185 CA3 16 DG 71 Eth 17 FF 3 HPF 16 HY 10 LGd 38 LP 12 MB 22 MGd 32 MGm 13 MGv 88 MRN 118 NB 39 PIL 52 POST 49 PoT 10 SNr 3 TH 74 VISam 89 VISl 54 VISp 170 VISpm 97 VISrl 161 ZI 33 Name: cluster_id, dtype: int64
Now we'll grab spike times and calculate the change response for 'good' units in V1. Note that how you filter units will depend on your analysis. Consult the unit metrics notebook for more details.
#first let's sort our units by depth
unit_channels = unit_channels.sort_values('probe_vertical_position', ascending=False)
#now we'll filter them
good_unit_filter = ((unit_channels['snr']>1)&
(units['isi_violations']<1)&
(units['firing_rate']>0.1))
good_units = unit_channels.loc[good_unit_filter]
spike_times = session.spike_times
We can get the times when the image changes occurred from the stimulus presentations table. For now, we'll only take the image changes shown during the active behavior block
stimulus_presentations = session.stimulus_presentations
change_times = stimulus_presentations[stimulus_presentations['active']&
stimulus_presentations['is_change']]['start_time'].values
#Convenience function to compute the PSTH
def makePSTH(spikes, startTimes, windowDur, binSize=0.001):
bins = np.arange(0,windowDur+binSize,binSize)
counts = np.zeros(bins.size-1)
for i,start in enumerate(startTimes):
startInd = np.searchsorted(spikes, start)
endInd = np.searchsorted(spikes, start+windowDur)
counts = counts + np.histogram(spikes[startInd:endInd]-start, bins)[0]
counts = counts/startTimes.size
return counts/binSize, bins
We'll include enough time in our plot to see three image responses: the pre-change image response, the change response and the post-change response
#Here's where we loop through the units in our area of interest and compute their PSTHs
area_of_interest = 'VISp'
area_change_responses = []
area_units = good_units[good_units['structure_acronym']==area_of_interest]
time_before_change = 1
duration = 2.5
for iu, unit in area_units.iterrows():
unit_spike_times = spike_times[iu]
unit_change_response, bins = makePSTH(unit_spike_times,
change_times-time_before_change,
duration, binSize=0.01)
area_change_responses.append(unit_change_response)
area_change_responses = np.array(area_change_responses)
#Plot the results
fig, ax = plt.subplots(1,2)
fig.set_size_inches([12,4])
clims = [np.percentile(area_change_responses, p) for p in (0.1,99.9)]
im = ax[0].imshow(area_change_responses, clim=clims)
ax[0].set_title('Active Change Responses for {}'.format(area_of_interest))
ax[0].set_ylabel('Unit number, sorted by depth')
ax[0].set_xlabel('Time from change (s)')
ax[0].set_xticks(np.arange(0, bins.size-1, 20))
_ = ax[0].set_xticklabels(np.round(bins[:-1:20]-time_before_change, 2))
ax[1].plot(bins[:-1]-time_before_change, np.mean(area_change_responses, axis=0), 'k')
ax[1].set_title('{} population active change response (n={})'\
.format(area_of_interest, area_change_responses.shape[0]))
ax[1].set_xlabel('Time from change (s)')
ax[1].set_ylabel('Firing Rate')
Text(0, 0.5, 'Firing Rate')
Now we'll plot receptive fields for these same units. First we need to get stimulus presentation data for the receptive field mapping stimulus (gabors)
rf_stim_table = stimulus_presentations[stimulus_presentations['stimulus_name'].str.contains('gabor')]
xs = np.sort(rf_stim_table.position_x.unique()) #positions of gabor along azimuth
ys = np.sort(rf_stim_table.position_y.unique()) #positions of gabor along altitude
area_rfs = []
for iu, unit in area_units.iterrows():
unit_spike_times = spike_times[iu]
unit_rf = np.zeros([ys.size, xs.size])
for ix, x in enumerate(xs):
for iy, y in enumerate(ys):
stim_times = rf_stim_table[(rf_stim_table.position_x==x)
&(rf_stim_table.position_y==y)]['start_time'].values
unit_response, bins = makePSTH(unit_spike_times,
stim_times+0.01,
0.2, binSize=0.001)
unit_rf[iy, ix] = unit_response.mean()
area_rfs.append(unit_rf)
fig, axes = plt.subplots(int(len(area_rfs)/10)+1, 10)
for irf, rf in enumerate(area_rfs):
ax_row = int(irf/10)
ax_col = irf%10
axes[ax_row][ax_col].imshow(rf, origin='lower')
for ax in axes.flat:
ax.axis('off')
Since this is an SST mouse, we should see putative SST+ interneurons that are activated during our optotagging protocol. Let's load the optotagging stimulus table and plot PSTHs triggered on the laser onset. For more examples, you can check out the Visual Coding Neuropixels Optagging notebook here
opto_table = session.optotagging_table
opto_table.head()
start_time | condition | level | stop_time | stimulus_name | duration | |
---|---|---|---|---|---|---|
id | ||||||
0 | 8803.77363 | half-period of a cosine wave | 0.78 | 8804.77363 | raised_cosine | 1.00 |
1 | 8805.63517 | a single square pulse | 1.00 | 8805.64517 | pulse | 0.01 |
2 | 8807.35557 | a single square pulse | 1.70 | 8807.36557 | pulse | 0.01 |
3 | 8809.04593 | a single square pulse | 1.00 | 8809.05593 | pulse | 0.01 |
4 | 8811.07608 | half-period of a cosine wave | 0.78 | 8812.07608 | raised_cosine | 1.00 |
If you check out this table, you'll see that we use 2 different laser waveforms: a half-period cosine wave that's 1 second long and a short square pulse that's 10 ms long. We drive each at three light levels, giving us 6 total conditions. Now let's plot how cortical neurons respond to the short pulse at high power.
duration = opto_table.duration.min()
level = opto_table.level.max()
cortical_units = good_units[good_units['structure_acronym'].str.contains('VIS')]
opto_times = opto_table.loc[(opto_table['duration']==duration)&
(opto_table['level']==level)]['start_time'].values
opto_response = []
unit_id = []
for iu, unit in cortical_units.iterrows():
unit_spike_times = spike_times[iu]
unit_response, bins = makePSTH(unit_spike_times,
opto_times-0.05, 0.1,
binSize=0.0005)
opto_response.append(unit_response)
unit_id.append(iu)
opto_response = np.array(opto_response)
fig, ax = plt.subplots()
fig.set_size_inches((5,10))
fig.suptitle('Optotagging: ' + str(session.metadata['ecephys_session_id'])
+ ' ' + session.metadata['full_genotype'])
im = ax.imshow(opto_response[:, 80:140], extent=[-10, 20, 0, len(unit_id)],
origin='lower', aspect='auto')
min_clim_val = 0
max_clim_val = 300
im.set_clim([min_clim_val, max_clim_val])
[ax.axvline(bound, linestyle=':', color='white', linewidth=1.0)\
for bound in [0.25, 10.25]]
ax.set_xlabel('Time from laser onset (ms)')
ax.set_ylabel('Unit number')
Text(0, 0.5, 'Unit number')
Here we can see that most units don't respond to the short laser pulse. But there is a population of ~20 units that do show elevated firing rates. Note that the activity occurring at the onset and offset of the laser (along the dotted lines) is likely artifactual and should be excluded from analysis.