{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 03 - Introduction to neural data types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Goals\n",
    "\n",
    "- Become aware of the processing steps typically applied to raw data acquired with a Neuralynx system\n",
    "- Learn the data types used for neural data analysis\n",
    "- Use the loading functions for all these files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Renaming\n",
    "The naming convention in our lab is:\n",
    "R###-YYYY-MM-DD-filetype.extension\n",
    "\n",
    "By default, Neuralynx will name files filetype.extension,\n",
    "so the ratID and date need to be added as a prefix.\n",
    "This can be done with a script,\n",
    "just be careful to double-check that your script only renames the files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Manual annotation\n",
    "Some common features that are manually annotated are:\n",
    "- experimenter name\n",
    "- subject ID\n",
    "- session ID\n",
    "- experiment time blocks\n",
    "\n",
    "These are stored for easy access in an `info.py` file within your project folder."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spike sorting\n",
    "Spike sorting is the process of automated and manual detection \n",
    "of individual units from a tetrode or probe recording.\n",
    "See Module08 for details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data files overview\n",
    "\n",
    "A typical recording session from a Neuralynx system contains:\n",
    "\n",
    "- ''.ncs'' (\"**N**euralynx **C**ontinuously **S**ampled\") files typically contain \n",
    "local field potentials (LFPs) sampled at 2kHz and filtered between 1 and 475 Hz.\n",
    "- ''.t' or ''.\\_t'' (\"**T**etrode\") files,\n",
    "which are generated by MClust from the raw ''.ntt'' (\"**N**euralynx **T**e**T**rode\") files\n",
    "contains a set of times from individual putative neurons.\n",
    "- ''.nvt'' (\"**N**euralynx **V**ideo **T**racking\") file\n",
    "contains the location of the rat as tracked by an overhead camera,\n",
    "typically sampled at 30 Hz.\n",
    "- ''.nev'' (\"**N**euralynx **EV**ents\") file \n",
    "contains timestamps and labels of events\n",
    "- ''info.py'' file contains experimenter-provided information that describes this data set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nept data types\n",
    "\n",
    "Documentation for all the data types we use in Nept can be found \n",
    "[here](http://nept.readthedocs.io/en/latest/core.html).\n",
    "But we'll also briefly describe them in this module."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### AnalogSignal\n",
    "Genrally, data acquisition systems work by _sampling_ some quantity of interest,\n",
    "so often only periodically measurements are taken.\n",
    "A sampled signal is essentially a list of values (data points), taken at specific times. \n",
    "Thus, what we need to fully describe such a signal is two arrays of the same length: \n",
    "one with the timestamps and the other with the corresponding data.\n",
    "For this, we use an `AnalogSignal`, which has `time` and `data` fields."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### LocalFieldPotential\n",
    "\n",
    "A `LocalFieldPotential` is a subclass of `AnalogSignal`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Position\n",
    "\n",
    "A `Position` is a subclass of `AnalogSignal`,\n",
    "with some unique methods that are useful when dealing with positions.\n",
    "This includes computing the `distance` between two positions, \n",
    "or computing the `speed` of the position."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Epoch\n",
    "\n",
    "Epochs describe occurrences that have start and stop times, \n",
    "such as a _trial_ of an experiment or the presence of a cue (e.g. a light or a tone)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SpikeTrain\n",
    "\n",
    "Much of neuroscience analyses attributes particular significance to action potentials, or \"spikes\", \n",
    "which are typically understood as all-or-none events that occur at a specific point in time. \n",
    "It is not necessary to state all the times at which there was no spike to describe a train of spikes;\n",
    "it suffices to maintain a list of those times (known as \"timestamps\") \n",
    "at which a spike was emitted."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading real data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T17:58:42.878857Z",
     "start_time": "2017-07-26T17:58:42.868331Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# import necessary packages\n",
    "%matplotlib inline\n",
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import nept\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Define where your data folder is located\n",
    "data_path = os.path.join(os.path.abspath('.'), 'data')\n",
    "data_folder = os.path.join(data_path, 'R016-2012-10-03')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T17:58:43.250323Z",
     "start_time": "2017-07-26T17:58:43.243801Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Load the info file, which contains experiment-specific information\n",
    "sys.path.append(data_folder)\n",
    "import r016d3 as info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T17:58:44.183824Z",
     "start_time": "2017-07-26T17:58:43.836078Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Load LFP (.ncs) from rat ventral striatum (same as Module01)\n",
    "lfp = nept.load_lfp(os.path.join(data_folder, info.lfp_theta_filename))\n",
    "\n",
    "# Slice the LFP to our time of interest\n",
    "start = 1261\n",
    "stop = 1262\n",
    "sliced_lfp = lfp.time_slice(start, stop)\n",
    "\n",
    "# Plot the sliced LFP\n",
    "plt.plot(sliced_lfp.time, sliced_lfp.data)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T17:59:31.471313Z",
     "start_time": "2017-07-26T17:59:31.116976Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Load position (.nvt) from this experiment\n",
    "position = nept.load_position(os.path.join(data_folder, info.position_filename), \n",
    "                              info.pxl_to_cm)\n",
    "\n",
    "# Plot the position\n",
    "plt.plot(position.x, position.y, 'k.', ms=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:03:11.162172Z",
     "start_time": "2017-07-26T18:03:11.147161Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Load events (.nev) from this experiment\n",
    "events = nept.load_events(os.path.join(data_folder, info.event_filename), \n",
    "                          info.event_labels)\n",
    "\n",
    "# Print what events we're working with\n",
    "print(events.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:03:11.896175Z",
     "start_time": "2017-07-26T18:03:11.691511Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Plot the events\n",
    "plt.plot(events['cue1or5'], np.ones(len(events['cue1or5'])), 'o')\n",
    "plt.plot(events['cue2or4'], np.ones(len(events['cue2or4']))+1, 'o')\n",
    "plt.plot(events['feeder0'], np.ones(len(events['feeder0']))+2, 'o')\n",
    "plt.legend(['cue1or5', 'cue2or4', 'feeder0'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:03:49.651891Z",
     "start_time": "2017-07-26T18:03:49.436739Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Load spikes (.t and ._t) from this experiment\n",
    "spikes = nept.load_spikes(data_folder)\n",
    "\n",
    "# Plot the spikes\n",
    "for idx, spiketrain in enumerate(spikes):\n",
    "    plt.plot(spiketrain.time, np.ones(len(spiketrain.time))+idx, '|')\n",
    "plt.show()\n",
    "\n",
    "# Print the number of spikes emitted from each neuron\n",
    "for idx, spiketrain in enumerate(spikes):\n",
    "    print('Neuron', idx+1, 'has', len(spiketrain.time), 'spikes')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary\n",
    "\n",
    "Together, these three types of data can describe most data sets encountered in neuroscience. \n",
    "Putting them together in a simple visualization might look something like this, with LFP in black, \n",
    "position in cyan, epochs in red and green, and spikes in blue:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:06:02.185371Z",
     "start_time": "2017-07-26T18:06:01.971513Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "# Construct LFP\n",
    "time = np.arange(100)\n",
    "data = np.array(random.sample(range(100), 100)) / 100\n",
    "\n",
    "toy_lfp = nept.LocalFieldPotential(data, time)\n",
    "\n",
    "# Plot the LFP\n",
    "plt.plot(toy_lfp.time, toy_lfp.data, 'k')\n",
    "\n",
    "# Construct 1D position\n",
    "time = np.arange(100)\n",
    "x = np.array(random.sample(range(100), 100)) / -100\n",
    "\n",
    "toy_position = nept.Position(x, time)\n",
    "\n",
    "# Plot the position\n",
    "plt.plot(toy_position.time, toy_position.x, 'c')\n",
    "\n",
    "# Construct some epochs, which might be something like a light cue and a sound cue\n",
    "toy_light = nept.Epoch(np.array([[10, 15],\n",
    "                                 [70, 75]]))\n",
    "\n",
    "toy_sound = nept.Epoch(np.array([[50, 60]]))\n",
    "\n",
    "# Plot the epochs\n",
    "for light_start, light_stop in zip(toy_light.starts, toy_light.stops):\n",
    "    plt.axvspan(light_start, light_stop, color='g', alpha=0.5, lw=0)\n",
    "\n",
    "for sound_start, sound_stop in zip(toy_sound.starts, toy_sound.stops):\n",
    "    plt.axvspan(sound_start, sound_stop, color='r', alpha=0.5, lw=0)\n",
    "\n",
    "# Construct some spikes\n",
    "toy_spikes = [nept.SpikeTrain(np.array(random.sample(range(100), 80))),\n",
    "              nept.SpikeTrain(np.array(random.sample(range(100), 30))),\n",
    "              nept.SpikeTrain(np.array(random.sample(range(100), 50)))]\n",
    "\n",
    "# Plot the spikes\n",
    "for idx, spiketrain in enumerate(toy_spikes):\n",
    "    plt.plot(spiketrain.time, np.ones(len(spiketrain.time))+idx+0.5, '|', color='b', ms=20, mew=2)\n",
    "\n",
    "# Turn off axis\n",
    "plt.axis('off')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:07:40.670960Z",
     "start_time": "2017-07-26T18:07:40.371248Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# With real data, this might look something more like this:\n",
    "\n",
    "# Define our time of interest\n",
    "start = 1240\n",
    "stop = 1250\n",
    "\n",
    "# Slice the LFP to our time of interest\n",
    "sliced_lfp = lfp.time_slice(start, stop)\n",
    "\n",
    "# Plot the LFP\n",
    "plt.plot(sliced_lfp.time, sliced_lfp.data*1000, 'k')\n",
    "\n",
    "# Slice the position to our time of interest\n",
    "sliced_position = position.time_slice(start, stop)\n",
    "\n",
    "# Plot the position\n",
    "plt.plot(sliced_position.time, sliced_position.x*-0.001-0.3, 'c')\n",
    "\n",
    "# Let's only work with one event, our event of interest\n",
    "event_of_interest = events['pellet1']\n",
    "\n",
    "# Slice our event to the time of interest\n",
    "sliced_event = event_of_interest[(start < event_of_interest) & (event_of_interest < stop)]\n",
    "if len(sliced_event) % 2 != 0: # remove last sample if odd number of events\n",
    "    sliced_event = sliced_event[:-1]\n",
    "\n",
    "# For visualization purposes, let's say this event was present for a certain amount of time\n",
    "# and can be represented as an epoch, with a start and a stop\n",
    "event_epoch = nept.Epoch([sliced_event[::2], sliced_event[1::2]])\n",
    "\n",
    "# Plot the epoch\n",
    "for epoch_start, epoch_stop in zip(event_epoch.starts, event_epoch.stops):\n",
    "    plt.axvspan(epoch_start, epoch_stop, color='r', alpha=0.5, lw=0)\n",
    "   \n",
    "# Slice our spikes to the time of interest\n",
    "sliced_spikes = [spiketrain.time_slice(start, stop) for spiketrain in spikes]\n",
    "\n",
    "# Plot the spikes\n",
    "for idx, spiketrain in enumerate(sliced_spikes):\n",
    "    plt.plot(spiketrain.time, np.ones(len(spiketrain.time))+idx, '|', color='b', ms=20, mew=2)\n",
    "\n",
    "# Turn off axis\n",
    "plt.axis('off')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Real world examples\n",
    "\n",
    "Here are two examples that illustrate some simple operations. \n",
    "You should run them and make sure you understand how the raw data is transformed.\n",
    "For this, we will use a more involved data set R042-2013-08-18,\n",
    "which has recordings from hippocampal CA1 neurons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:08:55.470749Z",
     "start_time": "2017-07-26T18:08:55.458722Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# import necessary packages\n",
    "%matplotlib inline\n",
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import nept\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats\n",
    "\n",
    "# define where your data folder is located\n",
    "data_path = os.path.join(os.path.abspath('.'), 'data')\n",
    "data_folder = os.path.join(data_path, 'R042-2013-08-18')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:08:31.176009Z",
     "start_time": "2017-07-26T18:08:31.171006Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# load the info file, which contains experiment-specific information\n",
    "sys.path.append(data_folder)\n",
    "import r042d3 as info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:08:32.236214Z",
     "start_time": "2017-07-26T18:08:31.707838Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Load position (.nvt) from this experiment\n",
    "position = nept.load_position(os.path.join(data_folder, info.position_filename), info.pxl_to_cm)\n",
    "\n",
    "# Plot the position\n",
    "plt.plot(position.x, position.y, 'k.', ms=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:08:34.817461Z",
     "start_time": "2017-07-26T18:08:34.643356Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "left = position.time_slice(info.experiment_times['left_trials'].starts, \n",
    "                           info.experiment_times['left_trials'].stops)\n",
    "\n",
    "plt.plot(left.x, left.y, 'b.', ms=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:08:36.469326Z",
     "start_time": "2017-07-26T18:08:36.299706Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# This looks like right trials because the camera reverses the image.\n",
    "# We can fix this by reversing the y_axis\n",
    "\n",
    "plt.plot(left.x, left.y, 'b.', ms=1)\n",
    "plt.gca().invert_yaxis()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:08:38.976079Z",
     "start_time": "2017-07-26T18:08:38.972546Z"
    },
    "collapsed": true
   },
   "source": [
    "Now let's look at how these data types can work together\n",
    "to detect potential artifacts in the LFP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-26T18:09:17.457800Z",
     "start_time": "2017-07-26T18:09:13.571343Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Load LFP with good sharp-wave ripples\n",
    "lfp_swr = nept.load_lfp(os.path.join(data_folder, info.lfp_swr_filename))\n",
    "\n",
    "# Normalize the LFP\n",
    "lfp_norm = nept.LocalFieldPotential(scipy.stats.zscore(lfp_swr.data), lfp_swr.time)\n",
    "\n",
    "# Detect possible artifacts, save as epochs\n",
    "thresh = -8\n",
    "below_thresh = np.squeeze(lfp_norm.data < thresh)\n",
    "below_thresh = np.hstack([[False], below_thresh, [False]]) # pad start and end to properly compute the diff\n",
    "\n",
    "get_diff = np.diff(below_thresh.astype(int))\n",
    "\n",
    "artifact_starts = lfp_norm.time[np.where(get_diff == 1)[0]]\n",
    "artifact_stops = lfp_norm.time[np.where(get_diff == -1)[0]]\n",
    "\n",
    "artifacts = nept.Epoch([artifact_starts, artifact_stops])\n",
    "\n",
    "# Plot the LFP and artifact epochs\n",
    "plt.plot(lfp_norm.time, lfp_norm.data)\n",
    "\n",
    "for epoch_start, epoch_stop in zip(artifacts.starts, artifacts.stops):\n",
    "    plt.axvspan(epoch_start, epoch_stop, color='r', alpha=0.5, lw=1)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.4.5"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "497px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_position": {
    "height": "485px",
    "left": "0px",
    "right": "1074px",
    "top": "106px",
    "width": "206px"
   },
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}