Brain Observatory Trace Analysis

This notebook demonstrates how to run the stimulus-specific tuning analysis code in the SDK. First let's instantiate a BrainObservatoryCache instance.

Download this notebook in .ipynb format here.

from pathlib import Path

output_dir = '.'
resources_dir = Path.cwd().parent / 'resources'
RUN_LOCALLY_SPARSE_NOISE = True

from allensdk.core.brain_observatory_cache import BrainObservatoryCache
from pathlib import Path

boc =  BrainObservatoryCache(
    manifest_file=str(Path(output_dir) / 'brain_observatory_manifest.json'))

Drifting Gratings

In this example, we'll show how you can plot a heatmap of a cell's response organized by orientation and temporal frequency. Here we start with a known experiment ID. Take a look at the other notebook to see how you can find experiments of interest. You can run the drifting grating analysis code on that experiment's NWB file as follows:

import matplotlib.pyplot as plt
%matplotlib inline

from allensdk.brain_observatory.drifting_gratings import DriftingGratings

data_set = boc.get_ophys_experiment_data(502376461)
dg = DriftingGratings(data_set)

2023-11-30 07:08:42,394 allensdk.api.api.retrieve_file_over_http INFO     Downloading URL: http://api.brain-map.org/api/v2/well_known_file_download/514609062

If you know which cell you're interested in, here's how you can find out where it is in the NWB File.

import numpy as np
specimen_id = 517425074
cell_loc = data_set.get_cell_specimen_indices([specimen_id])[0]

print("Specimen ID:", specimen_id)
print("Cell loc:", cell_loc)

Specimen ID: 517425074
Cell loc: 97

The response property of the stimulus-specific analysis objects is 4-D array organized with the following dimensions:

0: num. grating directions
1: num. grating temporal frequencies + 1 (0=blank sweep)
2: num. cells + 1 (running speed)
3: 0=response mean, 1=response standard error of the mean, 2=number of signficant trials

Dimension 2 of the response array has one index per cell in the experiment, plus one. The final index of that dimension is the running speed (response[:,:,-1,:]). This organization allows users to examine whether the mouse ran more for some specific stimulus conditions.

# skip the blank sweep column of the temporal frequency dimension
plt.imshow(dg.response[:,1:,cell_loc,0], cmap='hot', interpolation='none')
plt.xticks(range(5), dg.tfvals[1:])
plt.yticks(range(8), dg.orivals)
plt.xlabel("Temporal frequency (Hz)", fontsize=20)
plt.ylabel("Direction (deg)", fontsize=20)
plt.tick_params(labelsize=14)
cbar= plt.colorbar()
cbar.set_label("DF/F (%)")

The peak property of the analysis object is a Pandas DataFrame of peak conditions (direction and temporal frequency) as well as computed response metrics. For drifting gratings this includes:

ori_dg: preferred direction (index into dg.orivals)
tf_dg: preferred temporal frequency (index into tf.tfvals)
response_reliability_dg: response reliability 
osi_dg: orientation selectivity index
dsi_dg: direction selectivity index
ptest_dg: number of signficant cells
p_run_dg: K-S statistic comparing running trials to stationary trials
run_modulation_dg: ratio of mean fluorescence during running vs static
cv_dg: circular variance 

dg.peak.loc[cell_loc]

/home/runner/work/AllenSDK/AllenSDK/allensdk/brain_observatory/drifting_gratings.py:203: RuntimeWarning: invalid value encountered in double_scalars
  peak.cv_os_dg.iloc[nc] = np.abs(CV_top_os.sum()) / tuning.sum()
/home/runner/work/AllenSDK/AllenSDK/allensdk/brain_observatory/drifting_gratings.py:204: RuntimeWarning: invalid value encountered in double_scalars
  peak.cv_ds_dg.iloc[nc] = np.abs(CV_top_ds.sum()) / tuning.sum()

ori_dg                       4
tf_dg                        2
reliability_dg        0.186626
osi_dg                1.202353
dsi_dg                0.358697
peak_dff_dg           12.09474
ptest_dg                   0.0
p_run_dg              0.061258
run_modulation_dg    -0.620711
cv_os_dg                   1.0
cv_ds_dg              0.358697
tf_index_dg            0.30792
cell_specimen_id     517425074
Name: 97, dtype: object

Next let's plot all trials for a given cell's preferred condition.

pref_ori = dg.orivals[dg.peak.ori_dg[cell_loc]]
pref_tf = dg.tfvals[dg.peak.tf_dg[cell_loc]]
print("Preferred direction:", pref_ori)
print("Preferred temporal frequency:", pref_tf)

Preferred direction: 180
Preferred temporal frequency: 2

pref_trials = dg.stim_table[(dg.stim_table.orientation==pref_ori)&(dg.stim_table.temporal_frequency==pref_tf)]
pref_trials

	temporal_frequency	orientation	blank_sweep	start	end
1	2.0	180.0	0.0	837	897
71	2.0	180.0	0.0	7157	7217
73	2.0	180.0	0.0	7338	7398
100	2.0	180.0	0.0	9776	9835
141	2.0	180.0	0.0	13478	13537
175	2.0	180.0	0.0	16549	16609
271	2.0	180.0	0.0	55014	55074
291	2.0	180.0	0.0	56820	56879
323	2.0	180.0	0.0	59709	59769
435	2.0	180.0	0.0	97819	97879
449	2.0	180.0	0.0	99083	99143
488	2.0	180.0	0.0	102604	102664
509	2.0	180.0	0.0	104500	104560
518	2.0	180.0	0.0	105313	105372
606	2.0	180.0	0.0	113259	113319

sweep_response is a DataFrame that contains the DF/F of each cell during each stimulus trial. It shares its index with stim_table. Each cell contains a timeseries that extends from 1 second prior to the start of the trial to 1 second after the end of the trial. The final column of sweep_response, named dx, is the running speed of the mouse during each trial. The data in this DataFrame is used to create another DataFrame called mean_sweep_response that contains the mean DF/F during the trial for each cell (and the mean running speed in the last column).

subset = dg.sweep_response[(dg.stim_table.orientation==pref_ori)&(dg.stim_table.temporal_frequency==pref_tf)]

Here we look at the mean running speed during trials that presented the preferred condition.

subset_mean = dg.mean_sweep_response[(dg.stim_table.orientation==pref_ori)&(dg.stim_table.temporal_frequency==pref_tf)]
subset_mean['dx']

1       0.920868
71      0.060374
73      0.034393
100     4.897258
141    -0.000856
175     0.018537
271    38.012079
291     0.000437
323    36.196411
435    -0.012866
449    11.234307
488     5.972907
509     1.967057
518    -0.004841
606     6.747042
Name: dx, dtype: float64

Plot the response to each trial of the preferred condition, labeled with the mean running speed during the trial

trial_timestamps = np.arange(-1*dg.interlength, dg.interlength+dg.sweeplength, 1.)/dg.acquisition_rate
plt.figure(figsize=(8,20))
for i in range(len(subset)):
    plt.subplot(len(pref_trials),1,i+1)
    plt.plot(trial_timestamps, subset[str(cell_loc)].iloc[i], color='k', lw=2)
    plt.axvspan(0,2,color='red', alpha=0.3)
    plt.ylabel("DF/F (%)")
    plt.ylim(-10,600)
    plt.yticks(range(0,700,200))
    plt.text(2.5, 300, str(round(subset_mean['dx'].iloc[i],2))+" cm/s")
    if i<(len(subset)-1):
        plt.xticks([])
    else:
        plt.xticks([-1,0,1,2,3])
        plt.xlabel("Time (s)")

Static Gratings

The static gratings analysis object is quite similar to the drifting gratings analysis object. Here we'll just take a look at the peak table, which contains information about the preferred orientation, spatial frequency, phase, as well as a number of other metrics.

from allensdk.brain_observatory.static_gratings import StaticGratings

# example loading drifting grating data
data_set = boc.get_ophys_experiment_data(510938357)

sg = StaticGratings(data_set)

2023-11-30 07:09:48,431 allensdk.api.api.retrieve_file_over_http INFO     Downloading URL: http://api.brain-map.org/api/v2/well_known_file_download/516802541

sg.peak.head()

	ori_sg	sf_sg	phase_sg	reliability_sg	osi_sg	peak_dff_sg	ptest_sg	time_to_peak_sg	cell_specimen_id	p_run_sg	cv_os_sg	run_modulation_sg	sf_index_sg
0	3	5	1	-0.010099	0.749948	3.406017	0.516874	0.46438	517399188	0.286915	0.49489	-0.826111	0.14186
1	5	3	3	0.025564	0.893034	18.603113	0.0	0.43121	517399994	0.063211	0.816515	0.774433	0.280605
2	4	3	3	0.018425	0.956513	4.879125	0.753551	0.3317	517399857	0.110951	0.590487	0.886165	0.214098
3	5	5	1	-0.003168	1.03934	2.191308	0.098264	0.3317	517399727	0.285707	0.927832	1.042196	0.138216
4	3	5	0	-0.009008	0.314582	3.980433	0.068358	0.82925	517399442	0.120239	0.258181	0.958413	0.201612

Natural Scenes

The natural scenes analysis object is again similar to the others. In addition to computing the sweep_response and mean_sweep_response arrays, NaturalScenes reports the cell's preferred scene, running modulation, time to peak response, and other metrics.

from allensdk.brain_observatory.natural_scenes import NaturalScenes

data_set = boc.get_ophys_experiment_data(510938357)

ns = NaturalScenes(data_set)
print("done analyzing natural scenes")

done analyzing natural scenes

ns.peak.head()

	scene_ns	reliability_ns	peak_dff_ns	ptest_ns	p_run_ns	run_modulation_ns	time_to_peak_ns	cell_specimen_id	image_selectivity_ns
0	48	0.006932	4.961453	0.218335	0.060267	-1.151861	0.72974	517399188	0.325864
1	96	0.132598	8.396541	0.0	0.010665	-1.228357	0.3317	517399994	0.553542
2	103	-0.008523	2.640692	0.013945	0.096297	1.008399	0.29853	517399857	-0.00278
3	103	0.010471	1.747232	0.589355	0.239155	1.033275	0.3317	517399727	0.247847
4	21	-0.000713	1.511358	0.059133	0.427964	0.53565	0.53072	517399442	0.280356

Locally Sparse Noise

The locally sparse noise stimulus object is a bit different from the others. It does not have a peak condition table, instead providing a method to retrieve the "on" and "off" receptive fields of all cells. The receptive field of a cell is computed by averaging responses to trials in which a given sparse noise grid location is on/off.

from allensdk.brain_observatory.locally_sparse_noise import LocallySparseNoise
import allensdk.brain_observatory.stimulus_info as stim_info

specimen_id = 587179530
cell = boc.get_cell_specimens(ids=[specimen_id])[0]

exp = boc.get_ophys_experiments(experiment_container_ids=[cell['experiment_container_id']],
                                stimuli=[stim_info.LOCALLY_SPARSE_NOISE])[0]
                                         
data_set = boc.get_ophys_experiment_data(exp['id'])

2023-11-30 07:22:53,855 allensdk.api.api.retrieve_file_over_http INFO     Downloading URL: http://api.brain-map.org/api/v2/well_known_file_download/514409433

import matplotlib.pyplot as plt
from IPython.display import display
from PIL import Image
%matplotlib inline

if RUN_LOCALLY_SPARSE_NOISE:
    lsn = LocallySparseNoise(data_set)
    print("done analyzing locally sparse noise")

    cell_idx = data_set.get_cell_specimen_indices([specimen_id])[0]

    plt.imshow(lsn.receptive_field[:,:,cell_idx,0], interpolation='nearest', cmap='PuRd', origin='lower')
    plt.title("on receptive field")
    plt.show()
    plt.imshow(lsn.receptive_field[:,:,cell_idx,1], interpolation='nearest', cmap='Blues', origin='lower')
    plt.title("off receptive field")
    plt.show()
else:
    # read in LocallySparseNoise output and display it
    off_receptive_field = Image.open(Path(resources_dir) /
                                     'brain_observatory_analysis' /
                                     'off_receptive_field.png').convert('RGB')
    display(off_receptive_field)

    on_receptive_field = Image.open(Path(resources_dir) /
                                     'brain_observatory_analysis' /
                                     'on_receptive_field.png').convert('RGB')
    display(on_receptive_field)