01_raw_inspection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Raw Inspection\n",
    "\n",
    "Inspecting and screening the raw data, making annotations about:\n",
    "\n",
    "- bad segments across time\n",
    "- bad channels across the overall data\n",
    "\n",
    "We will also automatically set the block breaks, and the recorded time\n",
    "before the start and after the end of the experiment to \"BAD\".\n",
    "\n",
    "Note 2020-07-11: could have used [mne.preprocessing.mark_flat](https://github.com/mne-tools/mne-python/blob/76fc838e954e70697f33cab0edbca61872ca1c23/mne/preprocessing/flat.py#L14-L104https://github.com/mne-tools/mne-python/blob/76fc838e954e70697f33cab0edbca61872ca1c23/mne/preprocessing/flat.py#L14-L104) to enhance visual inspection with some automatic \"pre\"-marking.\n",
    "\n",
    "**THIS CODE IS MEANT TO BE RUN INTERACTIVELY**\n",
    "\n",
    "**VISUAL INSPECTION BY AN EXPERT**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib qt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "import os\n",
    "import os.path as op\n",
    "\n",
    "import mne\n",
    "import numpy as np\n",
    "\n",
    "from utils import BIDS_ROOT, task_not_present_for_subject"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# IO: Where to find the data\n",
    "eeg_path_template = op.join(\n",
    "    BIDS_ROOT, \"sub-{0:02}\", \"eeg\", \"sub-{0:02}_task-{1}_eeg.vhdr\"\n",
    ")\n",
    "\n",
    "# Where to save the annotations and bad channels\n",
    "fname_annots_template = op.join(\n",
    "    BIDS_ROOT, \"derivatives\", \"sub-{0:02}\", \"sub-{0:02}_task-{1}_annotations.txt\"\n",
    ")\n",
    "fname_channels_template = op.join(\n",
    "    BIDS_ROOT, \"derivatives\", \"sub-{0:02}\", \"sub-{0:02}_task-{1}_badchannels.txt\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a list of all data files to screen\n",
    "subjects = range(1, 41)\n",
    "tasks = [\"ActiveFixed\", \"ActiveVariable\", \"YokedFixed\", \"YokedVariable\", \"description\"]\n",
    "\n",
    "combinations = list(itertools.product(subjects, tasks))\n",
    "\n",
    "for combo in combinations:\n",
    "    if task_not_present_for_subject(*combo):\n",
    "        combinations.remove(combo)\n",
    "\n",
    "print(\"{} datsets to screen overall\".format(len(combinations)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# increment this index whenever one subject/task combination has been screened\n",
    "current_work_index = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the data for the current work index\n",
    "subject, task = combinations[current_work_index]\n",
    "fname_annots = fname_annots_template.format(subject, task)\n",
    "fname_channels = fname_channels_template.format(subject, task)\n",
    "if op.exists(fname_annots) or op.exists(fname_channels):\n",
    "    print('Data exists. You may want to increment \"current_work_index\".\\n\\n\\n')\n",
    "\n",
    "fpath = eeg_path_template.format(subject, task)\n",
    "raw = mne.io.read_raw_brainvision(fpath, preload=True)\n",
    "\n",
    "# Suppress an automatic \"average reference\"\n",
    "raw.set_eeg_reference(ref_channels=[])\n",
    "\n",
    "# Set the EOG and ECG channels to their type\n",
    "raw.set_channel_types({\"ECG\": \"ecg\", \"HEOG\": \"eog\", \"VEOG\": \"eog\"})\n",
    "\n",
    "# Set a standard montage for plotting later\n",
    "# NOTE: Could potentially set the \"true\" (i.e., measured, digitized) electrode\n",
    "# coordinates from the experiment here ... but template positions seem fine\n",
    "montage = mne.channels.make_standard_montage(\"standard_1020\")\n",
    "\n",
    "raw.set_montage(montage)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get segments for 1) block breaks and 2) pre/post experiment recording\n",
    "general_description = \"BAD_break\"\n",
    "buffer_plus_minus = 1.0  # buffer in seconds around breaks that won't be marked\n",
    "orig_time = raw.annotations.orig_time\n",
    "\n",
    "\n",
    "def get_stim_onset(raw, stim, nth_stim=0):\n",
    "    \"\"\"Help to find onset of a stimulus in the data.\"\"\"\n",
    "    idx = (raw.annotations.description == stim).nonzero()[0][nth_stim]\n",
    "    return idx, raw.annotations.onset[idx]\n",
    "\n",
    "\n",
    "# Get data for pre experiment\n",
    "_, start_onset = get_stim_onset(raw, \"Stimulus/S  1\")\n",
    "recording_onset = raw.first_samp / raw.info[\"sfreq\"]\n",
    "start_duration = (start_onset - recording_onset) - buffer_plus_minus\n",
    "if start_duration < 0:\n",
    "    start_duration = start_onset - recording_onset\n",
    "\n",
    "raw.annotations.append(recording_onset, start_duration, general_description)\n",
    "\n",
    "# Get data for block breaks ... and post experiment\n",
    "n_blocks = 5\n",
    "for block in range(n_blocks):\n",
    "    break_idx, break_onset = get_stim_onset(raw, \"Stimulus/S 23\", nth_stim=block)\n",
    "    # next event after a break marks the end of the break\n",
    "    if block < n_blocks - 1:\n",
    "        break_offset = raw.annotations.onset[break_idx + 1]\n",
    "        break_duration = (break_offset - break_onset) - buffer_plus_minus\n",
    "        if break_duration < 0:\n",
    "            break_duration = break_offset - break_onset\n",
    "\n",
    "    # unless for the last block, where we can extend the bad segment until\n",
    "    # the end of the recording\n",
    "    elif block == n_blocks - 1:\n",
    "        break_offset = raw.last_samp / raw.info[\"sfreq\"]\n",
    "        break_duration = break_offset - break_onset\n",
    "\n",
    "    if break_duration > buffer_plus_minus:\n",
    "        break_onset += buffer_plus_minus\n",
    "\n",
    "    raw.annotations.append(break_onset, break_duration, general_description)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Delete all other annotations except the BAD ones\n",
    "# ... for focusing on data quality during visual inspection\n",
    "# irrespective on current condition in the experiment\n",
    "bad_idxs = [descr.startswith(\"BAD\") for descr in raw.annotations.description]\n",
    "to_delete = np.arange(len(raw.annotations))[~np.array(bad_idxs)]\n",
    "raw.annotations.delete(to_delete)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Downsample the data to speed up plotting\n",
    "raw.resample(sfreq=250)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Inspect power spectral density of the data\n",
    "raw.plot_psd(reject_by_annotation=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the data, interactively annotating\n",
    "fig = raw.plot(n_channels=len(raw.ch_names), bad_color=(1, 0, 0), duration=20.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Whether to overwrite or not\n",
    "overwrite = False\n",
    "\n",
    "# Save the marked annotations and bad channels\n",
    "if not op.exists(fname_annots) or overwrite:\n",
    "    os.makedirs(op.split(fname_annots)[0], exist_ok=True)\n",
    "    raw.annotations.save(fname_annots)\n",
    "else:\n",
    "    print(\n",
    "        '{} already exists. Maybe increment \"current_work_index\".'.format(fname_annots)\n",
    "    )\n",
    "\n",
    "if not op.exists(fname_channels) or overwrite:\n",
    "    os.makedirs(op.split(fname_channels)[0], exist_ok=True)\n",
    "    with open(fname_channels, \"w\") as fout:\n",
    "        lines = \"\\n\".join(raw.info[\"bads\"])\n",
    "        fout.writelines(lines)\n",
    "else:\n",
    "    print(\n",
    "        '{} already exists. Maybe increment \"current_work_index\".'.format(\n",
    "            fname_channels\n",
    "        )\n",
    "    )"
   ]
  },
  {
   "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.7.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}