Brain Observatory¶
This notebook documents some classes and functions in the AllenSDK that help manipulate files and data structures in the Allen Brain Observatory. The main entry point in the BrainObservatoryCache class. This class is responsible for downloading any requested data or metadata on request and storing it in well known locations.
Download this file in .ipynb format here.
Experiment Containers¶
The experiment container describes a set of experiments performed with the same targeted area, imaging depth, and Cre line. The BrainObservatoryCache has a number of functions for figuring out what experiment containers are available at the moment.
from allensdk.core.brain_observatory_cache import BrainObservatoryCache
import pprint
# This class uses a 'manifest' to keep track of downloaded data and metadata.
# All downloaded files will be stored relative to the directory holding the manifest
# file. If 'manifest_file' is a relative path (as it is below), it will be
# saved relative to your working directory. It can also be an absolute path.
boc = BrainObservatoryCache(manifest_file='boc/manifest.json')
# Download a list of all targeted areas
targeted_structures = boc.get_all_targeted_structures()
print("all targeted structures: " + str(targeted_structures))
# Download experiment containers for VISp experiments
visp_ecs = boc.get_experiment_containers(targeted_structures=['VISp'])
print("all VISp experiment containers: %d" % len(visp_ecs))
# Download a list of all imaging depths
depths = boc.get_all_imaging_depths()
print("all imaging depths: " + str(depths))
# Download a list of all stimuli
stims = boc.get_all_stimuli()
print("all stimuli:\n")
pprint.pprint(stims)
# Download a list of all cre driver lines
cre_lines = boc.get_all_cre_lines()
print("all cre lines:\n")
pprint.pprint(cre_lines)
# Download experiment containers for Cux2 experiments
cux2_ecs = boc.get_experiment_containers(cre_lines=['Cux2-CreERT2'])
print("Cux2 experiments: %d\n" % len(cux2_ecs))
print("Example experiment container record:")
pprint.pprint(cux2_ecs[0])
Download Experiments for a Container¶
An experiment container is a group of experiments. Each experiment has a different stimulus protocol. For example, one experiment protocol contains the static gratings stimulus and another has the natural scenes stimulus. The BrainObservatoryCache helps you find out which experiment associated with a container has the stimuli you are interested in. First, let's see what experiments are available for a single container.
# Find all of the experiments for an experiment container
cux2_ec_id = cux2_ecs[0]['id']
exps = boc.get_ophys_experiments(experiment_container_ids=[cux2_ec_id])
print("Experiments for experiment_container_id %d: %d\n" % (cux2_ec_id, len(exps)))
pprint.pprint(exps)
The session_type field indicates which experimental protocol was used. If you just want to find which experiment contains the static gratings stimulus, you can do the following:
import allensdk.brain_observatory.stimulus_info as stim_info
# pick one of the cux2 experiment containers
cux2_ec_id = cux2_ecs[-1]['id']
# Find the experiment with the static static gratings stimulus
exp = boc.get_ophys_experiments(experiment_container_ids=[cux2_ec_id],
stimuli=[stim_info.STATIC_GRATINGS])[0]
print("Experiment with static gratings:")
pprint.pprint(exp)
Now we can download the NWB file for this experiment.
exp = boc.get_ophys_experiment_data(exp['id'])
# print out the metadata available in the NWB file
pprint.pprint(exp.get_metadata())
Find Cells of Interest¶
Another way to look for data is to search for cells with interesting tuning properties. We have pre-computed a set of cell metrics. Take a look at this page to see the full list.
You can use these properties to filter the full set of cells down to the set your are interested in.
import pandas as pd
# Download cells for a set of experiments and convert to DataFrame
cells = boc.get_cell_specimens()
cells = pd.DataFrame.from_records(cells)
print("total cells: %d" % len(cells))
# find direction selective cells in VISp
visp_ec_ids = [ ec['id'] for ec in visp_ecs ]
visp_cells = cells[cells['experiment_container_id'].isin(visp_ec_ids)]
print("VISp cells: %d" % len(visp_cells))
# significant response to drifting gratings stimulus
sig_cells = visp_cells[visp_cells['p_dg'] < 0.05]
print("cells with sig. response to drifting gratings: %d" % len(sig_cells))
# direction selective cells
dsi_cells = sig_cells[(sig_cells['g_dsi_dg'] > 0.9)]
print("direction-selective cells: %d" % len(dsi_cells))
Find Experiments for Cells¶
Once you have found a set of cells you would like to investigate, you can use the BrainObservatoryCache to find the experiments for those cells that contain the relevant stimulus.
import allensdk.brain_observatory.stimulus_info as stim_info
# find experiment containers for those cells
dsi_ec_ids = dsi_cells['experiment_container_id'].unique()
print("total dsi experiment containers: %d" % len(dsi_ec_ids))
# Download the ophys experiments containing the drifting gratings stimulus for VISp experiment containers
dsi_exps = boc.get_ophys_experiments(experiment_container_ids=dsi_ec_ids, stimuli=[stim_info.DRIFTING_GRATINGS])
print("VISp drifting gratings ophys experiments: %d" % len(dsi_exps))
print("Example ophys experiment:")
pprint.pprint(dsi_exps[0])
Download Experiment Data for a Cell¶
Once you have some experiments, you can download the NWB files that contain the fluorescence traces for segmented cells in those experiments.
# pick a direction-selective cell and find its NWB file
dsi_cell = dsi_cells.iloc[0]
# figure out which ophys experiment has the drifting gratings stimulus for that cell
cell_exp = boc.get_ophys_experiments(cell_specimen_ids=[dsi_cell['cell_specimen_id']],
stimuli=[stim_info.DRIFTING_GRATINGS])[0]
data_set = boc.get_ophys_experiment_data(cell_exp['id'])
print("Metadata from NWB file:")
pprint.pprint(data_set.get_metadata())
print("stimuli available in this file:")
print(data_set.list_stimuli())
Fluorescence Traces¶
Now that we have a data set, we can plot the traces for the cell we care about.
dsi_cell_id = dsi_cell['cell_specimen_id']
time, raw_traces = data_set.get_fluorescence_traces(cell_specimen_ids=[dsi_cell_id])
_, demixed_traces = data_set.get_demixed_traces(cell_specimen_ids=[dsi_cell_id])
_, neuropil_traces = data_set.get_neuropil_traces(cell_specimen_ids=[dsi_cell_id])
_, corrected_traces = data_set.get_corrected_fluorescence_traces(cell_specimen_ids=[dsi_cell_id])
_, dff_traces = data_set.get_dff_traces(cell_specimen_ids=[dsi_cell_id])
from matplotlib import pyplot as plt
%matplotlib inline
# plot raw and corrected ROI trace
plt.figure(figsize=(14,4))
plt.title("Raw Fluorescence Trace")
plt.plot(time, raw_traces[0])
plt.show()
plt.figure(figsize=(14,4))
plt.title("Demixed Fluorescence Trace")
plt.plot(time, demixed_traces[0])
plt.show()
plt.figure(figsize=(14,4))
plt.title("Neuropil-corrected Fluorescence Trace")
plt.plot(time, corrected_traces[0])
plt.show()
plt.figure(figsize=(14,4))
plt.title("dF/F Trace")
# warning: dF/F can occasionally be one element longer or shorter
# than the time stamps for the original traces.
plt.plot(time[:len(dff_traces[0])], dff_traces[0])
plt.show()
ROI Masks¶
If you want to take a look at the cell visually, you can open the NWB file and extract a pixel mask. You can also pull out the maximum intensity projection of the movie for context.
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
data_set = boc.get_ophys_experiment_data(510221121)
# get the specimen IDs for a few cells
cids = data_set.get_cell_specimen_ids()[:15:5]
# get masks for specific cells
roi_mask_list = data_set.get_roi_mask(cell_specimen_ids=cids)
# plot each mask
f, axes = plt.subplots(1, len(cids)+2, figsize=(15, 3))
for ax, roi_mask, cid in zip(axes[:-2], roi_mask_list, cids):
ax.imshow(roi_mask.get_mask_plane(), cmap='gray')
ax.set_title('cell %d' % cid)
# make a mask of all ROIs in the experiment
all_roi_masks = data_set.get_roi_mask_array()
combined_mask = all_roi_masks.max(axis=0)
axes[-2].imshow(combined_mask, cmap='gray')
axes[-2].set_title('all ROIs')
# show the movie max projection
max_projection = data_set.get_max_projection()
axes[-1].imshow(max_projection, cmap='gray')
axes[-1].set_title('max projection')
plt.show()
ROI Analysis¶
The code to compute all of the cell metrics available for download (along with several others we haven't put in our database yet) is in the AllenSDK. All you need is an NWB file. For example, this is how you compute metrics for the drifting gratings stimulus and plot some results.
from allensdk.brain_observatory.drifting_gratings import DriftingGratings
# example loading drifing grating data
data_set = boc.get_ophys_experiment_data(512326618)
dg = DriftingGratings(data_set)
dg_peak = dg.peak
print("done analyzing drifting gratings")
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
# filter for visually responding, selective cells
vis_cells = (dg_peak.ptest_dg < 0.05) & (dg_peak.peak_dff_dg > 3)
osi_cells = vis_cells & (dg_peak.osi_dg > 0.5) & (dg_peak.osi_dg <= 1.5)
dsi_cells = vis_cells & (dg_peak.dsi_dg > 0.5) & (dg_peak.dsi_dg <= 1.5)
# 2-d tf vs. ori histogram
# tfval = 0 is used for the blank sweep, so we are ignoring it here
os = np.zeros((len(dg.orivals), len(dg.tfvals)-1))
ds = np.zeros((len(dg.orivals), len(dg.tfvals)-1))
for i,trial in dg_peak[osi_cells].iterrows():
os[trial.ori_dg, trial.tf_dg-1] += 1
for i,trial in dg_peak[dsi_cells].iterrows():
ds[trial.ori_dg, trial.tf_dg-1] += 1
max_count = max(os.max(), ds.max())
fig, (ax1, ax2) = plt.subplots(1,2)
# plot direction selectivity
im = ax1.imshow(ds, clim=[0,max_count], cmap='hot', interpolation='nearest')
ax1.set_xlabel('temporal frequency')
ax1.set_ylabel('direction')
ax1.set_xticks(np.arange(len(dg.tfvals)-1))
ax1.set_xticklabels(dg.tfvals[1:])
ax1.set_yticks(np.arange(len(dg.orivals)))
ax1.set_yticklabels(dg.orivals)
ax1.set_title('direction selective cells')
# plot orientation selectivity
im = ax2.imshow(os, clim=[0,max_count], cmap='hot', interpolation='nearest')
ax2.set_xlabel('temporal frequency')
ax2.set_ylabel('orientation')
ax2.set_xticks(np.arange(len(dg.tfvals)-1))
ax2.set_xticklabels(dg.tfvals[1:])
ax2.set_yticks(np.arange(len(dg.orivals)))
ax2.set_yticklabels(dg.orivals)
ax2.set_title('orientation selective cells')
# plot a colorbar
fig.subplots_adjust(right=0.9)
cbar_ax = fig.add_axes([0.95, 0.05, 0.05, 0.85])
cbar = fig.colorbar(im, cax=cbar_ax)
cbar.set_ticks(np.arange(0, max_count, 2)+0.5)
cbar.set_ticklabels(np.arange(0, max_count, 2, dtype=int))
plt.show()
Neuropil Correction¶
All of the raw fluorescence traces are available in NWB files, but some of these signals are contaminated by nearby neuropil signal. Neuropil correction is performed after the signal from overlapping traces has been demixed. The code to perform neuropil correction is available in the AllenSDK and can be used as follows:
from allensdk.brain_observatory.r_neuropil import estimate_contamination_ratios
data_set = boc.get_ophys_experiment_data(569407590)
csid = data_set.get_cell_specimen_ids()[0]
time, demixed_traces = data_set.get_demixed_traces(cell_specimen_ids=[csid])
_, neuropil_traces = data_set.get_neuropil_traces(cell_specimen_ids=[csid])
results = estimate_contamination_ratios(demixed_traces[0], neuropil_traces[0])
correction = demixed_traces[0] - results['r'] * neuropil_traces[0]
print("r = %f" % results['r'])
print("max error = %f" % results['err'])
The NWB files contain the neuropil traces and r values needed to perform this correction, so you don't need to recompute them. The corrected trace can be computed on the fly when you use the get_corrected_fluorescence_traces method:
_, corrected_traces = data_set.get_corrected_fluorescence_traces(cell_specimen_ids=[csid])
plt.figure(figsize=(14,4))
plt.title("Neuropil-corrected Fluorescence Trace")
plt.plot(time, corrected_traces[0])
plt.show()
Compute dF/F¶
You can compute dF/F for yourself using the allensdk.brain_observatory.dff module.
from allensdk.brain_observatory.dff import compute_dff
plt.figure(figsize=(14,4))
plt.title("dF/F Trace")
dff = compute_dff(np.array(corrected_traces))
plt.plot(time, dff[0,:])
plt.show()
Running Speed¶
We recorded the animal's running speed during the course of the experiment and made it available in the NWB file.
from matplotlib import pyplot as plt
%matplotlib inline
data_set = boc.get_ophys_experiment_data(501940850)
dxcm, dxtime = data_set.get_running_speed()
plt.figure(figsize=(14,4))
plt.plot(dxtime, dxcm)
plt.show()
Motion Correction¶
X and Y translation values in pixels required to correct for motion artifacts during the experiment are available as well.
from matplotlib import pyplot as plt
%matplotlib inline
data_set = boc.get_ophys_experiment_data(501940850)
mc = data_set.get_motion_correction()
plt.figure(figsize=(14,4))
plt.plot(mc.timestamp, mc.x_motion)
plt.plot(mc.timestamp, mc.y_motion)
plt.legend(['x motion','y motion'])
plt.show()
Eye Tracking¶
Many experiments contain pupil position and pupil size from eye tracking. Extreme outliers from the tracking have been removed and replaced with NaN, which will appear as gaps in the plotted data. If an experiment does not have eye tracking data, a NoEyeTrackingException will be raised.
from allensdk.brain_observatory.brain_observatory_exceptions import NoEyeTrackingException
from matplotlib import pyplot as plt
%matplotlib inline
# example with no eye tracking data
data_set = boc.get_ophys_experiment_data(501940850)
try:
timestamps, locations = data_set.get_pupil_location()
except NoEyeTrackingException:
print("No eye tracking for experiment %s." % data_set.get_metadata()["ophys_experiment_id"])
data_set = boc.get_ophys_experiment_data(569407590)
# looking at azimuth and altitude over time
# by default locations returned are (azimuth, altitude)
# passing as_spherical=False to get_pupil_location will return (x,y) in cm
timestamps, locations = data_set.get_pupil_location()
plt.figure(figsize=(14,4))
plt.plot(timestamps, locations.T[0])
plt.plot(timestamps, locations.T[1])
plt.title("Eye position over time")
plt.xlabel("time (s)")
plt.ylabel("angle (deg)")
plt.legend(['azimuth', 'altitude'])
plt.show()
#pupil size over time
timestamps, area = data_set.get_pupil_size()
plt.figure(figsize=(14,4))
plt.plot(timestamps, area)
plt.title("Pupil size over time")
plt.xlabel("time (s)")
plt.ylabel("area (px)")
plt.ylim(0, 20000)
plt.show()
# scatter of gaze positions over approximate screen area
plt.figure()
plt.scatter(locations.T[0], locations.T[1], s=2, c="m", edgecolor="")
plt.title("Eye position scatter plot")
plt.xlim(-70, 70)
plt.ylim(-60, 60)
plt.xlabel("azimuth (deg)")
plt.ylabel("altitude (deg)")
plt.show()
