notebooks/generate_test_nwb.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate test data for plotting functionalities based on existing dandisets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import modules and define functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import modules\n",
    "import hdmf._version\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import platform\n",
    "import pynwb\n",
    "import seaborn as sns\n",
    "import uuid\n",
    "\n",
    "from datetime import datetime\n",
    "from dateutil.tz import tzlocal\n",
    "from nwbwidgets import nwb2widget\n",
    "from pynwb import NWBHDF5IO, NWBFile, TimeSeries\n",
    "from pynwb.ophys import TwoPhotonSeries, OpticalChannel, ImageSegmentation, Fluorescence\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define functions\n",
    "def load_nwb_file(filename):\n",
    "    directory = '..//test_data//'\n",
    "    read_io = NWBHDF5IO((directory + filename), 'r')\n",
    "    nwbfile_in = read_io.read()\n",
    "\n",
    "    return nwbfile_in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate initial NWB file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate the main nwb file\n",
    "start_time = datetime.now(tz=tzlocal())\n",
    "create_date = datetime.now(tz=tzlocal())\n",
    "\n",
    "unique_identifier = str(uuid.uuid4())\n",
    "\n",
    "hdmf_ver = 'v%s'%hdmf._version.get_versions()['version']\n",
    "info = 'An example NWB file created with pynwb v%s (hdmf %s), Python v%s'%(pynwb.__version__,hdmf_ver,platform.python_version())\n",
    "\n",
    "nwbfile_out = pynwb.NWBFile(session_description='Example with various datatypes for testing plotting functionality',\n",
    "                        identifier=unique_identifier,\n",
    "                        session_start_time=start_time,\n",
    "                        file_create_date=create_date,\n",
    "                        notes=info,\n",
    "                        experimenter='No one',\n",
    "                        experiment_description='We used a python script to synthesize this data.',\n",
    "                        institution='Institute X',\n",
    "                        lab='No Lab.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add example data modules\n",
    "### Multi-dimensional time series data + two photon image series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load test file \n",
    "time_series_multi_dim_fname = 'sub-661968859_ses-681698752_behavior+ophys.nwb'\n",
    "nwbfile_in = load_nwb_file(time_series_multi_dim_fname)\n",
    "proc_example = nwbfile_in.processing['brain_observatory_pipeline']\n",
    "\n",
    "#create device, imaging plan, optical channel used\n",
    "device = pynwb.file.Device('CAM2P_1')\n",
    "nwbfile_out.add_device(device)\n",
    "optical_channel = OpticalChannel('optchan', 'description', 480.)\n",
    "imaging_plane = nwbfile_out.create_imaging_plane('imaging_plane_1', optical_channel,\n",
    "                                             description='imaging_plane_1',\n",
    "                                             device=device, excitation_lambda=910.,\n",
    "                                             imaging_rate=30., indicator='GCaMP6f',\n",
    "                                             location='PFC', reference_frame='home')\n",
    "\n",
    "#make imaging series \n",
    "two_p_example = proc_example['max_project'].data[:,:,:]\n",
    "image_series = TwoPhotonSeries(name='image_series', \n",
    "                               dimension=[2], \n",
    "                               data = two_p_example,                                 # TO DO - add synthetic data here\n",
    "                               imaging_plane=imaging_plane,\n",
    "                               starting_frame=[0], \n",
    "                               starting_time=0.0, rate=1.0,\n",
    "                               unit='seconds')\n",
    "nwbfile_out.add_acquisition(image_series)\n",
    "\n",
    "\n",
    "# store image segmentation output\n",
    "mod = nwbfile_out.create_processing_module('ophys', 'contains optical physiology processed data')\n",
    "img_seg = ImageSegmentation()\n",
    "mod.add(img_seg)\n",
    "ps = img_seg.create_plane_segmentation('output from image plane', imaging_plane, 'my_planeseg', image_series)\n",
    "\n",
    "#add image mask\n",
    "mask_example = proc_example['ImageSegmentation']['PlaneSegmentation'].image_mask     # TO DO - add synthetic data here\n",
    "ps.add_roi(image_mask=mask_example[0])\n",
    "ps.add_roi(image_mask=mask_example[1])\n",
    "\n",
    "#add fluorescence and roi response - TO DO - add synthetic data here \n",
    "data_example = proc_example['Fluorescence']['DfOverF'].data[:,:2]                    # TO DO - add synthetic data here\n",
    "timestamps_example = proc_example['Fluorescence']['DfOverF'].timestamps[:]\n",
    "\n",
    "fl = Fluorescence()\n",
    "mod.add(fl)\n",
    "rt_region = ps.create_roi_table_region('the first of two ROIs', region=[0,1])\n",
    "rrs = fl.create_roi_response_series('DfOverF', data_example, rt_region, unit='lumens', timestamps=timestamps_example)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Behavioral events time series + trials table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# behavioral events time series\n",
    "behavior_events_markers_fname = 'sub-anm369962_ses-20170309.nwb'\n",
    "nwbfile_in = load_nwb_file(behavior_events_markers_fname)\n",
    "\n",
    "lick_data_example = nwbfile_in.acquisition['lick_times'].time_series['lick_left_times'].data[:] \n",
    "lick_timestamps_example = nwbfile_in.acquisition['lick_times'].time_series['lick_left_times'].timestamps[:]\n",
    "\n",
    "ts = TimeSeries(name='lick_left_times',\n",
    "                data=lick_data_example,                                              # TO DO - add synthetic data here\n",
    "                timestamps=lick_timestamps_example,\n",
    "                unit='a.u.')\n",
    "lick_times = pynwb.behavior.BehavioralEvents(name='lick_times', time_series=ts)\n",
    "\n",
    "nwbfile_out.add_acquisition(lick_times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add trials table\n",
    "trials_example = nwbfile_in.trials[:]\n",
    "\n",
    "for index, data in trials_example.iterrows():\n",
    "    nwbfile_out.add_trial(start_time=data.start_time, stop_time=data.stop_time)      # TO DO - add synthetic data here\n",
    "\n",
    "for col_name, col_data in trials_example.iteritems():\n",
    "    if col_name not in nwbfile_out.trials.colnames:\n",
    "        nwbfile_out.add_trial_column(name=col_name, description='',data=col_data.to_list())   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spatial series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# spatial series with x vs. y instead of data vs. timestamps\n",
    "N = 2001\n",
    "timestamps = np.linspace(0, N, N) #np.arange(N)\n",
    "ss1 = pynwb.behavior.SpatialSeries(name='spatial_series_1D',\n",
    "                                   reference_frame='Zero is origin..?',\n",
    "                                   data=np.cos(timestamps/4),\n",
    "                                   timestamps=timestamps,\n",
    "                                   comments='A wave...',\n",
    "                                   description='Description of this...')\n",
    "\n",
    "nwbfile_out.add_acquisition(ss1)\n",
    "\n",
    "x = np.cos(timestamps/4)\n",
    "y = np.sin(timestamps/5)\n",
    "ss2 = pynwb.behavior.SpatialSeries(name='spatial_series_2D',\n",
    "                                   reference_frame='Zero is origin..?',\n",
    "                                   data=np.array([x,y]).T,\n",
    "                                   timestamps=timestamps,\n",
    "                                   comments='A wave...',\n",
    "                                   description='Description of this...')\n",
    "\n",
    "nwbfile_out.add_acquisition(ss2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stimulus image series data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#optical series\n",
    "imaging_series_fname = 'sub-P10HMH_ses-20060901_ecephys+image.nwb'\n",
    "nwbfile_in = load_nwb_file(imaging_series_fname)\n",
    "\n",
    "optical_series_example = nwbfile_in.stimulus['StimulusPresentation'].data[:,:,:,:]   # TO DO - add synthetic data here\n",
    "timestamps_example=nwbfile_in.stimulus['StimulusPresentation'].timestamps[:]\n",
    "\n",
    "os = pynwb.image.OpticalSeries(name='StimulusPresentation',\n",
    "                               data=optical_series_example,\n",
    "                               timestamps=timestamps_example,\n",
    "                               distance=0.5,\n",
    "                               field_of_view=[0.2,0.3,0.5],\n",
    "                               orientation='lower left',\n",
    "                               unit='meters')\n",
    "nwbfile_out.add_stimulus(os)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Time stamps alignment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# time stamps with whole session vs. all aligned to one location\n",
    "time_stamps_aligned_fname = '10_03_19-1.nwb'\n",
    "nwbfile_in = load_nwb_file(time_stamps_aligned_fname)\n",
    "\n",
    "# add stimulus\n",
    "vcss_example_1 = nwbfile_in.stimulus['index_000'].data[:]                             #TO DO - add synthetic data here\n",
    "vcss_start_1 = nwbfile_in.stimulus['index_000'].starting_time\n",
    "vcss_example_2 = nwbfile_in.stimulus['index_001'].data[:]\n",
    "vcss_start_2 = nwbfile_in.stimulus['index_001'].starting_time\n",
    "\n",
    "# add acquisition\n",
    "vcs_example_1 = nwbfile_in.acquisition['index_000'].data[:]\n",
    "vcs_start_1 = nwbfile_in.acquisition['index_000'].starting_time\n",
    "vcs_example_2 = nwbfile_in.acquisition['index_001'].data[:]\n",
    "vcs_start_2 = nwbfile_in.acquisition['index_001'].starting_time\n",
    "\n",
    "# add devices and electrodes\n",
    "device = nwbfile_out.create_device(name='Heka EPC10')\n",
    "elec = nwbfile_out.create_icephys_electrode(name=\"elec0\",\n",
    "                                        description='an intracellular electrode',\n",
    "                                        device=device)\n",
    "\n",
    "# add vclamp stimulus series\n",
    "vcss1 = pynwb.icephys.VoltageClampStimulusSeries(name=\"vclamp_stim_series1\",\n",
    "                                        data=vcss_example_1,\n",
    "                                        starting_time=vcss_start_1,\n",
    "                                        rate=10000.0,\n",
    "                                        electrode=elec,\n",
    "                                        gain=1.)\n",
    "vcss2 = pynwb.icephys.VoltageClampStimulusSeries(name=\"vclamp_stim_series2\",\n",
    "                                        data=vcss_example_2,\n",
    "                                        starting_time=vcss_start_2,\n",
    "                                        rate=10000.0,\n",
    "                                        electrode=elec,\n",
    "                                        gain=1.)\n",
    "# add vclamp seriess\n",
    "vcs1 = pynwb.icephys.VoltageClampSeries(name=\"vclamp_series1\",\n",
    "                                        data=vcs_example_1,\n",
    "                                        starting_time=vcs_start_1,\n",
    "                                        rate=10000.0,\n",
    "                                        electrode=elec,\n",
    "                                        gain=100000000.0)\n",
    "vcs2 = pynwb.icephys.VoltageClampSeries(name=\"vclamp_series2\",\n",
    "                                        data=vcs_example_2,\n",
    "                                        starting_time=vcs_start_2,\n",
    "                                        rate=10000.0,\n",
    "                                        electrode=elec,\n",
    "                                        gain=100000000.0)\n",
    "\n",
    "nwbfile_out.add_stimulus(vcss1)\n",
    "nwbfile_out.add_stimulus(vcss2)\n",
    "nwbfile_out.add_acquisition(vcs1)\n",
    "nwbfile_out.add_acquisition(vcs2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Single units"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#single unit raster plots\n",
    "single_units_fname = 'sub-anm184389_ses-20130207_behavior+ecephys.nwb'\n",
    "nwbfile_in = load_nwb_file(single_units_fname)\n",
    "\n",
    "#add electrodes info\n",
    "device = nwbfile_out.create_device(name='Neuronexus_32chan')\n",
    "electrode_name = 'Neuronexus1'\n",
    "description = \"an example probe\"\n",
    "location = \"somewhere in the hippocampus\"\n",
    "electrode_group = nwbfile_out.create_electrode_group(electrode_name,\n",
    "                                                 description=description,\n",
    "                                                 location=location,\n",
    "                                                 device=device)\n",
    "\n",
    "electrode_example = nwbfile_in.electrodes[:]\n",
    "for idx, elec in electrode_example.iterrows():\n",
    "    nwbfile_out.add_electrode(id=idx,                                                #TO DO - add synthetic data here\n",
    "                          x=elec.x, y=elec.y, z=elec.z,\n",
    "                          imp=elec.imp,\n",
    "                          location='CA1', filtering='none',\n",
    "                          group=electrode_group)\n",
    "\n",
    "# add single units info\n",
    "units_example = nwbfile_in.units[:]\n",
    "for idx, unit in units_example.iterrows():\n",
    "    nwbfile_out.add_unit(id=idx,                                                     #TO DO - add synthetic data here\n",
    "                     spike_times=unit.spike_times,\n",
    "                     waveform_mean=unit.waveform_mean,\n",
    "                     waveform_sd = unit.waveform_sd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save the data file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/stephanieprince/opt/anaconda3/envs/gsoc-project/lib/python3.8/site-packages/hdmf/build/objectmapper.py:233: DtypeConversionWarning: Spec 'TimeSeries/timestamps': Value with data type float64 is being converted to data type float64 as specified.\n",
      "  warnings.warn(full_warning_msg, DtypeConversionWarning)\n",
      "/Users/stephanieprince/opt/anaconda3/envs/gsoc-project/lib/python3.8/site-packages/hdmf/build/objectmapper.py:233: DtypeConversionWarning: Spec 'RoiResponseSeries/timestamps': Value with data type float64 is being converted to data type float64 as specified.\n",
      "  warnings.warn(full_warning_msg, DtypeConversionWarning)\n",
      "/Users/stephanieprince/opt/anaconda3/envs/gsoc-project/lib/python3.8/site-packages/hdmf/build/objectmapper.py:233: DtypeConversionWarning: Spec 'OpticalSeries/timestamps': Value with data type float64 is being converted to data type float64 as specified.\n",
      "  warnings.warn(full_warning_msg, DtypeConversionWarning)\n",
      "/Users/stephanieprince/opt/anaconda3/envs/gsoc-project/lib/python3.8/site-packages/hdmf/build/objectmapper.py:233: DtypeConversionWarning: Spec 'Units/waveform_mean': Value with data type float32 is being converted to data type float32 as specified.\n",
      "  warnings.warn(full_warning_msg, DtypeConversionWarning)\n",
      "/Users/stephanieprince/opt/anaconda3/envs/gsoc-project/lib/python3.8/site-packages/hdmf/build/objectmapper.py:233: DtypeConversionWarning: Spec 'Units/waveform_sd': Value with data type float32 is being converted to data type float32 as specified.\n",
      "  warnings.warn(full_warning_msg, DtypeConversionWarning)\n"
     ]
    }
   ],
   "source": [
    "with NWBHDF5IO('..//test_data//test-plotting-datatypes.nwb', 'w') as export_io:\n",
    "    export_io.write(nwbfile_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inspect the output file with NWB widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "nwbfile_out = load_nwb_file('..//test_data//test_plotting_datatypes.nwb')\n",
    "nwb2widget(nwbfile_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}