notebooks/CodingLab1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Neural Data Analysis_\n",
    "\n",
    "Lecturer: Prof. Dr. Philipp Berens\n",
    "\n",
    "Tutors: Sarah Strauss, Ziwei Huang\n",
    "\n",
    "Summer term 2021\n",
    "\n",
    "Due date: 2021-04-26, 9am\n",
    "\n",
    "Student names: *FILL IN YOUR NAMES HERE*\n",
    "\n",
    "# Coding Lab 1\n",
    "\n",
    "Download the data file ```nda_ex1.csv``` from ILIAS and save it in a subfolder ```../data/```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import signal\n",
    "from sklearn.decomposition import PCA\n",
    "import scipy as sp\n",
    "import itertools as it\n",
    "sns.set_style('whitegrid')\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The recording contains extracellular recordings from 4 channels. Details of the experimental procedure were given during the lecture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fs = 30000     # sampling rate of the signal in Hz\n",
    "dt = 1/Fs\n",
    "gain = .5      # gain of the signal\n",
    "x = pd.read_csv('../data/nda_ex_1.csv', header=0, names=('Ch1', 'Ch2', 'Ch3', 'Ch4'))  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Ch1</th>\n",
       "      <th>Ch2</th>\n",
       "      <th>Ch3</th>\n",
       "      <th>Ch4</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>1.920000e+07</td>\n",
       "      <td>1.920000e+07</td>\n",
       "      <td>1.920000e+07</td>\n",
       "      <td>1.920000e+07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>3.600331e+00</td>\n",
       "      <td>-8.850918e-01</td>\n",
       "      <td>2.864284e-01</td>\n",
       "      <td>2.210982e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>5.824474e+02</td>\n",
       "      <td>6.014818e+02</td>\n",
       "      <td>6.464363e+02</td>\n",
       "      <td>6.126105e+02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>-3.607000e+03</td>\n",
       "      <td>-3.739000e+03</td>\n",
       "      <td>-3.871000e+03</td>\n",
       "      <td>-3.750000e+03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>-3.460000e+02</td>\n",
       "      <td>-3.610000e+02</td>\n",
       "      <td>-3.950000e+02</td>\n",
       "      <td>-3.640000e+02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>1.200000e+01</td>\n",
       "      <td>8.000000e+00</td>\n",
       "      <td>-1.000000e+00</td>\n",
       "      <td>1.000000e+01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>3.650000e+02</td>\n",
       "      <td>3.720000e+02</td>\n",
       "      <td>4.010000e+02</td>\n",
       "      <td>3.810000e+02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>2.873000e+03</td>\n",
       "      <td>3.004000e+03</td>\n",
       "      <td>3.099000e+03</td>\n",
       "      <td>3.017000e+03</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                Ch1           Ch2           Ch3           Ch4\n",
       "count  1.920000e+07  1.920000e+07  1.920000e+07  1.920000e+07\n",
       "mean   3.600331e+00 -8.850918e-01  2.864284e-01  2.210982e+00\n",
       "std    5.824474e+02  6.014818e+02  6.464363e+02  6.126105e+02\n",
       "min   -3.607000e+03 -3.739000e+03 -3.871000e+03 -3.750000e+03\n",
       "25%   -3.460000e+02 -3.610000e+02 -3.950000e+02 -3.640000e+02\n",
       "50%    1.200000e+01  8.000000e+00 -1.000000e+00  1.000000e+01\n",
       "75%    3.650000e+02  3.720000e+02  4.010000e+02  3.810000e+02\n",
       "max    2.873000e+03  3.004000e+03  3.099000e+03  3.017000e+03"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1: Filter Signal\n",
    "\n",
    "In order to detect action potentials, the first step is to filter out low frequency fluctuations (LFP) and high frequency noise. Determine appropriate filter settings and implement the filtering in the function ```filterSignal()```. A typical choice for this task would be a butterworth filter. Plot a segment of the raw signal and the filtered signal for all four channels with matching y-axis. The segment you choose should contain spikes. When you apply the function also test different filter settings. Hint: Make sure to use a filter that compensates for phase shifting.\n",
    "\n",
    "*Grading: 2 pts*\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def filterSignal(x, Fs, low, high):\n",
    "# Filter raw signal\n",
    "#   y = filterSignal(x, Fs, low, high) filters the signal x. Each column in x is one\n",
    "#   recording channel. Fs is the sampling frequency. low and high specify the passband in Hz.\n",
    "#   The filter delay is compensated in the output y.\n",
    "\n",
    "        \n",
    "\n",
    "    return y\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'y' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-5-3de16f8722e2>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mxf\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfilterSignal\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mFs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m500\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m4000\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32m<ipython-input-4-9706818744d6>\u001b[0m in \u001b[0;36mfilterSignal\u001b[1;34m(x, Fs, low, high)\u001b[0m\n\u001b[0;32m      7\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      8\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 9\u001b[1;33m     \u001b[1;32mreturn\u001b[0m \u001b[0my\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'y' is not defined"
     ]
    }
   ],
   "source": [
    "xf = filterSignal(x, Fs, 500, 4000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 8))\n",
    "\n",
    "T = 100000\n",
    "t = np.arange(0,T) * dt \n",
    "\n",
    "for i, col in enumerate(xf):\n",
    "    plt.subplot(4,2,2*i+1)\n",
    "    plt.plot(t,x[col][0:T],linewidth=.5)\n",
    "    plt.ylim((-1000, 1000))\n",
    "    plt.xlim((0,3))\n",
    "    plt.ylabel('Voltage')\n",
    "    \n",
    "    \n",
    "    plt.subplot(4,2,2*i+2)\n",
    "    plt.plot(t,xf[col][0:T],linewidth=.5)\n",
    "    plt.ylim((-400, 250))\n",
    "    plt.xlim((0,3))\n",
    "    plt.ylabel('Voltage')\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2: Detect action potentials\n",
    "\n",
    "Action potentials are usually detected by finding large-amplitude deflections in the continuous signal. A good choice of threshold for detecting spikes is important. If it is too low, you will detect too many low amplitude events (noise); if it is too high, you run the risk of missing good spikes. Implement an automatic procedure to obtain a reasonable threshold and detect the times when spikes occurred in the function ```detectSpikes()``` . Plot a segment of the filtered signal for all four channels with matching y-axis and indicate the time points where you detected spikes. Plot the threshold. Are the detected time points well aligned with peaks in the signal? Hint: Pay attention to the sign of the signal, what does it mean that it's an extracellular recording?\n",
    "\n",
    "*Grading: 3 pts*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def detectSpikes(x,Fs):\n",
    "# Detect spikes\n",
    "# s, t = detectSpikes(x,Fs) detects spikes in x, where Fs the sampling\n",
    "#   rate (in Hz). The outputs s and t are column vectors of spike times in\n",
    "#   samples and ms, respectively. By convention the time of the zeroth\n",
    "#   sample is 0 ms.\n",
    "\n",
    "\n",
    "\n",
    "    return (s, t)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "T = xf.shape[0]\n",
    "s, t = detectSpikes(xf.values(),Fs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(7, 8))\n",
    "\n",
    "tt = np.arange(0,T) * dt \n",
    "\n",
    "for i, col in enumerate(xf):\n",
    "    plt.subplot(4,1,i+1)\n",
    "    plt.plot(tt,xf[col],linewidth=.5)\n",
    "    plt.plot(tt[s],xf[col][s],'r.')\n",
    "    plt.ylim((-400, 400))\n",
    "    plt.xlim((0.025,0.075))\n",
    "    plt.ylabel('Voltage')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 3: Extract waveforms\n",
    "For later spike sorting we need the waveforms of all detected spikes. Extract the waveforms segments (1 ms) on all four channels for each spike time (as a result each spike is represented by a 4x30 element matrix). Implement this procedure in the function ```extractWaveforms()```. Plot (a) the first 100 spikes you detected and (b) the 100 largest spikes you detected. Are there a lot of very small spikes (likely noise) among your detected spikes? If so your threshold may be too low. Can you see obvious artifacts, not looking like spikes at all?\n",
    "\n",
    "*Grading: 2 pts*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extractWaveforms(x, s):\n",
    "# Extract spike waveforms.\n",
    "#   w = extractWaveforms(x, s) extracts the waveforms at times s (given in\n",
    "#   samples) from the filtered signal x using a fixed window around the\n",
    "#   times of the spikes. The return value w is a 3d array of size\n",
    "#   length(window) x #spikes x #channels.\n",
    "\n",
    "\n",
    "    return w\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "w = extractWaveforms(xf.as_matrix(),s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot first 100 spike waveforms, include 10 ms before the peak, and 20 ms after."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.arange(-10,20) * dt * 1000 \n",
    "\n",
    "plt.figure(figsize=(11, 8))\n",
    "\n",
    "for i, col in enumerate(xf):\n",
    "    plt.subplot(2,2,i+1)\n",
    "    plt.plot(t,w[:,1:100,i],'k', linewidth=1)\n",
    "    plt.ylim((-500, 250))\n",
    "    plt.xlim((-0.33,0.66))\n",
    "    plt.ylabel('Voltage')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot largest 100 spike waveforms, include 10 ms before the peak, and 20 ms after."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = np.argsort(np.min(np.min(w,axis=2),axis=0))\n",
    "\n",
    "\n",
    "t = np.arange(-10,20) * dt * 1000 \n",
    "\n",
    "plt.figure(figsize=(11, 8))\n",
    "for i, col in enumerate(xf):\n",
    "    plt.subplot(2,2,i+1)\n",
    "    plt.plot(t,w[:,idx[0:100],i],'k', linewidth=1)\n",
    "    plt.ylim((-1000, 500))\n",
    "    plt.xlim((-0.33,0.66))\n",
    "    plt.ylabel('Voltage')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 4: Extract features using PCA\n",
    "Compute the first three PCA features on each channel separately in ```extractFeatures()```. You can use a available PCA implementation or implement it yourself. After that, each spike is represented by a 12 element vector. Compute the fraction of variance captured by these three PCs.\n",
    "Plot scatter plots for all pairwise combinations of 1st PCs. Do you see clusters visually? \n",
    "\n",
    "*Grading: 2+1 pts*\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def extractFeatures(w):\n",
    "# Extract features for spike sorting.\n",
    "#   b = extractFeatures(w) extracts features for spike sorting from the\n",
    "#   waveforms in w, which is a 3d array of size length(window) x #spikes x\n",
    "#   #channels. The output b is a matrix of size #spikes x #features.\n",
    "#   The implementation should do PCA on the waveforms of each channel\n",
    "#   separately and uses the first three principal components. Thus, we get\n",
    "#   a total of 12 features. Also, the varianced explained by the 3 features per channel\n",
    "#   should be computed.\n",
    "    \n",
    "\n",
    "    \n",
    "\n",
    "    return b\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = extractFeatures(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 6))\n",
    "plt.suptitle('Scatter plots',fontsize=20)\n",
    "\n",
    "idx = [0, 3, 6, 9]\n",
    "p = 1\n",
    "labels = ['Ch1','Ch2','Ch3','Ch4']\n",
    "for i in np.arange(0,4):\n",
    "    for j in np.arange(i+1,4):\n",
    "        ax = plt.subplot(2,3,p, aspect='equal')\n",
    "        plt.plot(b[:,idx[i]],b[:,idx[j]],'.k', markersize=.7) \n",
    "        plt.xlabel(labels[i])\n",
    "        plt.ylabel(labels[j])\n",
    "        plt.xlim((-1500,1500))\n",
    "        plt.ylim((-1500,1500))\n",
    "        ax.set_xticks([])\n",
    "        ax.set_yticks([])\n",
    "        p = p+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('../data/nda_ex_1_features',b)\n",
    "np.save('../data/nda_ex_1_spiketimes',s)\n",
    "np.save('../data/nda_ex_1_waveforms',w)"
   ]
  }
 ],
 "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.7.3"
  },
  "nteract": {
   "version": "0.13.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}