Tuto_MNE_Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "N6GLbwRVhwgm"
   },
   "source": [
    "# ERF and spectral analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we start by importing the packages that we will use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "PDDNgQbRhwgo"
   },
   "outputs": [],
   "source": [
    "import mne\n",
    "import numpy as np # This is a basic python package that includes tools for manipulating the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "fWlRmx2Vhwgs"
   },
   "outputs": [],
   "source": [
    "#data_path = 'C:\\\\Users\\\\yhare\\\\mne_data\\\\MNE-brainstorm-data\\\\bst_raw\\\\'\n",
    "data_path = '/home/hyruuk/mne_data/MNE-brainstorm-data/bst_raw/'\n",
    "fname = 'clean_dataset_epochs.fif.gz'\n",
    "epochs = mne.read_epochs(data_path + fname, preload=True); # Here we load the cleaned data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "RpI-zP15hwhA",
    "outputId": "d608e9e1-73fb-4bdb-a6b5-96790e4eb763"
   },
   "source": [
    "### Plotting evoked data\n",
    "Now, we can try to see if there is a difference between the evoked activity in our two conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evoked = epochs.average() # The function average() creates an Evoked MNE object, which is the average of all the epochs.\n",
    "\n",
    "\n",
    "evoked.plot(); # First, we can make a \"butterfly\" view of the evoked activity (across both conditions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We clearly see a response that is phase-locked to the stimulus. Let's look at the topographic distribution of this evoked activity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evoked.plot_topomap();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the topography at 33ms. We can visually identify at least two generators. Are they somehow link to our experimental conditions ? Let's display them separately to confirm !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs['Right'].average().plot_topomap(title='Right');\n",
    "epochs['Left'].average().plot_topomap(title='Left');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like it ! That makes sense because we know that a stimulation of the right hand will induce activity in the left sensorimotor cortices, and conversely."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we could also plot the Event-Related Field (ERF), to compare the two conditions across time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "RxWbkuFGhwhT",
    "outputId": "c8d5bee2-dd18-4b2f-d8b7-20bf07f06c27"
   },
   "outputs": [],
   "source": [
    "# First, we have to build a dict containing the two conditions in separate entries\n",
    "conditions = [\"Right\", \"Left\"]\n",
    "evoked_dict = dict()\n",
    "\n",
    "# Now we run a loop that will collect and average epochs separately for each condition\n",
    "for condition in conditions:\n",
    "    evoked_dict[condition] = epochs[condition].average().filter(l_freq=None, h_freq=30, method='iir') # Note that we apply a\n",
    "                        # low-pass filter to the data (after averaging). This is a common practice to visualize data efficiently.\n",
    "                                                                                \n",
    "print(evoked_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now plot an ERF of each condition !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pick = epochs.ch_names.index('MLC21-4408') # I choose this sensor since it is just above sensorimotor areas.\n",
    "colors = dict(Right='green', Left='red')\n",
    "mne.viz.plot_compare_evokeds(evoked_dict, picks=pick, colors=colors);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's do the same for the same channel but in the other hemisphere\n",
    "pick = epochs.ch_names.index('MRC21-4408')\n",
    "colors = dict(Right='green', Left='red')\n",
    "\n",
    "mne.viz.plot_compare_evokeds(evoked_dict, picks=pick, colors=colors);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The activity looked like it is reversed ! Again, this is exactly what we expected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "Qk8yn-z-hwhb"
   },
   "source": [
    "### Plotting spectral (PSD) data\n",
    "Now that we did a quick overlook on the evoked activity, we should also perform some kind of spectral analysis.\n",
    "Quantifying Power Spectrum Densities (PSD) in various frequency bands can inform us about the power distribution of cannonical oscillations. Alternatively, by using Time-Frequency Maps (TF maps) we can distinguish between evoked (phase-locked to the stimuli) and induced (not phase-locked) activity, which may behave differently."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "Fw2hAE6Mhwhf",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "epochs.plot_psd(fmax=200); # First, let's have a look to the whole power spectrum. This is not very informative in our case \n",
    "                           # but it can help you to understand how your signal is structured (1/f, cut-off frequency etc...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs['Right'].plot_psd_topomap(); # You can also plot the topographies of frequency bands.\n",
    "epochs['Left'].plot_psd_topomap();  # Let's compare both conditions !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Meh, we don't see anything really convincing... Maybe a more lateralised activity in the 'Right' condition ?\n",
    "Let's explore the TF representation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "freqs = np.logspace(*np.log10([6, 35]), num=12) # Define the frequency range (i.e. frequential definition)\n",
    "n_cycles = freqs / 4.  # Different number of cycle per frequency (i.e. temporal definition)\n",
    "\n",
    "power, itc = mne.time_frequency.tfr_morlet(epochs, freqs=freqs, n_cycles=n_cycles, use_fft=True, return_itc=True, decim=3, n_jobs=1)\n",
    "\n",
    "power.plot_topo(baseline=baseline, mode='logratio', title='Average power', picks=None); # Plot the TF maps for each sensor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we want to compare Left and Right stimulations, we should build a TF map for each condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "power_right, itc_right = mne.time_frequency.tfr_morlet(epochs['Right'], freqs=freqs, n_cycles=n_cycles, use_fft=True, return_itc=True, decim=3, n_jobs=1)\n",
    "power_left, itc_left = mne.time_frequency.tfr_morlet(epochs['Left'], freqs=freqs, n_cycles=n_cycles, use_fft=True, return_itc=True, decim=3, n_jobs=1)\n",
    "# Note that here we only selected epochs from the Left or Right conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pick = power.ch_names.index('MLC21-4408')\n",
    "power_right.plot([pick], baseline=baseline, mode='logratio', title='Right ' + power.ch_names[pick]);\n",
    "power_left.plot([pick], baseline=baseline, mode='logratio', title='Left ' + power.ch_names[pick]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pick = power.ch_names.index('MRC21-4408')\n",
    "power_right.plot([pick], baseline=baseline, mode='logratio', title='Right ' + power.ch_names[pick]);\n",
    "power_left.plot([pick], baseline=baseline, mode='logratio', title='Left ' + power.ch_names[pick]);"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "default_view": {},
   "name": "Tuto_MNE_Analysis.ipynb",
   "provenance": [],
   "version": "0.3.2",
   "views": {}
  },
  "kernelspec": {
   "display_name": "electrophy",
   "language": "python",
   "name": "electrophy"
  },
  "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}