.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/domain/ophys.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_tutorials_domain_ophys.py: Calcium Imaging Data ==================== This tutorial will demonstrate how to write calcium imaging data. The workflow demonstrated here involves five main steps: 1. Create imaging plane 2. Add acquired two-photon images 3. Add motion correction (optional) 4. Add image segmentation 5. Add fluorescence and dF/F responses It is recommended to cover :ref:`basics` before this tutorial. .. note:: It is recommended to check if your source data is supported by `NeuroConv Optical Physiology Gallery `_. If it is supported, it is recommended to use NeuroConv to convert your data. The following examples will reference variables that may not be defined within the block they are used in. For clarity, we define them here: .. GENERATED FROM PYTHON SOURCE LINES 24-45 .. code-block:: Python from datetime import datetime from uuid import uuid4 import matplotlib.pyplot as plt import numpy as np from dateutil.tz import tzlocal from pynwb import NWBHDF5IO, NWBFile, TimeSeries from pynwb.image import ImageSeries from pynwb.ophys import ( CorrectedImageStack, Fluorescence, ImageSegmentation, MotionCorrection, OnePhotonSeries, OpticalChannel, RoiResponseSeries, TwoPhotonSeries, ) .. GENERATED FROM PYTHON SOURCE LINES 47-51 Creating the NWB file --------------------- When creating an NWB file, the first step is to create the :py:class:`~pynwb.file.NWBFile` object. .. GENERATED FROM PYTHON SOURCE LINES 51-66 .. code-block:: Python nwbfile = NWBFile( session_description="my first synthetic recording", identifier=str(uuid4()), session_start_time=datetime.now(tzlocal()), experimenter=[ "Baggins, Bilbo", ], lab="Bag End Laboratory", institution="University of Middle Earth at the Shire", experiment_description="I went on an adventure to reclaim vast treasures.", keywords=["ecephys", "exploration", "wanderlust"], related_publications="doi:10.1016/j.neuron.2016.12.011", ) .. GENERATED FROM PYTHON SOURCE LINES 67-95 Imaging Plane ------------- First, we must create an :py:class:`~pynwb.ophys.ImagingPlane` object, which will hold information about the area and method used to collect the optical imaging data. This first requires creation of a :py:class:`~pynwb.device.Device` object for the microscope and an :py:class:`~pynwb.ophys.OpticalChannel` object. .. only:: html .. image:: ../../_static/ImagingPlane.svg :width: 50% :alt: imaging plane UML diagram :align: center .. only:: latex .. image:: ../../_static/ImagingPlane.png :width: 50% :alt: imaging plane UML diagram :align: center Create a :py:class:`~pynwb.device.Device` named ``"Thorlabs Bergamo II"`` in the :py:class:`~pynwb.file.NWBFile` object. The fields ``description``, ``serial_number``, and ``model`` are optional, but recommended. The :py:class:`~pynwb.device.DeviceModel` object stores information about the device model, which can be useful when searching a set of NWB files or a data archive for all files that use a specific device model (e.g., specific microscope model). Then create an :py:class:`~pynwb.ophys.OpticalChannel` named ``"OpticalChannel"``. .. GENERATED FROM PYTHON SOURCE LINES 95-116 .. code-block:: Python device_model = nwbfile.create_device_model( name="Thorlabs Bergamo II Model", description="Two-photon microscope for in vivo imaging", manufacturer="Thorlabs", model_number="Bergamo II", ) device = nwbfile.create_device( name="Thorlabs Bergamo II", description="Two-photon microscope for in vivo imaging", model=device_model, serial_number="SN-123456789", ) optical_channel = OpticalChannel( name="OpticalChannel", description="an optical channel", emission_lambda=500.0, ) .. GENERATED FROM PYTHON SOURCE LINES 117-119 Now, create a :py:class:`~pynwb.ophys.ImagingPlane` named ``"ImagingPlane"``, passing in the :py:class:`~pynwb.ophys.OpticalChannel` object and the :py:class:`~pynwb.device.Device` object. .. GENERATED FROM PYTHON SOURCE LINES 119-135 .. code-block:: Python imaging_plane = nwbfile.create_imaging_plane( name="ImagingPlane", optical_channel=optical_channel, imaging_rate=30.0, description="a very interesting part of the brain", device=device, excitation_lambda=600.0, indicator="GFP", location="V1", grid_spacing=[0.01, 0.01], grid_spacing_unit="meters", origin_coords=[1.0, 2.0, 3.0], origin_coords_unit="meters", ) .. GENERATED FROM PYTHON SOURCE LINES 136-140 One-photon Series ----------------- Now that we have our :py:class:`~pynwb.ophys.ImagingPlane`, we can create a :py:class:`~pynwb.ophys.OnePhotonSeries` object to store raw one-photon imaging data. .. GENERATED FROM PYTHON SOURCE LINES 140-151 .. code-block:: Python # the image data will be stored inside the NWB file one_p_series = OnePhotonSeries( name="OnePhotonSeries", description="Raw 1p data", data=np.ones((1000, 100, 100)), imaging_plane=imaging_plane, rate=1.0, unit="normalized amplitude", ) .. GENERATED FROM PYTHON SOURCE LINES 152-155 Since these one-photon data are acquired data, we will add the :py:class:`~pynwb.ophys.OnePhotonSeries` objects to the :py:class:`~pynwb.file.NWBFile` as acquired data. .. GENERATED FROM PYTHON SOURCE LINES 155-158 .. code-block:: Python nwbfile.add_acquisition(one_p_series) .. GENERATED FROM PYTHON SOURCE LINES 159-178 Two-photon Series ----------------- :py:class:`~pynwb.ophys.TwoPhotonSeries` objects store acquired two-photon imaging data. This class behaves similarly to :py:class:`~pynwb.ophys.OnePhotonSeries`. .. only:: html .. image:: ../../_static/TwoPhotonSeries.svg :width: 80% :alt: two-photon series UML diagram :align: center .. only:: latex .. image:: ../../_static/TwoPhotonSeries.png :width: 80% :alt: two-photon series UML diagram :align: center .. GENERATED FROM PYTHON SOURCE LINES 178-191 .. code-block:: Python # the image data will be stored inside the NWB file two_p_series = TwoPhotonSeries( name="TwoPhotonSeries", description="Raw 2p data", data=np.ones((1000, 100, 100)), imaging_plane=imaging_plane, rate=1.0, unit="normalized amplitude", ) nwbfile.add_acquisition(two_p_series) .. GENERATED FROM PYTHON SOURCE LINES 192-198 Motion Correction (optional) ---------------------------- You can also store the result of motion correction using a :py:class:`~pynwb.ophys.MotionCorrection` object, which is a :py:class:`~pynwb.core.MultiContainerInterface` (similar to :py:class:`~pynwb.behavior.Position`) which holds 1 or more :py:class:`~pynwb.ophys.CorrectedImageStack` objects. .. GENERATED FROM PYTHON SOURCE LINES 198-226 .. code-block:: Python corrected = ImageSeries( name="corrected", # this must be named "corrected" description="A motion corrected image stack", data=np.ones((1000, 100, 100)), unit="na", format="raw", starting_time=0.0, rate=1.0, ) xy_translation = TimeSeries( name="xy_translation", description="x,y translation in pixels", data=np.ones((1000, 2)), unit="pixels", starting_time=0.0, rate=1.0, ) corrected_image_stack = CorrectedImageStack( corrected=corrected, original=one_p_series, xy_translation=xy_translation, ) motion_correction = MotionCorrection(corrected_image_stacks=[corrected_image_stack]) .. GENERATED FROM PYTHON SOURCE LINES 227-230 We will create a :py:class:`~pynwb.base.ProcessingModule` named "ophys" to store optical physiology data and add the motion correction data to the :py:class:`~pynwb.file.NWBFile`. .. GENERATED FROM PYTHON SOURCE LINES 230-237 .. code-block:: Python ophys_module = nwbfile.create_processing_module( name="ophys", description="optical physiology processed data" ) ophys_module.add(motion_correction) .. GENERATED FROM PYTHON SOURCE LINES 238-284 Plane Segmentation ------------------ The :py:class:`~pynwb.ophys.PlaneSegmentation` class stores the detected regions of interest in the :py:class:`~pynwb.ophys.TwoPhotonSeries` data. :py:class:`~pynwb.ophys.PlaneSegmentation` is a subclass of :py:class:`~hdmf.common.table.DynamicTable`, where each row represents a single region of interest (ROI). .. only:: html .. image:: ../../_static/PlaneSegmentation.svg :width: 60% :alt: plane segmentation UML diagram :align: center .. only:: latex .. image:: ../../_static/PlaneSegmentation.png :width: 60% :alt: plane segmentation UML diagram :align: center The :py:class:`~pynwb.ophys.ImageSegmentation` class can contain multiple :py:class:`~pynwb.ophys.PlaneSegmentation` tables, so that we can store results of different segmentation algorithms or different segmentation classes. .. only:: html .. image:: ../../_static/ImageSegmentation.svg :width: 100% :alt: image segmentation UML diagram :align: center .. only:: latex .. image:: ../../_static/ImageSegmentation.png :width: 100% :alt: image segmentation UML diagram :align: center First, we create an :py:class:`~pynwb.ophys.ImageSegmentation` object, then from that object we create a :py:class:`~pynwb.ophys.PlaneSegmentation` table with a link to the :py:class:`~pynwb.ophys.ImagingPlane` created earlier. Then we will add the :py:class:`~pynwb.ophys.ImageSegmentation` object to the previously created :py:class:`~pynwb.base.ProcessingModule`. .. GENERATED FROM PYTHON SOURCE LINES 284-296 .. code-block:: Python img_seg = ImageSegmentation() ps = img_seg.create_plane_segmentation( name="PlaneSegmentation", description="output from segmenting my favorite imaging plane", imaging_plane=imaging_plane, reference_images=one_p_series, # optional ) ophys_module.add(img_seg) .. GENERATED FROM PYTHON SOURCE LINES 297-308 Regions Of Interest (ROIs) -------------------------- Image masks ^^^^^^^^^^^ You can add ROIs to the :py:class:`~pynwb.ophys.PlaneSegmentation` table using an image mask or a pixel mask. An image mask is an array that is the same size as a single frame of the :py:class:`~pynwb.ophys.TwoPhotonSeries` that indicates the mask weight of each pixel in the image. Image mask values (weights) may be boolean or continuous between 0 and 1. .. GENERATED FROM PYTHON SOURCE LINES 308-324 .. code-block:: Python for _ in range(30): image_mask = np.zeros((100, 100)) # randomly generate example image masks x = np.random.randint(0, 95) y = np.random.randint(0, 95) image_mask[x:x + 5, y:y + 5] = 1 # add image mask to plane segmentation ps.add_roi(image_mask=image_mask) # show one of the image masks plt.imshow(image_mask) .. GENERATED FROM PYTHON SOURCE LINES 325-331 Pixel masks ^^^^^^^^^^^ Alternatively, you could define ROIs using a pixel mask, which is an array of triplets (x, y, weight) that have a non-zero weight. All undefined pixels are assumed to be 0. .. GENERATED FROM PYTHON SOURCE LINES 331-354 .. code-block:: Python ps2 = img_seg.create_plane_segmentation( name="PlaneSegmentation2", description="output from segmenting my favorite imaging plane", imaging_plane=imaging_plane, reference_images=one_p_series, # optional ) for _ in range(30): # randomly generate example starting points for region x = np.random.randint(0, 95) y = np.random.randint(0, 95) # define an example 4 x 3 region of pixels of weight '1' pixel_mask = [] for ix in range(x, x + 4): for iy in range(y, y + 3): pixel_mask.append((ix, iy, 1)) # add pixel mask to plane segmentation ps2.add_roi(pixel_mask=pixel_mask) .. GENERATED FROM PYTHON SOURCE LINES 355-366 Voxel masks ^^^^^^^^^^^ When storing the segmentation of volumetric imaging, you can use 3D imaging masks. Alternatively, you could define ROIs using a voxel mask, which is an array of triplets (x, y, z, weight) that have a non-zero weight. All undefined voxels are assumed to be 0. .. note:: You need to be consistent within a :py:class:`~pynwb.ophys.PlaneSegmentation` table. You can add ROIs either using image masks, pixel masks, or voxel masks. .. GENERATED FROM PYTHON SOURCE LINES 366-395 .. code-block:: Python ps3 = img_seg.create_plane_segmentation( name="PlaneSegmentation3", description="output from segmenting my favorite imaging plane", imaging_plane=imaging_plane, reference_images=one_p_series, # optional ) from itertools import product for _ in range(30): # randomly generate example starting points for region x = np.random.randint(0, 95) y = np.random.randint(0, 95) z = np.random.randint(0, 15) # define an example 4 x 3 x 2 voxel region of weight '0.5' voxel_mask = [] for ix, iy, iz in product( range(x, x + 4), range(y, y + 3), range(z, z + 2) ): voxel_mask.append((ix, iy, iz, 0.5)) # add voxel mask to plane segmentation ps3.add_roi(voxel_mask=voxel_mask) .. GENERATED FROM PYTHON SOURCE LINES 396-398 We can view the :py:class:`~pynwb.ophys.PlaneSegmentation` table with pixel masks in tabular form by converting it to a :py:class:`~pandas.DataFrame`. .. GENERATED FROM PYTHON SOURCE LINES 398-401 .. code-block:: Python ps2.to_dataframe() .. GENERATED FROM PYTHON SOURCE LINES 402-433 Storing Fluorescence Measurements --------------------------------- Now that the regions of interest are stored, you can store fluorescence data for these ROIs. This type of data is stored using the :py:class:`~pynwb.ophys.RoiResponseSeries` classes. .. only:: html .. image:: ../../_static/RoiResponseSeries.svg :width: 80% :alt: RoiResponseSeries UML diagram :align: center .. only:: latex .. image:: ../../_static/RoiResponseSeries.png :width: 80% :alt: RoiResponseSeries UML diagram :align: center To create a :py:class:`~pynwb.ophys.RoiResponseSeries` object, we will need to reference a set of rows from a :py:class:`~pynwb.ophys.PlaneSegmentation` table to indicate which ROIs correspond to which rows of your recorded data matrix. This is done using a :py:class:`~hdmf.common.table.DynamicTableRegion`, which is a type of link that allows you to reference specific rows of a :py:class:`~hdmf.common.table.DynamicTable`, such as a :py:class:`~pynwb.ophys.PlaneSegmentation` table by row indices. First, we create a :py:class:`~hdmf.common.table.DynamicTableRegion` that references the first two ROIs of the :py:class:`~pynwb.ophys.PlaneSegmentation` table. .. GENERATED FROM PYTHON SOURCE LINES 433-438 .. code-block:: Python rt_region = ps.create_roi_table_region( region=[0, 1], description="the first of two ROIs" ) .. GENERATED FROM PYTHON SOURCE LINES 439-441 Then we create a :py:class:`~pynwb.ophys.RoiResponseSeries` object to store fluorescence data for those two ROIs. .. GENERATED FROM PYTHON SOURCE LINES 441-451 .. code-block:: Python roi_resp_series = RoiResponseSeries( name="RoiResponseSeries", description="Fluorescence responses for two ROIs", data=np.ones((50, 2)), # 50 samples, 2 ROIs rois=rt_region, unit="lumens", rate=30.0, ) .. GENERATED FROM PYTHON SOURCE LINES 452-472 To help data analysis and visualization tools know that this :py:class:`~pynwb.ophys.RoiResponseSeries` object represents fluorescence data, we will store the :py:class:`~pynwb.ophys.RoiResponseSeries` object inside of a :py:class:`~pynwb.ophys.Fluorescence` object. Then add the :py:class:`~pynwb.ophys.Fluorescence` object into the same :py:class:`~pynwb.base.ProcessingModule` named ``"ophys"`` that we created earlier. .. only:: html .. image:: ../../_static/Fluorescence.svg :width: 100% :alt: fluorescence UML diagram :align: center .. only:: latex .. image:: ../../_static/Fluorescence.png :width: 100% :alt: fluorescence UML diagram :align: center .. GENERATED FROM PYTHON SOURCE LINES 472-476 .. code-block:: Python fl = Fluorescence(roi_response_series=roi_resp_series) ophys_module.add(fl) .. GENERATED FROM PYTHON SOURCE LINES 477-483 .. tip:: If you want to store dF/F data instead of fluorescence data, then store the :py:class:`~pynwb.ophys.RoiResponseSeries` object in a :py:class:`~pynwb.ophys.DfOverF` object, which works the same way as the :py:class:`~pynwb.ophys.Fluorescence` class. .. GENERATED FROM PYTHON SOURCE LINES 485-490 Write the file --------------------------------- Once we have finished adding all of our data to our :py:class:`~pynwb.file.NWBFile`, make sure to write the file. IO operations are carried out using :py:class:`~pynwb.NWBHDF5IO`. .. GENERATED FROM PYTHON SOURCE LINES 490-494 .. code-block:: Python with NWBHDF5IO("ophys_tutorial.nwb", "w") as io: io.write(nwbfile) .. GENERATED FROM PYTHON SOURCE LINES 495-511 Read the NWBFile ---------------- We can access the raw data by indexing ``nwbfile.acquisition`` with a name of the :py:class:`~pynwb.ophys.TwoPhotonSeries`, e.g., ``"TwoPhotonSeries1"``. We can also access the fluorescence responses by indexing ``nwbfile.processing`` with the name of the processing module, ``"ophys"``. Then, we can access the :py:class:`~pynwb.ophys.Fluorescence` object inside of the ``"ophys"`` processing module by indexing it with the name of the :py:class:`~pynwb.ophys.Fluorescence` object. The default name of :py:class:`~pynwb.ophys.Fluorescence` objects is ``"Fluorescence"``. Finally, we can access the :py:class:`~pynwb.ophys.RoiResponseSeries` object inside the :py:class:`~pynwb.ophys.Fluorescence` object by indexing it with the name of the :py:class:`~pynwb.ophys.RoiResponseSeries` object, which we named ``"RoiResponseSeries"``. .. GENERATED FROM PYTHON SOURCE LINES 511-519 .. code-block:: Python with NWBHDF5IO("ophys_tutorial.nwb", "r") as io: read_nwbfile = io.read() print(read_nwbfile.acquisition["TwoPhotonSeries"]) print(read_nwbfile.processing["ophys"]) print(read_nwbfile.processing["ophys"]["Fluorescence"]) print(read_nwbfile.processing["ophys"]["Fluorescence"]["RoiResponseSeries"]) .. GENERATED FROM PYTHON SOURCE LINES 520-530 Accessing your data ------------------------------ Data arrays are read passively from the file. Calling the data attribute on a :py:class:`~pynwb.base.TimeSeries` such as a :py:class:`~pynwb.ophys.RoiResponseSeries` does not read the data values, but presents an ``h5py`` object that can be indexed to read data. You can use the ``[:]`` operator to read the entire data array into memory. Load and print all the data values of the :py:class:`~pynwb.ophys.RoiResponseSeries` object representing the fluorescence data. .. GENERATED FROM PYTHON SOURCE LINES 530-537 .. code-block:: Python with NWBHDF5IO("ophys_tutorial.nwb", "r") as io: read_nwbfile = io.read() print(read_nwbfile.acquisition["TwoPhotonSeries"]) print(read_nwbfile.processing["ophys"]["Fluorescence"]["RoiResponseSeries"].data[:]) .. GENERATED FROM PYTHON SOURCE LINES 538-547 Accessing data regions ------------------------------ It is often preferable to read only a portion of the data. To do this, index or slice into the data attribute just like if you were indexing or slicing a :py:mod:`numpy` array. The following code prints elements ``0:10`` in the first dimension (time) and ``0:3`` (ROIs) in the second dimension from the fluorescence data we have written. .. GENERATED FROM PYTHON SOURCE LINES 547-558 .. code-block:: Python with NWBHDF5IO("ophys_tutorial.nwb", "r") as io: read_nwbfile = io.read() roi_resp_series = read_nwbfile.processing["ophys"]["Fluorescence"][ "RoiResponseSeries" ] print("section of fluorescence responses:") print(roi_resp_series.data[0:10, 0:3]) .. _sphx_glr_download_tutorials_domain_ophys.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ophys.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: ophys.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: ophys.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_