1-basic-sleep-spindle-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "79767ed2",
   "metadata": {},
   "source": [
    "## Tutorial on sleep spindle analysis using YASA\n",
    "(c) Simon Kern 2022\n",
    "\n",
    "Welcome to this short tutorial on getting started with Jupyter Notebook and YASA. This tutorial will try to provide a beginner-friendly instruction on how to use Jupyter notebooks to run sleep spindle analysis. It is be recommended to watch a short introduction on Jupyter notebooks before running this notebook to understand the basics. For German speakers this one seems well suited: https://www.youtube.com/watch?v=1S4Cgtkxqhs \n",
    "\n",
    "Generally, Jupyter notebooks allow to subdivide Python code into different cells that can be run individually and directly show the output of the given code. Additionally, cells with human-readable explanations (such as this cell itself) are possible. You can execute a cell by pressing SHIFT-ENTER. If it's a code cell, the code inside will be executed, if it's a human readable cell (i.e. a MarkDown cell) it will render the code. \n",
    "1. Try double-clicking this text to enter the edit mode. \n",
    "2. As you can see, you now can edit the text and formatting \n",
    "3. Press STRG+ENTER or SHIFT-ENTER to render the text again"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6da8614",
   "metadata": {},
   "source": [
    "### Getting started\n",
    "\n",
    "Before running the analysis we need to tell Python which packages it needs import. Select the cell below and press \"STRG+ENTER\" to execute it. You will see a star appearing left to it `In [*]` which will indicate the current cell is running. After a successful run it will display `In [1]` or any other number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38971eb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import mne # MNE is a package for processing EEG data in Python (e.g. loading EDF files)\n",
    "import yasa # the package we previously installed, which contains most of the code we needy\n",
    "import sleep_utils # this is a local package which is part of the GitHub repository\n",
    "print('importing done!')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "742d0579",
   "metadata": {},
   "source": [
    "#### Loading a sleep recording\n",
    "Now as an example we load a single EDF file and run a sleep spindle analysis on it using YASA.\n",
    "\n",
    "For this you first need to first adapt the `data_dir` where the EDFs can be found and then execute the cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7359116a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading data, this might take a while...\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'mne' 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~\\AppData\\Local\\Temp\\ipykernel_30912\\3239857560.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Loading data, this might take a while...'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 6\u001b[1;33m \u001b[0mraw\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmne\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mio\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mread_raw_edf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0medf_file\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpreload\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;31m# we load the file into memory using the function `mne.io.read_raw_edf`\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      7\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Done Loading'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'mne' is not defined"
     ]
    }
   ],
   "source": [
    "data_dir = 'Z:/Exercise_Sleep_Project/EDF Export EEG/' # adapt the path on the left to where the data can be found on your computer\n",
    "\n",
    "edf_file = data_dir + '/AA3_EX_AA3_EX_(1).edf'  # this is simply the first EDF file in the list\n",
    "\n",
    "print('Loading data, this might take a while...')\n",
    "raw = mne.io.read_raw_edf(edf_file, preload=True) # we load the file into memory using the function `mne.io.read_raw_edf`\n",
    "print('Done Loading')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "853cdd2b",
   "metadata": {},
   "source": [
    "Now using the `print` command, we can display some information about this file\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "701d6ac9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# the object with the name `raw` is now a reference to the EEG containing all data, channels and other information\n",
    "print(raw.info) # in the info object for example we will find e.g the sampling frequency and channel list"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbcadba1",
   "metadata": {},
   "source": [
    "### Using YASA for spindle detection\n",
    "\n",
    "Now we can use YASA to run some basic spindle analysis on one of the channels. We do this by simply calling this function:\n",
    "\n",
    "```Python\n",
    "yasa.spindles_detect(data_of_one_channel, sf=sampling_frequency)\n",
    "```\n",
    "\n",
    "To do this, we\n",
    "1. Extract the data of a channel of interest from the RAW object, rescaled to uV (microvolt)\n",
    "2. Apply a bipolar reference to this channel (assuming no reference had been applied previously)\n",
    "3. Substract the reference channel from the channel data\n",
    "4. Feed in this data into the spindle detect algorithm\n",
    "\n",
    "\n",
    "Additionally, there are many more parameters that the `spindle_detect` function takes. However it has some sensible defaults that we just use for now. For a list of parameters you can look here: [documentation:yasa.spindles_detect](https://raphaelvallat.com/yasa/build/html/generated/yasa.spindles_detect.html)\n",
    "\n",
    "These are the parameters that it uses by default\n",
    "```Python\n",
    "freq_sp      = (12,15)  # look for spindles between 12-15 Hz\n",
    "freq_broad   = (1,30)   # filter signal between 1 to 30 Hz\n",
    "duration     = (0.5, 2) # spindles should be between 0.5-2 seconds long\n",
    "min_distance = 0.5      # minimum distance between two spindles is 0.5 seconds\n",
    "\n",
    "```\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a683356",
   "metadata": {},
   "outputs": [],
   "source": [
    "ch = 'C4' # let's look at C4\n",
    "ref = 'M1' # take M1 as reference\n",
    "\n",
    "data_ch  = raw.get_data(ch, units='uV')  # retrieve data from the channel\n",
    "data_ref = raw.get_data(ref, units='uV') # retrieve data from the reference channel\n",
    "\n",
    "data = data_ch-data_ref # subtract the two signals (bipolar EEG reference)\n",
    "\n",
    "sampling_frequency = raw.info['sfreq'] # get the sampling frequency for these channels\n",
    "\n",
    "# run the spindle detection with the above mentioned defaults\n",
    "spindles = yasa.spindles_detect(data, sf=sampling_frequency) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0a36d85",
   "metadata": {},
   "source": [
    "The `spindles` variable now contains the results of the simple spindle detection, let's have a look at it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7cd5191",
   "metadata": {},
   "outputs": [],
   "source": [
    "spindles_summary = spindles.summary()\n",
    "spindles_summary # simple statement for printing the output nicely"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e05cdeaa",
   "metadata": {},
   "source": [
    "It found 793 spindles in total, and gives us quite a few parameters for these spindles (duration, amplitude, RMS, number of oscillations and so on! Looks promising.). Now lets calculate the mean of some of the values.\n",
    "\n",
    "The variable `spindles_summary` is a so called `DataFrame`. You can access the columns simply by putting their names in the the squared brackets. On the columns dozens of functions are defined that you can easily run, such as `max(), min(), mean(), median(), sum(), quantile()` and many many more. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f435c1da",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_duration = spindles_summary['Duration'].mean()\n",
    "median_amplitude = spindles_summary['Amplitude'].median()\n",
    "freq_quant_90 = spindles_summary['Frequency'].quantile(0.9)\n",
    "\n",
    "# let's calculate the spindle density by dividing the total number of spindles by the minutes\n",
    "minutes_in_data = raw.n_times / sampling_frequency / 60 # calculate minutes in data by number of samples divided by sampling_frequency and then by 60\n",
    "spindle_density = len(spindles_summary) / minutes_in_data # length of spindles_summary=792 divided by minutes_in_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5b144df",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets print out these results\n",
    "\n",
    "print('Mean spindle duration', mean_duration, 'seconds')\n",
    "print('Median spindle amplitude', median_amplitude, 'uV')\n",
    "print('90% quantile spindle frequency', freq_quant_90, 'Hz')\n",
    "print('spindle density', spindle_density, 'spindles/seconds')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d0fca14",
   "metadata": {},
   "source": [
    "#### Now it's your turn! \n",
    "\n",
    "in the cell below try to fill in the code to calculate and print\n",
    "- the mean abs power of all spindles\n",
    "- the mean spindle frequency\n",
    "- the 75 quantile of the absolute power"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f84e4451",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_abs_power = ....\n",
    "mean_spindle_freq = ...\n",
    "quant_75 = ...\n",
    "\n",
    "\n",
    "print(...)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "918f1ea8",
   "metadata": {},
   "source": [
    "## Run analysis on only certain sleep stages\n",
    "\n",
    "Of course now the detection algorithm will find a lot of spindles that are not in the are of interest, e.g. REM or other artefacts that are no spindles at all. Therefore we want to limit our analysis to just run on certain segments of the data. For this we are loading the hypnogram of the specific file into memory.\n",
    "\n",
    "- Make sure the hypnogram file (.TXT) is in the same folder as the EDF file and has the same name as the EDF file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5eae859",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make sure the hypnogram of this participant is in the same folder as the EDF\n",
    "hypnogram_file = data_dir + '/AA3_EX_AA3_EX_(1).edf.txt' # the hypnogram file\n",
    "hypnogram = sleep_utils.read_hypno(hypnogram_file)\n",
    "\n",
    "yasa.plot_hypnogram(hypnogram)\n",
    "yasa.plot_spectrogram(data[0], sampling_frequency)\n",
    "\n",
    "print('Now lets look at the hypnogram and the spectrogram of the data')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6176391f",
   "metadata": {},
   "source": [
    "#### running YASA with specific sleep stages only\n",
    "\n",
    "We can tell the YASA sleep spindle detection algorithm to only run on a certain segment of the data. We are basically doing the same analysis as above, we just pass some additional parameters to `yasa.spindles_detect()`. \n",
    "\n",
    "For this we first need to fit the hypnogram, which currently is in steps of 30 seconds such that each data point (i.e. 256 points per second) has one annotation. Then we provide the hypnogram to the `yasa.spindles_detect` function and tell it to only include certain stages.\n",
    "\n",
    "Here the stages are defined as following : \n",
    "```\n",
    "-2 = Unscored\n",
    "-1 = Artefact / Movement\n",
    " 0 = Wake\n",
    " 1 = N1 sleep\n",
    " 2 = N2 sleep\n",
    " 3 = N3 sleep\n",
    " 4 = REM sleep\n",
    "```\n",
    "IMPORTANT: Always make sure that these match your hypnogram! There is no standard in sleep research, and sometimes stages get annotated entirely differently. For example YASA only works with AASM, so there is no S4. In other annotation systems, `4` would annotate S4 and `5` would be REM, etc. So before analysing other data sets, make sure you have the data in the right format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fd0ca45",
   "metadata": {},
   "outputs": [],
   "source": [
    "stages_of_interest = [2, 3] # only run on NREM Stages 2 and 3 (SWS).\n",
    "hypno = yasa.hypno_upsample_to_data(hypnogram, sf_hypno=1/30, data=data, sf_data=sampling_frequency)\n",
    "\n",
    "spindles_NREM = yasa.spindles_detect(data, \n",
    "                                sf=sampling_frequency, \n",
    "                                hypno=hypno, \n",
    "                                include=stages_of_interest,\n",
    "                                freq_sp=(12, 15),\n",
    "                                verbose='INFO') \n",
    "\n",
    "print('Spindles found in S2 and S3')\n",
    "spindles_NREM.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8f6253d",
   "metadata": {},
   "source": [
    "As you can see, there were aroun ~180 spindles that were detected outside of NREM sleep. Now you could calculate the mean/median etc just as already done above, but this time using the DataFrame `spindles_NREM` as a basis."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75c47dfa",
   "metadata": {},
   "source": [
    "#### Now it's your turn (2)! \n",
    "Now as an exercise, in the cell below try to calculate \n",
    "- the difference between the mean spindle frequency in S2 and S3!\n",
    "- If you want, as a \"bonus exercise\", you can try to find a way to make another calculation:\n",
    "   1. Calculate \"slow\" spindles by limiting freq_sp to (12,14) \n",
    "   2. Calculate \"fast\" spindles by limiting freq_sp to (14,16)\n",
    "   3. See if there are more fast or slow spindles in the signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b932ca2",
   "metadata": {},
   "outputs": [],
   "source": [
    "spindles_summary_s2 = ...\n",
    "spindles_summary_s3 = ...\n",
    "\n",
    "mean_spindle_freq_s2 = ...\n",
    "mean_spindle_freq_s3 = ...\n",
    "\n",
    "print('The mean S2 spindle freq is', mean_spindle_freq_s2, 'Hz')\n",
    "print('The mean S3 spindle freq is', mean_spindle_freq_s3, 'Hz')\n",
    "\n",
    "print('There is a absolute difference of ', abs(mean_spindle_freq_s2 - mean_spindle_freq_s3), 'Hz of spindles in S2 and S3')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44ef6e59",
   "metadata": {},
   "source": [
    "### Summary\n",
    "\n",
    "In this notebook you saw how to perform the spindle analysis on one participant. However you are of course interested in not only one participant, but in several ones. Additionally, you want the data in a format that can be read by other programs, e.g. XLS or CSV. Obviously you do not need to run each participant manually and extract values by copy&paste, but you can do so with a simple script.\n",
    "\n",
    "The next notebook will provide the code to load and run data for several participants in one go."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "676b30e6",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}