# install cell skipped during CI (deps preinstalled into system Python)

import os
import json
import pathlib
import warnings
from typing import Union, Tuple, Literal

import h5py
import matplotlib
import numpy
import pandas
import pynwb
import remfile
import scipy.stats as scipy_stats
import matplotlib.pyplot as matplotlib_pyplot
from dandi.dandiapi import DandiAPIClient, RemoteAsset, RemoteDandiset

%matplotlib inline

# You may need to change these paths to the files depending on your current working directory
# These defaults assume the current working directory is the folder containing this notebook 

asset_metadata_folder_path = pathlib.Path("asset_metadata")
experiment_containers_metadata_file_path = asset_metadata_folder_path / "experiment_containers_metadata.json"
visual_coding_sessions_metadata_file_path = asset_metadata_folder_path / "visual_coding_sessions_metadata.json"

with open(file=experiment_containers_metadata_file_path, mode="r") as fp:
    experiment_containers_metadata = json.load(fp=fp)
with open(file=visual_coding_sessions_metadata_file_path, mode="r") as fp:
    visual_coding_sessions_metadata = json.load(fp=fp)

all_targeted_structures = list(set([experiment_container["targeted_structure"] for experiment_container in experiment_containers_metadata]))
all_targeted_structures

['VISam', 'VISp', 'VISl', 'VISpm', 'VISrl', 'VISal']

all_cre_lines = list(set([experiment_container["cre_line"] for experiment_container in experiment_containers_metadata]))
all_cre_lines

['Ntsr1-Cre_GN220',
 'Slc17a7-IRES2-Cre',
 'Pvalb-IRES-Cre',
 'Nr5a1-Cre',
 'Emx1-IRES-Cre',
 'Rbp4-Cre_KL100',
 'Scnn1a-Tg3-Cre',
 'Cux2-CreERT2',
 'Fezf2-CreER',
 'Vip-IRES-Cre',
 'Tlx3-Cre_PL56',
 'Rorb-IRES2-Cre',
 'Sst-IRES-Cre']

all_reporter_lines = list(set([experiment_container["reporter_line"] for experiment_container in experiment_containers_metadata]))
all_reporter_lines

['Ai148(TIT2L-GC6f-ICL-tTA2)',
 'Ai94(TITL-GCaMP6s)',
 'Ai93(TITL-GCaMP6f)-hyg',
 'Ai162(TIT2L-GC6s-ICL-tTA2)',
 'Ai93(TITL-GCaMP6f)']

all_session_types = list(set([session["session_type"] for session in visual_coding_sessions_metadata]))
all_session_types

['three_session_B', 'three_session_C2', 'three_session_C', 'three_session_A']

experiments_subset = [
    experiment_container
    for experiment_container in experiment_containers_metadata
    if experiment_container["targeted_structure"] == "VISp" and experiment_container["cre_line"] == "Cux2-CreERT2"
]
experiment_subset_table = pandas.DataFrame(experiments_subset)
experiment_subset_table

number_of_experiment_containers = len(experiment_containers_metadata)
number_of_experiment_containers

540

experiment_container_id = experiments_subset[0]["id"]
experiment_container_id

511510736

sessions = [
    visual_coding_session for visual_coding_session in visual_coding_sessions_metadata
    if visual_coding_session["experiment_container_id"] == experiment_container_id
]
sessions_table = pandas.DataFrame(sessions)
sessions_table

session_id = [session for session in sessions if session["session_type"] == "three_session_B"][0]["id"]
session_id

501559087

def read_dandi_file(*, asset: RemoteAsset) -> pynwb.NWBFile:
    """Remotely load the NWB file from a DANDI asset object."""
    # Harmless warnings from PyNWB
    warnings.filterwarnings(action="ignore", message="No cached namespaces found in .*")
    warnings.filterwarnings(action="ignore", message="Ignoring cached namespace .*")
    
    s3_url = asset.get_content_url()
    byte_stream = remfile.File(url=s3_url)
    file = h5py.File(name=byte_stream)
    io = pynwb.NWBHDF5IO(file=file, mode="r", load_namespaces=True)
    nwbfile = io.read()

    return nwbfile

def read_visual_coding_ophys_file(
    *, dandiset: RemoteDandiset, session_id: str, raw_or_processed: Literal["raw", "processed"] = "processed"
) -> pynwb.NWBFile:
    """Load the NWB file from DANDI given only the unique session ID in the Visual Coding - Optical Physiology dataset."""
    if raw_or_processed == "processed":
        asset = next(asset for asset in dandiset.get_assets() if f"ses-{session_id}" in asset.path and "behavior" in asset.path)
    elif raw_or_processed == "raw":
        asset = next(asset for asset in dandiset.get_assets() if f"ses-{session_id}" in asset.path and "behavior" not in asset.path)
    return read_dandi_file(asset=asset)

dandiset_id = "000728"

client = DandiAPIClient()
dandiset = client.get_dandiset(dandiset_id=dandiset_id)

nwbfile = read_visual_coding_ophys_file(dandiset=dandiset, session_id=session_id)

nwbfile

def plot_image(
    *,  # "*" character forces keyword argument usage, which makes code much more readable
    image: numpy.ndarray,
    size: int = 8,
    color_map: Union[str, None] = None
) -> None:
    """
    Pretty rendering of a simple square grayscale image.
    
    Axis labels and ticks are stripped since their size is viewable from the image array shape.
    Also, display the color bar at a good size if color is not gray.
    """
    matplotlib_pyplot.figure(figsize=(size, size))

    frame = matplotlib_pyplot.gca()
    for label in frame.axes.get_xticklabels():
        label.set_visible(False)
        label.set_fontsize(0.0)
    for label in frame.axes.get_yticklabels():
        label.set_fontsize(0.0)
        label.set_visible(False)
    for tick in frame.axes.get_xticklines():
        tick.set_visible(False)
    for tick in frame.axes.get_yticklines():
        tick.set_visible(False)
    
    plot = matplotlib_pyplot.imshow(image, cmap=color_map);

    if color_map != "gray":
        color_bar = matplotlib_pyplot.colorbar(plot, fraction=0.046, pad=0.04)
        color_bar.solids.set_edgecolor("face")

def plot_traces(
    *,
    timestamps: numpy.ndarray,
    traces: numpy.ndarray,
    hide_y_scale: bool = True,
    x_axis_label: Union[str, None] = "Time (s)",  # Since timestamps in seconds are so common in NWB, just make it the default
    y_axis_label: Union[str, None] = None,
    padding: float = 2.0,
    trace_color: Union[str, None] = None,
    reset_figure: bool = True,
    figure_shape: Tuple[int, int] = (10, 4),
    return_frame: bool = False,
) -> Union[matplotlib.axes.Axes, None]:
    """
    Pretty rendering of traces of activity (such as ROI or running speed) over time.

    The Y axis of ROI traces is rarely known in practice, so numerical scale doesn't really matter for plotting.
    Also, timestamps from NWB will always be in seconds.
    """
    if reset_figure:
        matplotlib_pyplot.figure(figsize=figure_shape)

    if len(traces.shape) == 1:
        matplotlib_pyplot.plot(timestamps, traces, color=trace_color)
    else:
        for trace_index in range(traces.shape[1]):
            matplotlib_pyplot.plot(
                timestamps,
                traces[:, trace_index] + (trace_index * padding), # Add vertical padding to separate each trace
                color=trace_color,
            )

    frame = matplotlib_pyplot.gca()
    frame.axes.spines['top'].set_visible(False)
    frame.axes.spines['right'].set_visible(False)
    frame.axes.spines['bottom'].set_visible(False)
    frame.axes.spines['left'].set_visible(False)

    if hide_y_scale:
        for label in frame.axes.get_yticklabels():
            label.set_fontsize(0.0)
            label.set_visible(False)
        for tick in frame.axes.get_yticklines():
            tick.set_visible(False)

    matplotlib_pyplot.xlabel(x_axis_label, fontsize=16)
    if y_axis_label is not None:
        matplotlib_pyplot.ylabel(y_axis_label, fontsize=16)

    if return_frame:
        return frame

def apply_epoch_shading(
    *,
    frame: matplotlib.axes.Axes,
    epoch_table: pynwb.file.TimeIntervals,
) -> None:
    """For each epoch, shade the plot based on the stimulus type."""
    # Clear previous legend to prevent appending
    frame.legend([])
    
    epoch_table_data_frame = epoch_table.to_dataframe()

    stimulus_type_to_color = dict(
        spontaneous="grey",
        static_gratings="fuchsia",
        drifting_gratings="lavender",
        natural_scenes="orange",
        natural_movie_one="green",
        natural_movie_two="red",
        natural_movie_three="blue",
    )
    stimulus_type_to_label = dict(
        spontaneous="Spontaneous",
        static_gratings="Static Gratings",
        drifting_gratings="Drifting Gratings",
        natural_scenes="Natural Scenes",
        natural_movie_one="Natural Movie (1)",
        natural_movie_two="Natural Movie (2)",
        natural_movie_three="Natural Movie (3)",
    )

    patches= list()
    for stimulus_name in epoch_table_data_frame.stimulus_type.unique():
        stimulus = epoch_table_data_frame[epoch_table_data_frame.stimulus_type==stimulus_name]

        if len(stimulus) > 0:
            patches.append(
                matplotlib.patches.Patch(
                    color=stimulus_type_to_color[stimulus_name],
                    label=stimulus_type_to_label[stimulus_name],
                    alpha=0.1,
                )
            )

        for _, row in stimulus.iterrows():
            frame.axvspan(
                xmin=row.start_time, xmax=row.stop_time, color=stimulus_type_to_color[stimulus_name], alpha=0.1
            )

    frame.legend(title="Epochs", loc="center left", handles=patches, shadow=True, fancybox=True, bbox_to_anchor=(1, 0.5));

nwbfile_with_raw_imaging = read_visual_coding_ophys_file(dandiset=dandiset, session_id=session_id, raw_or_processed="raw")
nwbfile_with_raw_imaging

motion_corrected_imaging = nwbfile_with_raw_imaging.acquisition["MotionCorrectedTwoPhotonSeries"]
motion_corrected_imaging

plot_image(image=motion_corrected_imaging.data[3_000, ...], color_map="gray")

maximum_intensity_projection = nwbfile.processing["ophys"]["SummaryImages"]["maximum_intensity_projection"][:]
maximum_intensity_projection

array([[ 296,  308,  306, ..., 2760, 2164,  888],
       [ 312,  298,  303, ..., 2741, 1716, 1031],
       [ 278,  295,  298, ..., 1426,  929,  574],
       ...,
       [ 347,  348,  366, ...,  407,  403,  410],
       [ 335,  357,  398, ...,  317,  354,  362],
       [ 308,  288,  306, ...,  293,  316,  318]], dtype=uint16)

maximum_intensity_projection.shape

(512, 512)

plot_image(image=maximum_intensity_projection, color_map="gray")

plane_segmentation = nwbfile.processing["ophys"]["ImageSegmentation"]["PlaneSegmentation"]
plane_segmentation_table = plane_segmentation.to_dataframe()
plane_segmentation_table

pixel_masks = plane_segmentation["pixel_mask"][:]
pixel_masks[0][:10]

array([(55, 291, 1.), (55, 292, 1.), (56, 290, 1.), (56, 291, 1.),
       (56, 292, 1.), (56, 293, 1.), (56, 294, 1.), (56, 295, 1.),
       (56, 296, 1.), (57, 289, 1.)],
      dtype=[('x', '<u4'), ('y', '<u4'), ('weight', '<f4')])

combined_image_masks = numpy.zeros(shape=maximum_intensity_projection.shape)
for pixel_mask in plane_segmentation["pixel_mask"][:]:
    for x, y, w in pixel_mask:
        combined_image_masks[x,y] += w

plot_image(image=combined_image_masks)

nwbfile.processing["ophys"]["DfOverF"]

nwbfile.processing["ophys"]["Fluorescence"]

df_over_f = nwbfile.processing["ophys"]["DfOverF"]["DfOverF"].data
df_over_f_timestamps = nwbfile.processing["ophys"]["DfOverF"]["DfOverF"].timestamps[...]

roi_index = 0

plot_traces(timestamps=df_over_f_timestamps, traces=df_over_f[:, roi_index], y_axis_label = "$\Delta F \ / \ F$ (a.u.)")

<>:3: SyntaxWarning: invalid escape sequence '\D'
<>:3: SyntaxWarning: invalid escape sequence '\D'
/tmp/ipykernel_2280/934175159.py:3: SyntaxWarning: invalid escape sequence '\D'
  plot_traces(timestamps=df_over_f_timestamps, traces=df_over_f[:, roi_index], y_axis_label = "$\Delta F \ / \ F$ (a.u.)")

df_over_f_timestamps.shape

(113888,)

frame_slice = slice(10_000, 30_000)

plot_traces(timestamps=df_over_f_timestamps[frame_slice], traces=df_over_f[frame_slice, roi_index], y_axis_label = "$\Delta F \ / \ F$ (a.u.)")

<>:3: SyntaxWarning: invalid escape sequence '\D'
<>:3: SyntaxWarning: invalid escape sequence '\D'
/tmp/ipykernel_2280/933709331.py:3: SyntaxWarning: invalid escape sequence '\D'
  plot_traces(timestamps=df_over_f_timestamps[frame_slice], traces=df_over_f[frame_slice, roi_index], y_axis_label = "$\Delta F \ / \ F$ (a.u.)")

plot_traces(
    timestamps=df_over_f_timestamps,
    traces=df_over_f[:, 0:50],
    figure_shape=(10, 8),
    y_axis_label = "$\Delta F \ / \ F$ (a.u.)",
    trace_color="gray",
)

<>:5: SyntaxWarning: invalid escape sequence '\D'
<>:5: SyntaxWarning: invalid escape sequence '\D'
/tmp/ipykernel_2280/2496006190.py:5: SyntaxWarning: invalid escape sequence '\D'
  y_axis_label = "$\Delta F \ / \ F$ (a.u.)",

epoch_table = nwbfile.epochs
epoch_table

epoch_table.stimulus_type.data[:]

array(['static_gratings', 'natural_scenes', 'spontaneous',
       'natural_scenes', 'static_gratings', 'natural_movie_one',
       'natural_scenes', 'static_gratings'], dtype=object)

frame = plot_traces(
    timestamps=df_over_f_timestamps,
    traces=df_over_f[:, 0:50],
    figure_shape=(10, 8),
    y_axis_label = "$\Delta F \ / \ F$ (a.u.)",
    trace_color="gray",
    return_frame=True,
)

apply_epoch_shading(frame=frame, epoch_table=epoch_table)

<>:5: SyntaxWarning: invalid escape sequence '\D'
<>:5: SyntaxWarning: invalid escape sequence '\D'
/tmp/ipykernel_2280/3447667781.py:5: SyntaxWarning: invalid escape sequence '\D'
  y_axis_label = "$\Delta F \ / \ F$ (a.u.)",

running_speed = nwbfile.processing["behavior"]["BehavioralTimeSeries"]["running_speed"]
running_speed

plot_traces(timestamps=running_speed.timestamps[:], traces=running_speed.data[:], hide_y_scale=False, y_axis_label="Running speed (cm/s)")

plot_traces(
    timestamps=df_over_f_timestamps,
    traces=df_over_f[:, 0:50],
    figure_shape=(14,10),
    trace_color="gray",
)
frame = plot_traces(
    timestamps=running_speed.timestamps[:],
    traces=(0.2 * running_speed.data[:]) - 20,  # Apply some custom scaling to get both traces on the same x-axis
    reset_figure=False,
    return_frame=True,
)
    
apply_epoch_shading(frame=frame, epoch_table=epoch_table)

plot_traces(
    timestamps=df_over_f_timestamps,
    traces=df_over_f[:, 49],
    figure_shape=(14,4),
    trace_color="gray",
)
frame = plot_traces(
    timestamps=running_speed.timestamps[:],
    traces=(0.2 * running_speed.data[:]) - 20,  # Apply some custom scaling to get both traces on the same x-axis
    reset_figure=False,
    return_frame=True,
)
    
apply_epoch_shading(frame=frame, epoch_table=epoch_table)

plot_traces(
    timestamps=df_over_f_timestamps,
    traces=df_over_f[:, 4],
    figure_shape=(14,4),
    trace_color="gray",
)
frame = plot_traces(
    timestamps=running_speed.timestamps[:],
    traces=(0.2 * running_speed.data[:]) - 20,  # Apply some custom scaling to get both traces on the same x-axis
    reset_figure=False,
    return_frame=True,
)
    
apply_epoch_shading(frame=frame, epoch_table=epoch_table)

plot_traces(
    timestamps=df_over_f_timestamps,
    traces=df_over_f[:, 35],
    figure_shape=(14,4),
    trace_color="gray",
)
frame = plot_traces(
    timestamps=running_speed.timestamps[:],
    traces=(0.2 * running_speed.data[:]) - 20,  # Apply some custom scaling to get both traces on the same x-axis
    reset_figure=False,
    return_frame=True,
)
    
apply_epoch_shading(frame=frame, epoch_table=epoch_table)

df_over_f_events = nwbfile.processing["ophys"]["DfOverF"]["DfOverFEvents"]
df_over_f_events

plot_traces(
    timestamps=df_over_f_events.timestamps[:],
    traces=2 * df_over_f_events.data[:, 1],  # Arbitrarily scale by 2 for better appearance
)

plot_traces(
    timestamps=df_over_f_events.timestamps[:],
    traces=df_over_f[:, 1] - 5,  # Arbitrarily shift by 5 for better appearance
    reset_figure=False,
)

plot_traces(
    timestamps=df_over_f_events.timestamps[:],
    traces=2 * df_over_f_events.data[:, 1],  # Arbitrarily scale by 2 for better appearance
)

plot_traces(
    timestamps=df_over_f_events.timestamps[:],
    traces=df_over_f[:, 1] - 5,  # Arbitrarily shift by 5 for better appearance
    reset_figure=False,
)
matplotlib_pyplot.xlim(690, 720);  # in units seconds

plot_traces(
    timestamps=df_over_f_events.timestamps[:],
    traces=df_over_f_events.data[:, 0:50],
    figure_shape=(14,10),
    trace_color="gray",
)
frame = plot_traces(
    timestamps=running_speed.timestamps[:],
    traces=(0.2 * running_speed.data[:]) - 20,  # Apply some custom scaling to get both traces on the same x-axis
    reset_figure=False,
    return_frame=True,
)
    
apply_epoch_shading(frame=frame, epoch_table=epoch_table)

natural_scene_presentation = nwbfile.stimulus["natural_scenes_stimulus"]
natural_scene_presentation

natural_scene_presentation = nwbfile.stimulus["natural_scenes_stimulus"]
natural_scene_presentation.data[:]

array([81, 33, 76, ..., 34,  0, 87], dtype=uint32)

natural_scene_presentation.timestamps[:]

array([ 561.98893,  562.25491,  562.48764, ..., 3221.9665 , 3222.23249,
       3222.46524])

nwbfile.stimulus_template["natural_scenes_template"]

natural_scene_template = nwbfile.stimulus_template["natural_scenes_template"]
natural_scene_template

len(natural_scene_template.images)

118

scene_number = 101

plot_image(image=natural_scene_template[f"NaturalScene{scene_number}"], color_map="gray")

natural_movie_one_presentation = nwbfile.stimulus["natural_movie_one_stimulus"]
natural_movie_one_presentation

natural_movie_one_template = nwbfile.stimulus_template["natural_movie_one"]
natural_movie_one_template

natural_movie_one_template.data.shape

(900, 304, 608)

frame_index = 101

plot_image(image=natural_movie_one_template.data[frame_index, ...], color_map="gray")

static_gratings_table = nwbfile.stimulus["static_gratings"]
static_gratings_table.to_dataframe()

static_gratings_column_descriptions = {column.name: column.description for column in static_gratings_table.columns}
static_gratings_column_descriptions

{'start_time': 'Start time of epoch, in seconds',
 'stop_time': 'Stop time of epoch, in seconds',
 'orientation_in_degrees': 'Angle of the grating in degrees. NaN values correspond to a blank sweep.',
 'spatial_frequency_in_cycles_per_degree': 'Period of the grating in cycles/degree. NaN values correspond to a blank sweep.',
 'phase': 'Relative position of the grating. Phase 0 and Phase 0.5 are 180° apart so that the peak of the grating of phase 0 lines up with the trough of phase 0.5. NaN values correspond to a blank sweep.',
 'is_blank_sweep': 'Mean luminance gray image.'}

session_id_with_drifting_gratings = "501940850"

nwbfile_with_drifting_gratings = read_visual_coding_ophys_file(dandiset=dandiset, session_id=session_id_with_drifting_gratings)

drifting_gratings = nwbfile_with_drifting_gratings.stimulus["drifting_gratings"]
drifting_gratings

drifting_gratings_column_descriptions = {column.name: column.description for column in drifting_gratings.columns}
drifting_gratings_column_descriptions

{'start_time': 'Start time of epoch, in seconds',
 'stop_time': 'Stop time of epoch, in seconds',
 'orientation_in_degrees': 'Angle of the grating in degrees. NaN values correspond to a blank sweep.',
 'spatial_frequency_in_cycles_per_degree': 'Period of the grating in cycles/degree. NaN values correspond to a blank sweep.',
 'temporal_frequency_in_hz': 'The speed at which the grating moves in Hz. NaN values correspond to a blank sweep.',
 'is_blank_sweep': 'Mean luminance gray image.'}

df_over_f = df_over_f[...]

scene_number = 101
cell_index = 1

times_when_scene_was_presented = natural_scene_presentation.timestamps[natural_scene_presentation.data[:] == scene_number]

# Applying 30 frames of buffer room before and 50 after
df_over_f_times_to_frames = numpy.searchsorted(df_over_f_timestamps, times_when_scene_was_presented)

for df_over_f_time_to_frame in df_over_f_times_to_frames[:]:
    subset_df_over_f_frames = slice(df_over_f_time_to_frame - 30, df_over_f_time_to_frame + 50)

    timestamps = df_over_f_timestamps[subset_df_over_f_frames]
    aligned_timestamps = timestamps - timestamps[30]
    
    frame = plot_traces(
        timestamps=aligned_timestamps,
        traces=df_over_f[subset_df_over_f_frames, cell_index],
        x_axis_label="Time aligned to stimulus onset (s)",
        y_axis_label="$\Delta F \ / \ F$ (a.u.)",
        reset_figure=False,
        return_frame=True,
    )

# Here, we use the fact about duration of the stimulusu learned from the description of the NaturalScene objects
frame.axvspan(aligned_timestamps[30], aligned_timestamps[37], color="gray", alpha=0.3);

<>:19: SyntaxWarning: invalid escape sequence '\D'
<>:19: SyntaxWarning: invalid escape sequence '\D'
/tmp/ipykernel_2280/3822388140.py:19: SyntaxWarning: invalid escape sequence '\D'
  y_axis_label="$\Delta F \ / \ F$ (a.u.)",

scene_number = 101
cell_index = 1

times_when_scene_was_presented = natural_scene_presentation.timestamps[natural_scene_presentation.data[:] == scene_number]

# Applying 30 frames of buffer room before and 50 after
df_over_f_times_to_frames = numpy.searchsorted(df_over_f_timestamps, times_when_scene_was_presented)
timestamps = df_over_f_timestamps[subset_df_over_f_frames]
aligned_timestamps = timestamps - timestamps[30]

for presentation_index, df_over_f_time_to_frame in enumerate(df_over_f_times_to_frames[:], start=1):
    subset_df_over_f_frames = slice(df_over_f_time_to_frame - 30, df_over_f_time_to_frame + 50)

    frame = plot_traces(
        timestamps=aligned_timestamps,
        traces=2 * df_over_f_events.data[subset_df_over_f_frames, cell_index] + presentation_index,
        reset_figure=False,
        x_axis_label="Time aligned to stimulus onset (s)",
        y_axis_label="$\Delta F \ / \ F$ (a.u.)",
        return_frame=True,
    )

# Here, we use the fact about duration of the stimulusu learned from the description of the NaturalScene objects
frame.axvspan(aligned_timestamps[30], aligned_timestamps[37], color="gray", alpha=0.3);

<>:19: SyntaxWarning: invalid escape sequence '\D'
<>:19: SyntaxWarning: invalid escape sequence '\D'
/tmp/ipykernel_2280/3689056126.py:19: SyntaxWarning: invalid escape sequence '\D'
  y_axis_label="$\Delta F \ / \ F$ (a.u.)",

number_of_rois = df_over_f.shape[1]
number_of_natural_scenes = len(natural_scene_template.images)

response = numpy.empty(shape=(number_of_rois, number_of_natural_scenes))
for roi_index in range(number_of_rois):
    for scene_number in range(number_of_natural_scenes):
        times_when_scene_was_presented = natural_scene_presentation.timestamps[natural_scene_presentation.data[:] == scene_number]
        df_over_f_times_to_frames = numpy.searchsorted(df_over_f_timestamps, times_when_scene_was_presented)
        
        # Take the mean over a very short period just after stimulus onset (15 frames)
        frames_before_onset = 0
        frames_after_onset = 15
        temp = numpy.empty(shape=(frames_after_onset + frames_before_onset, len(times_when_scene_was_presented)))
        for temp_index, df_over_f_time_to_frame in enumerate(df_over_f_times_to_frames):
            subset_df_over_f_frames = slice(df_over_f_time_to_frame - frames_before_onset, df_over_f_time_to_frame + frames_after_onset)

            temp[:, temp_index] = df_over_f[subset_df_over_f_frames, roi_index]
        response[roi_index, scene_number] = temp.mean()

matplotlib_pyplot.plot(response[0,:], 'o-');

signal_correlations = numpy.empty(shape=(number_of_rois, number_of_rois))
for roi_index_1 in range(number_of_rois):
    for roi_index_2 in range(number_of_rois):
        r, p = scipy_stats.pearsonr(response[roi_index_1, :], response[roi_index_2, :])
        signal_correlations[roi_index_1, roi_index_2] = r

matplotlib_pyplot.imshow(signal_correlations);

pixel_masks = plane_segmentation["pixel_mask"][:]
centroid_x = numpy.zeros(shape=number_of_rois)
centroid_y = numpy.zeros(shape=number_of_rois)

for roi_index, pixel_mask in enumerate(pixel_masks):
    centroid_x[roi_index] = numpy.mean([pixel[0] for pixel in pixel_mask])
    centroid_y[roi_index] = numpy.mean([pixel[1] for pixel in pixel_mask])

distance = numpy.zeros(shape=(number_of_rois, number_of_rois))
for roi_index_1 in range(number_of_rois):
    for roi_index_2 in range(number_of_rois):
        distance[roi_index_1, roi_index_2] = numpy.sqrt(
            (centroid_x[roi_index_1] - centroid_x[roi_index_2])**2 + (centroid_y[roi_index_1] - centroid_y[roi_index_2])**2 
        )

matplotlib_pyplot.imshow(distance);

# Only use the upper half of the matrix due to symmetry of correlation
indices = numpy.triu_indices(number_of_rois, k=1)

matplotlib_pyplot.figure(figsize=(8,6))
matplotlib_pyplot.plot(distance[indices[0], indices[1]], signal_correlations[indices[0], indices[1]], ".", markersize=3.0)
matplotlib_pyplot.xlabel("Distance (pixels)", fontsize=16)
matplotlib_pyplot.ylabel("Signal correlation", fontsize=16);

Data type	uint32
Shape	(9000,)
Array size	35.16 KiB
Chunk shape	(9000,)
Compression	gzip
Compression opts	4
Uncompressed size (bytes)	36000
Compressed size (bytes)	1595
Compression ratio	22.570532915360502

Data type	float64
Shape	(9000,)
Array size	70.31 KiB
Chunk shape	(9000,)
Compression	gzip
Compression opts	4
Uncompressed size (bytes)	72000
Compressed size (bytes)	43200
Compression ratio	1.6666666666666667

Data type	uint8
Shape	(900, 304, 608)
Array size	158.64 MiB
Chunk shape	(270, 135, 270)
Compression	gzip
Compression opts	4
Uncompressed size (bytes)	166348800
Compressed size (bytes)	56192052
Compression ratio	2.9603617251777883

Data type	uint32
Shape	(5950,)
Array size	23.24 KiB
Chunk shape	(5950,)
Compression	gzip
Compression opts	4
Uncompressed size (bytes)	23800
Compressed size (bytes)	8381
Compression ratio	2.839756592292089

Data type	float64
Shape	(5950,)
Array size	46.48 KiB
Chunk shape	(5950,)
Compression	gzip
Compression opts	4
Uncompressed size (bytes)	47600
Compressed size (bytes)	35597
Compression ratio	1.3371913363485688

Installing requirements¶

Navigating dandiset 000728: Visual Coding - Optical Physiology

Part 1: Getting oriented to the dataset

1.1 Experiment containers & sessions¶

Part 2: Streaming all the data for a single session

root (NWBFile)

Part 3: Loading specific fields of the data for a single session

3.0 Raw Motion Corrected Movies¶

root (NWBFile)

MotionCorrectedTwoPhotonSeries (TwoPhotonSeries)

3.1 Maximum projection¶

3.2 PlaneSegmentation¶

3.3 Segmented and Processed Traces¶

DfOverF

Fluorescence

3.4 Stimulus epochs¶

epochs (TimeIntervals)

3.5 Running speed¶

running_speed (TimeSeries)

Interesting things¶

3.6 Extracted events from $\Delta F / F$¶

DfOverFEvents (RoiResponseSeries)

3.7 Stimulus Module¶

3.7a Natural Scenes¶

natural_scenes_stimulus (IndexSeries)

natural_scenes_template (Images)

natural_scenes_template (Images)

3.7b Natural Movies¶

natural_movie_one_stimulus (IndexSeries)

natural_movie_one (ImageSeries)

3.7c Static Gratings¶

3.7 Drifting Gratings¶

drifting_gratings (TimeIntervals)

Part 4: Compute activity aligned to stimulus presentations

Part 5: Compute the signal correlations as a function of distance between neurons

Conclusion

	id	imaging_depth	targeted_structure	cre_line	reporter_line	donor_name	specimen_name	tags	failed
0	511510736	175	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	222426	Cux2-CreERT2;Camk2a-tTA;Ai93-222426	[]	False
1	511510855	175	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	229106	Cux2-CreERT2;Camk2a-tTA;Ai93-229106	[]	False
2	511509529	175	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	222420	Cux2-CreERT2;Camk2a-tTA;Ai93-222420	[]	False
3	511507650	175	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	222424	Cux2-CreERT2;Camk2a-tTA;Ai93-222424	[]	False
4	702934962	275	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	382421	Cux2-CreERT2;Camk2a-tTA;Ai93-382421	[]	False
5	645413757	275	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	348262	Cux2-CreERT2;Camk2a-tTA;Ai93-348262	[]	False
6	659767480	275	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	360565	Cux2-CreERT2;Camk2a-tTA;Ai93-360565	[]	False
7	511510650	175	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	222425	Cux2-CreERT2;Camk2a-tTA;Ai93-222425	[]	False
8	712178509	275	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	390323	Cux2-CreERT2;Camk2a-tTA;Ai93-390323	[]	False
9	511510667	275	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	222420	Cux2-CreERT2;Camk2a-tTA;Ai93-222420	[]	False
10	524691282	275	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	243293	Cux2-CreERT2;Camk2a-tTA;Ai93-243293	[]	False
11	701412138	175	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	382421	Cux2-CreERT2;Camk2a-tTA;Ai93-382421	[]	False
12	511510718	175	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	231584	Cux2-CreERT2;Camk2a-tTA;Ai93-231584	[]	False
13	511510699	275	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	225037	Cux2-CreERT2;Camk2a-tTA;Ai93-225037	[]	False
14	511510779	275	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	225036	Cux2-CreERT2;Camk2a-tTA;Ai93-225036	[]	False
15	511510670	175	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	225037	Cux2-CreERT2;Camk2a-tTA;Ai93-225037	[]	False

	id	imaging_depth	targeted_structure	cre_line	reporter_line	acquisition_age_days	experiment_container_id	session_type	donor_name	specimen_name	fail_eye_tracking
0	501559087	175	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	103	511510736	three_session_B	222426	Cux2-CreERT2;Camk2a-tTA;Ai93-222426	True
1	501704220	175	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	104	511510736	three_session_A	222426	Cux2-CreERT2;Camk2a-tTA;Ai93-222426	True
2	501474098	175	VISp	Cux2-CreERT2	Ai93(TITL-GCaMP6f)	102	511510736	three_session_C	222426	Cux2-CreERT2;Camk2a-tTA;Ai93-222426	True

	start_time	stop_time	orientation_in_degrees	spatial_frequency_in_cycles_per_degree	phase	is_blank_sweep
id
0	51.55381	51.80381	90.0	0.04	0.5	False
1	51.78653	52.03653	150.0	0.04	0.5	False
2	52.05250	52.30250	30.0	0.02	0.0	False

Data type	float32
Shape	(113723,)
Array size	444.23 KiB
Chunk shape	(113723,)
Compression	gzip
Compression opts	4
Uncompressed size (bytes)	454892
Compressed size (bytes)	421117
Compression ratio	1.0802033639107422

	global_roi_id	pixel_mask
id
517473350	517473350	[[55, 291, 1.0], [55, 292, 1.0], [56, 290, 1.0], [56, 291, 1.0], [56, 292, 1.0], [56, 293, 1.0], [56, 294, 1.0], [56, 295, 1.0], [56, 296, 1.0], [57, 289, 1.0], [57, 290, 1.0], [57, 291, 1.0], [57, 292, 1.0], [57, 293, 1.0], [57, 294, 1.0], [57, 295, 1.0], [57, 296, 1.0], [57, 297, 1.0], [57, 298, 1.0], [58, 288, 1.0], [58, 289, 1.0], [58, 290, 1.0], [58, 291, 1.0], [58, 292, 1.0], [58, 293, 1.0], [58, 294, 1.0], [58, 295, 1.0], [58, 296, 1.0], [58, 297, 1.0], [58, 298, 1.0], [58, 299, 1.0], [59, 288, 1.0], [59, 289, 1.0], [59, 290, 1.0], [59, 291, 1.0], [59, 292, 1.0], [59, 293, 1.0], [59, 294, 1.0], [59, 295, 1.0], [59, 296, 1.0], [59, 297, 1.0], [59, 298, 1.0], [59, 299, 1.0], [60, 288, 1.0], [60, 289, 1.0], [60, 290, 1.0], [60, 291, 1.0], [60, 292, 1.0], [60, 293, 1.0], [60, 294, 1.0], [60, 295, 1.0], [60, 296, 1.0], [60, 297, 1.0], [60, 298, 1.0], [60, 299, 1.0], [61, 287, 1.0], [61, 288, 1.0], [61, 289, 1.0], [61, 290, 1.0], [61, 291, 1.0], [61, 292, 1.0], [61, 293, 1.0], [61, 294, 1.0], [61, 295, 1.0], [61, 296, 1.0], [61, 297, 1.0], [61, 298, 1.0], [61, 299, 1.0], [62, 287, 1.0], [62, 288, 1.0], [62, 289, 1.0], [62, 290, 1.0], [62, 291, 1.0], [62, 292, 1.0], [62, 293, 1.0], [62, 294, 1.0], [62, 295, 1.0], [62, 296, 1.0], [62, 297, 1.0], [62, 298, 1.0], [62, 299, 1.0], [63, 287, 1.0], [63, 288, 1.0], [63, 289, 1.0], [63, 290, 1.0], [63, 291, 1.0], [63, 292, 1.0], [63, 293, 1.0], [63, 294, 1.0], [63, 295, 1.0], [63, 296, 1.0], [63, 297, 1.0], [63, 298, 1.0], [63, 299, 1.0], [64, 287, 1.0], [64, 288, 1.0], [64, 289, 1.0], [64, 290, 1.0], [64, 291, 1.0], [64, 292, 1.0], ...]
517473341	517473341	[[58, 272, 1.0], [58, 273, 1.0], [58, 274, 1.0], [58, 275, 1.0], [58, 276, 1.0], [58, 277, 1.0], [59, 271, 1.0], [59, 272, 1.0], [59, 273, 1.0], [59, 274, 1.0], [59, 275, 1.0], [59, 276, 1.0], [59, 277, 1.0], [59, 278, 1.0], [59, 279, 1.0], [59, 280, 1.0], [59, 281, 1.0], [59, 282, 1.0], [60, 271, 1.0], [60, 272, 1.0], [60, 273, 1.0], [60, 274, 1.0], [60, 275, 1.0], [60, 276, 1.0], [60, 277, 1.0], [60, 278, 1.0], [60, 279, 1.0], [60, 280, 1.0], [60, 281, 1.0], [60, 282, 1.0], [60, 283, 1.0], [61, 270, 1.0], [61, 271, 1.0], [61, 272, 1.0], [61, 273, 1.0], [61, 274, 1.0], [61, 275, 1.0], [61, 276, 1.0], [61, 277, 1.0], [61, 278, 1.0], [61, 279, 1.0], [61, 280, 1.0], [61, 281, 1.0], [61, 282, 1.0], [62, 270, 1.0], [62, 271, 1.0], [62, 272, 1.0], [62, 273, 1.0], [62, 274, 1.0], [62, 275, 1.0], [62, 276, 1.0], [62, 277, 1.0], [62, 278, 1.0], [62, 279, 1.0], [62, 280, 1.0], [62, 281, 1.0], [62, 282, 1.0], [63, 269, 1.0], [63, 270, 1.0], [63, 271, 1.0], [63, 272, 1.0], [63, 273, 1.0], [63, 274, 1.0], [63, 275, 1.0], [63, 276, 1.0], [63, 277, 1.0], [63, 278, 1.0], [63, 279, 1.0], [63, 280, 1.0], [63, 281, 1.0], [63, 282, 1.0], [64, 268, 1.0], [64, 269, 1.0], [64, 270, 1.0], [64, 271, 1.0], [64, 272, 1.0], [64, 273, 1.0], [64, 274, 1.0], [64, 275, 1.0], [64, 276, 1.0], [64, 277, 1.0], [64, 278, 1.0], [64, 279, 1.0], [64, 280, 1.0], [64, 281, 1.0], [64, 282, 1.0], [65, 268, 1.0], [65, 269, 1.0], [65, 270, 1.0], [65, 271, 1.0], [65, 272, 1.0], [65, 273, 1.0], [65, 274, 1.0], [65, 275, 1.0], [65, 276, 1.0], [65, 277, 1.0], [65, 278, 1.0], [65, 279, 1.0], [65, 280, 1.0], [65, 281, 1.0], ...]
517473313	517473313	[[33, 459, 1.0], [33, 460, 1.0], [33, 461, 1.0], [33, 462, 1.0], [33, 463, 1.0], [34, 458, 1.0], [34, 459, 1.0], [34, 460, 1.0], [34, 461, 1.0], [34, 462, 1.0], [34, 463, 1.0], [34, 464, 1.0], [35, 457, 1.0], [35, 458, 1.0], [35, 459, 1.0], [35, 460, 1.0], [35, 461, 1.0], [35, 462, 1.0], [35, 463, 1.0], [35, 464, 1.0], [35, 465, 1.0], [35, 466, 1.0], [36, 456, 1.0], [36, 457, 1.0], [36, 458, 1.0], [36, 459, 1.0], [36, 460, 1.0], [36, 461, 1.0], [36, 462, 1.0], [36, 463, 1.0], [36, 464, 1.0], [36, 465, 1.0], [36, 466, 1.0], [36, 467, 1.0], [37, 456, 1.0], [37, 457, 1.0], [37, 458, 1.0], [37, 459, 1.0], [37, 460, 1.0], [37, 461, 1.0], [37, 462, 1.0], [37, 463, 1.0], [37, 464, 1.0], [37, 465, 1.0], [37, 466, 1.0], [37, 467, 1.0], [38, 455, 1.0], [38, 456, 1.0], [38, 457, 1.0], [38, 458, 1.0], [38, 459, 1.0], [38, 460, 1.0], [38, 461, 1.0], [38, 462, 1.0], [38, 463, 1.0], [38, 464, 1.0], [38, 465, 1.0], [38, 466, 1.0], [38, 467, 1.0], [39, 455, 1.0], [39, 456, 1.0], [39, 457, 1.0], [39, 458, 1.0], [39, 459, 1.0], [39, 460, 1.0], [39, 461, 1.0], [39, 462, 1.0], [39, 463, 1.0], [39, 464, 1.0], [39, 465, 1.0], [39, 466, 1.0], [39, 467, 1.0], [40, 455, 1.0], [40, 456, 1.0], [40, 457, 1.0], [40, 458, 1.0], [40, 459, 1.0], [40, 460, 1.0], [40, 461, 1.0], [40, 462, 1.0], [40, 463, 1.0], [40, 464, 1.0], [40, 465, 1.0], [40, 466, 1.0], [41, 455, 1.0], [41, 456, 1.0], [41, 457, 1.0], [41, 458, 1.0], [41, 459, 1.0], [41, 460, 1.0], [41, 461, 1.0], [41, 462, 1.0], [41, 463, 1.0], [41, 464, 1.0], [41, 465, 1.0], [41, 466, 1.0], [42, 455, 1.0], [42, 456, 1.0], [42, 457, 1.0], [42, 458, 1.0], ...]
517473255	517473255	[[41, 493, 1.0], [41, 494, 1.0], [41, 495, 1.0], [41, 496, 1.0], [42, 492, 1.0], [42, 493, 1.0], [42, 494, 1.0], [42, 495, 1.0], [42, 496, 1.0], [42, 497, 1.0], [43, 489, 1.0], [43, 490, 1.0], [43, 491, 1.0], [43, 492, 1.0], [43, 493, 1.0], [43, 494, 1.0], [43, 495, 1.0], [43, 496, 1.0], [43, 497, 1.0], [43, 498, 1.0], [44, 488, 1.0], [44, 489, 1.0], [44, 490, 1.0], [44, 491, 1.0], [44, 492, 1.0], [44, 493, 1.0], [44, 494, 1.0], [44, 495, 1.0], [44, 496, 1.0], [44, 497, 1.0], [44, 498, 1.0], [45, 487, 1.0], [45, 488, 1.0], [45, 489, 1.0], [45, 490, 1.0], [45, 491, 1.0], [45, 492, 1.0], [45, 493, 1.0], [45, 494, 1.0], [45, 495, 1.0], [45, 496, 1.0], [45, 497, 1.0], [45, 498, 1.0], [46, 486, 1.0], [46, 487, 1.0], [46, 488, 1.0], [46, 489, 1.0], [46, 490, 1.0], [46, 491, 1.0], [46, 492, 1.0], [46, 493, 1.0], [46, 494, 1.0], [46, 495, 1.0], [46, 496, 1.0], [46, 497, 1.0], [46, 498, 1.0], [47, 486, 1.0], [47, 487, 1.0], [47, 488, 1.0], [47, 489, 1.0], [47, 490, 1.0], [47, 491, 1.0], [47, 492, 1.0], [47, 493, 1.0], [47, 494, 1.0], [47, 495, 1.0], [47, 496, 1.0], [47, 497, 1.0], [47, 498, 1.0], [47, 499, 1.0], [48, 486, 1.0], [48, 487, 1.0], [48, 488, 1.0], [48, 489, 1.0], [48, 490, 1.0], [48, 491, 1.0], [48, 492, 1.0], [48, 493, 1.0], [48, 494, 1.0], [48, 495, 1.0], [48, 496, 1.0], [48, 497, 1.0], [48, 498, 1.0], [48, 499, 1.0], [49, 486, 1.0], [49, 487, 1.0], [49, 488, 1.0], [49, 489, 1.0], [49, 490, 1.0], [49, 491, 1.0], [49, 492, 1.0], [49, 493, 1.0], [49, 494, 1.0], [49, 495, 1.0], [49, 496, 1.0], [49, 497, 1.0], [49, 498, 1.0], [49, 499, 1.0], [50, 486, 1.0], [50, 487, 1.0], ...]

	start_time	stop_time	stimulus_type
id
0	51.55381	531.93378	static_gratings
1	561.98893	1042.44259	natural_scenes
2	1047.42969	1342.66703	spontaneous
3	1342.70028	1823.72944	natural_scenes

Data type	uint16
Shape	(113888, 512, 512)
Array size	55.61 GiB
Chunk shape	(19, 512, 512)
Compression	gzip
Compression opts	4
Uncompressed size (bytes)	59710111744
Compressed size (bytes)	40604349646
Compression ratio	1.4705348630028394

	start_time	stop_time	orientation_in_degrees	spatial_frequency_in_cycles_per_degree	temporal_frequency_in_hz	is_blank_sweep
id
0	34.35150	36.35150	135.0	0.04	1.0	False
1	37.33713	39.33713	0.0	0.04	2.0	False
2	40.35595	42.35595	0.0	0.04	8.0	False
3	43.34159	45.34159	225.0	0.04	8.0	False