03_run_ica.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run independent components analysis (ICA)\n",
    "\n",
    "We clean the data using ICA and then apply some further preprocessing steps to the ICA cleaned data.\n",
    "\n",
    "**Performed steps**\n",
    "\n",
    "- load concatenated raw data\n",
    "- preprocessing the data for ICA\n",
    "    - highpass filtering the data at 1Hz\n",
    "    - downsampling the data to 250Hz (applying appropriate lowpass filter to prevent aliasing)\n",
    "- running ICA on non-broken EEG channels only\n",
    "- Using the EOG and ECG channels to automatically mark bad components in the ICA\n",
    "- manually inspect all components, using automatically marked ones as guidance\n",
    "    - select additional bad components\n",
    "    - disselect false positive bad components from automatic marking\n",
    "- Save the ICA object\n",
    "- Apply the ICA to the concatenated raw data \n",
    "    - load the data fresh!\n",
    "    - this is NOT the data preprocessed for ICA\n",
    "- preprocess ICA cleaned data with the following steps:\n",
    "    - bandpass filtering\n",
    "    - interpolation of bad channels\n",
    "    - re-referencing to average\n",
    "- Save the ICA cleaned, preprocessed data\n",
    "\n",
    "**Script can be run in INTERACTIVE or NON-INTERACTIVE mode!**\n",
    "\n",
    "- --> set the `INTERACTIVE` variable to True or False\n",
    "- --> you need to run with `INTERACTIVE=True` to screen the data and decide which components to reject. After running interactively once, these components to be rejected are saved and the script can be run non-interactively if you wish so."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-04T20:19:52.783866Z",
     "start_time": "2020-11-04T20:19:52.479576Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib qt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-04T20:19:54.905779Z",
     "start_time": "2020-11-04T20:19:52.785416Z"
    }
   },
   "outputs": [],
   "source": [
    "import itertools\n",
    "import multiprocessing\n",
    "import os\n",
    "import os.path as op\n",
    "import sys\n",
    "\n",
    "import mne\n",
    "import numpy as np\n",
    "\n",
    "from utils import BIDS_ROOT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-04T20:19:54.938768Z",
     "start_time": "2020-11-04T20:19:54.931514Z"
    }
   },
   "outputs": [],
   "source": [
    "# IO: Where to find the data\n",
    "# Where to find concatenated raw data\n",
    "fname_rawconcat_template = op.join(\n",
    "    BIDS_ROOT, \"derivatives\", \"sub-{0:02}\", \"sub-{0:02}_concat_eeg-raw.fif.gz\"\n",
    ")\n",
    "\n",
    "# Where to save ica-cleaned data and ICA object\n",
    "fname_rawclean_template = op.join(\n",
    "    BIDS_ROOT, \"derivatives\", \"sub-{0:02}\", \"sub-{0:02}_clean_eeg-raw.fif.gz\"\n",
    ")\n",
    "fname_ica_unscreened_template = op.join(\n",
    "    BIDS_ROOT,\n",
    "    \"derivatives\",\n",
    "    \"sub-{0:02}\",\n",
    "    \"sub-{0:02}_concat_eeg-unscreened-ica.fif.gz\",\n",
    ")\n",
    "fname_ica_screened_template = op.join(\n",
    "    BIDS_ROOT, \"derivatives\", \"sub-{0:02}\", \"sub-{0:02}_concat_eeg-screened-ica.fif.gz\"\n",
    ")\n",
    "fname_excluded_comps_template = op.join(\n",
    "    BIDS_ROOT,\n",
    "    \"derivatives\",\n",
    "    \"sub-{0:02}\",\n",
    "    \"sub-{0:02}_concat_eeg-excluded-ica-comps.txt\",\n",
    ")\n",
    "# Pack all names in a dict\n",
    "name_templates = dict()\n",
    "name_templates[\"rawconcat\"] = fname_rawconcat_template\n",
    "name_templates[\"rawclean\"] = fname_rawclean_template\n",
    "name_templates[\"ica_unscreened\"] = fname_ica_unscreened_template\n",
    "name_templates[\"ica_screened\"] = fname_ica_screened_template\n",
    "name_templates[\"excluded_comps\"] = fname_excluded_comps_template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-04T20:19:54.955036Z",
     "start_time": "2020-11-04T20:19:54.941371Z"
    }
   },
   "outputs": [],
   "source": [
    "# Set here whether to run the script interactively or not\n",
    "INTERACTIVE = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-04T20:19:56.416443Z",
     "start_time": "2020-11-04T20:19:56.403416Z"
    }
   },
   "outputs": [],
   "source": [
    "subjects = range(1, 41)\n",
    "\n",
    "# How many subjects to run over in parallel\n",
    "# NOTE: The ICA per subject itself is not parallellized ...\n",
    "# however, several sub-calls of the ICA may do multithreaded\n",
    "# MKL computations. By exporting `MKL_NUM_THREADS=1`\n",
    "# as an environment variable, that can be controlled.\n",
    "NJOBS = max(1, multiprocessing.cpu_count() - 8)\n",
    "\n",
    "# Whether or not to overwrite existing files\n",
    "overwrite = False\n",
    "\n",
    "# Using autoreject to automatically determine a threshold for\n",
    "# segment rejection prior to ICA fitting\n",
    "use_autoreject = False\n",
    "\n",
    "# For filtering later on\n",
    "LOW_CUTOFF = 0.1\n",
    "HIGH_CUTOFF = 40.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-04T20:19:58.142280Z",
     "start_time": "2020-11-04T20:19:58.129593Z"
    }
   },
   "outputs": [],
   "source": [
    "def preproc_data_for_ica(subj, name_templates):\n",
    "    \"\"\"Preprocess raw data for ICA.\n",
    "    \n",
    "    We are using fixed settings of highpassing at 1Hz and\n",
    "    downsampling to 100 Hz sampling frequency, using an \n",
    "    appropriate anti-aliasing filter before.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    subj : int\n",
    "        The subject id to work on.\n",
    "    name_templates : dict\n",
    "        A dictionary of string templates. Needs the following keys:\n",
    "        \"rawconcat\"\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    raw : mne.io.fiff.raw.Raw\n",
    "        The raw data preprocessed for ICA.\n",
    "        \n",
    "    \"\"\"\n",
    "    # Get the data\n",
    "    fname_rawconcat = name_templates[\"rawconcat\"].format(subj)\n",
    "    raw = mne.io.read_raw_fif(fname_rawconcat, preload=True)\n",
    "\n",
    "    # Preprocessing for ICA\n",
    "    # highpass filter\n",
    "    raw.filter(l_freq=1, h_freq=None)\n",
    "\n",
    "    # downsample\n",
    "    raw.resample(sfreq=100)\n",
    "\n",
    "    return raw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-04T20:20:00.003972Z",
     "start_time": "2020-11-04T20:19:59.955677Z"
    }
   },
   "outputs": [],
   "source": [
    "def run_ica(raw, subj, name_templates, use_autoreject):\n",
    "    \"\"\"Run the ica on raw data.\n",
    "\n",
    "    For running the ICA, we divide the raw data into fixed length epochs,\n",
    "    drop all epochs that overlap with \"BAD\" annotations, and then estimate\n",
    "    an amplitude threshold to further reject bad epochs (optional, see\n",
    "    `use_autoreject`). The ICA is then fit using extended infomax for the\n",
    "    EEG channels only. Lastly, we use the ECG and EOG channels from the\n",
    "    raw data to automatically mark presumably bad components.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    raw : mne.io.fiff.raw.Raw\n",
    "        The raw data preprocessed for ICA.\n",
    "    subj : int\n",
    "        The subject id to work on.\n",
    "    name_templates : dict\n",
    "        A dictionary of string templates. Needs the following keys:\n",
    "        \"ica\"\n",
    "    use_autoreject : bool\n",
    "        Whether or not to use \"autoreject\" to further clean the data\n",
    "        prior to ICA fitting.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    ica : mne.preprocessing.ica.ICA\n",
    "        The ica object.\n",
    "\n",
    "    \"\"\"\n",
    "    # Length of data chunks for artifact rejection in seconds\n",
    "    tstep = 1.0\n",
    "\n",
    "    # Automatically estimate an amplitude threshold beyond which epochs\n",
    "    # should be classified as bad: \"autoreject\"\n",
    "    if use_autoreject:\n",
    "        # make sure you ran pip install autoreject\n",
    "        import autoreject\n",
    "\n",
    "        # Make even length epochs out of the continuous data\n",
    "        events = mne.make_fixed_length_events(raw, duration=tstep)\n",
    "        epochs = mne.Epochs(\n",
    "            raw, events, tmin=0.0, tmax=tstep, reject_by_annotation=True, baseline=None\n",
    "        )\n",
    "\n",
    "        # Drop epochs overlapping with \"BAD\" annotation\n",
    "        epochs.drop_bad()\n",
    "\n",
    "        # Automatically find a rejection threshold\n",
    "        reject = autoreject.get_rejection_threshold(epochs, ch_types=[\"eeg\"])\n",
    "    else:\n",
    "        reject = None\n",
    "\n",
    "    # Initialize an ICA object, using extended infomax\n",
    "    ica = mne.preprocessing.ICA(\n",
    "        random_state=42, method=\"infomax\", fit_params=dict(extended=True)\n",
    "    )\n",
    "\n",
    "    # Get the channel indices of all channels that are *clean* and of type *eeg*\n",
    "    ica_chn_idxs = list(range(len(raw.ch_names)))\n",
    "\n",
    "    bad_idxs = [raw.ch_names.index(ii) for ii in raw.info[\"bads\"]]\n",
    "    eog_idxs = [raw.ch_names.index(ii) for ii in raw.ch_names if \"EOG\" in ii]\n",
    "    ecg_idxs = [raw.ch_names.index(ii) for ii in raw.ch_names if \"ECG\" in ii]\n",
    "\n",
    "    for chn_to_exclude in [bad_idxs, eog_idxs, ecg_idxs]:\n",
    "        for chn in chn_to_exclude:\n",
    "            if chn in ica_chn_idxs:\n",
    "                ica_chn_idxs.remove(chn)\n",
    "\n",
    "    # Fit our raw (high passed, downsampled) data to our ica object\n",
    "    # we pass the estimated rejection threshold to further exclude\n",
    "    # bad data from the estimation\n",
    "    ica.fit(\n",
    "        raw, picks=ica_chn_idxs, reject=reject, tstep=tstep, reject_by_annotation=True\n",
    "    )\n",
    "\n",
    "    # Automatically find artifact components using the EOG and ECG data\n",
    "    # NOTE: for find_bads_ecg, the default threshold of 0.25 is improved upon\n",
    "    # in MNE 0.21 by adapting it automatically to the sampling rate. See:\n",
    "    # https://github.com/mne-tools/mne-python/blob/ca3595f0b073ed4c75470eae634d046236007838/mne/preprocessing/ica.py#L1110-L1136\n",
    "    # Manually running that code yields a threshold of `0.5` for sf=100Hz\n",
    "    veog_idx, veog_scores = ica.find_bads_eog(raw, \"VEOG\")\n",
    "    heog_idx, heog_scores = ica.find_bads_eog(raw, \"HEOG\")\n",
    "    ecg_idx, ecg_scores = ica.find_bads_ecg(raw, \"ECG\", method=\"ctps\", threshold=0.5)\n",
    "\n",
    "    # exclude the automatically identified components\n",
    "    ica.exclude = list(set(np.concatenate([veog_idx, heog_idx, ecg_idx]).astype(int)))\n",
    "\n",
    "    return ica"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run ICA for all subjects\n",
    "\n",
    "Save the results as \"unscreened\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-10T18:33:20.737149Z",
     "start_time": "2020-07-10T18:33:20.724706Z"
    }
   },
   "outputs": [],
   "source": [
    "def run_parallel(subj, name_templates, overwrite, use_autoreject):\n",
    "    \"\"\"Run workflow for one subj and save.\"\"\"\n",
    "\n",
    "    # Handle existing files\n",
    "    fname_ica = name_templates[\"ica_unscreened\"].format(subj)\n",
    "    if op.exists(fname_ica):\n",
    "        if overwrite:\n",
    "            os.remove(fname_ica)\n",
    "        else:\n",
    "            # if file exists and we do not overwrite, return early\n",
    "            return\n",
    "\n",
    "    # compute\n",
    "    raw = preproc_data_for_ica(subj, name_templates)\n",
    "    ica = run_ica(raw, subj, name_templates, use_autoreject)\n",
    "\n",
    "    # save\n",
    "    fname_ica_unscreened = name_templates[\"ica_unscreened\"].format(subj)\n",
    "    ica.save(fname_ica_unscreened)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T08:40:00.449878Z",
     "start_time": "2020-07-02T08:33:46.166511Z"
    }
   },
   "outputs": [],
   "source": [
    "# Run the pipeline in parallel over subjects\n",
    "pool_inputs = itertools.product(\n",
    "    subjects, [name_templates], [overwrite], [use_autoreject]\n",
    ")\n",
    "\n",
    "with multiprocessing.Pool(NJOBS) as pool:\n",
    "    pool.starmap(run_parallel, pool_inputs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Screen ICA results, apply to non-processed raw data\n",
    "\n",
    "- Save the updated results as \"screened\"\n",
    "- Also save the to-be-rejected ICA components in a TXT file (each component\n",
    "  index on a separate line)\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non-interactive way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-09T13:58:44.666919Z",
     "start_time": "2020-07-09T13:58:44.650276Z"
    }
   },
   "outputs": [],
   "source": [
    "# NON-INTERACTIVE way to run the script\n",
    "if not INTERACTIVE:\n",
    "\n",
    "    for subj in subjects:\n",
    "\n",
    "        # Read components to exclude\n",
    "        fname_comps_to_exclude = name_templates[\"excluded_comps\"].format(subj)\n",
    "        if not op.exists(fname_comps_to_exclude):\n",
    "            raise RuntimeError(\n",
    "                \"The components to exlude file does not exist yet. \"\n",
    "                \"You first have to run the script interactively and save the files.\"\n",
    "            )\n",
    "\n",
    "        with open(fname_comps_to_exclude, \"r\") as fin:\n",
    "            lines = fin.readlines()\n",
    "            comps_to_exclude = [int(i.strip()) for i in lines]\n",
    "        \n",
    "        # Process the unscreened ICA, apply components to exclude -> screened ICA\n",
    "        fname_ica_unscreened = name_templates[\"ica_unscreened\"].format(subj)\n",
    "        ica = mne.preprocessing.read_ica(fname_ica_unscreened)\n",
    "\n",
    "        ica.exclude = comps_to_exclude\n",
    "\n",
    "        fname_ica_screened = name_templates[\"ica_screened\"].format(subj)\n",
    "        ica.save(fname_ica_screened)\n",
    "\n",
    "        # Apply screened ICA to unprocessed raw data -> clean raw\n",
    "        # NOTE: we call ica.apply on the full, unprocessed raw data with\n",
    "        # all channels (64 EEG, 2 EOG, 1 ECG). Once could think that this\n",
    "        # does not work, because to fit the ICA, we only passed the *clean*\n",
    "        # EEG channels (so 64 or less EEG). However, it *does* work, because\n",
    "        # MNE-Python's ica object stores the channels used for fitting, and\n",
    "        # when we call apply, the ica mixing is only applied to those channels\n",
    "        # that were also used for fitting. All other channels will remain as\n",
    "        # they are: which is good for EOG and ECG ... and the bad EEG channels\n",
    "        # will later be interpolated using the clean EEG channels.        \n",
    "        unproc_raw = mne.io.read_raw_fif(\n",
    "            name_templates[\"rawconcat\"].format(subj), preload=True\n",
    "        )\n",
    "        ica.apply(unproc_raw)\n",
    "\n",
    "        # Finish with filtering + interpolation + re-reference\n",
    "        # Highpass filtering\n",
    "        unproc_raw = unproc_raw.filter(l_freq=LOW_CUTOFF, h_freq=None)\n",
    "\n",
    "        # Lowpass filtering\n",
    "        unproc_raw = unproc_raw.filter(l_freq=None, h_freq=HIGH_CUTOFF)\n",
    "\n",
    "        # Interpolation\n",
    "        unproc_raw = unproc_raw.interpolate_bads()\n",
    "\n",
    "        # Re-referencing\n",
    "        unproc_raw = unproc_raw.set_eeg_reference(ref_channels=\"average\",\n",
    "                                                  projection=False, ch_type=\"eeg\")\n",
    "\n",
    "        # Save as cleaned data\n",
    "        unproc_raw.save(name_templates[\"rawclean\"].format(subj))\n",
    "\n",
    "    # This will terminate the script, preventing it to go to the interactive\n",
    "    # part.\n",
    "    print(\"\\n\\nDone! ... quitting early to prevent interactive part.\")\n",
    "    sys.exit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interactive way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-04T20:20:09.152909Z",
     "start_time": "2020-11-04T20:20:09.141066Z"
    }
   },
   "outputs": [],
   "source": [
    "## Increment the subj variable by ID 1 to 40\n",
    "subj = 1\n",
    "\n",
    "fname = name_templates[\"excluded_comps\"].format(subj)\n",
    "if op.exists(fname):\n",
    "    print(\n",
    "        \"This data has already been screened. \"\n",
    "        \"You may want to increase the `subj` variable.\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-04T20:20:11.613630Z",
     "start_time": "2020-11-04T20:20:11.605649Z"
    }
   },
   "outputs": [],
   "source": [
    "fname_ica_screened = name_templates[\"ica_screened\"].format(subj)\n",
    "\n",
    "if op.exists(fname_ica_screened):\n",
    "    print(\"Screened ICA file already exists. Increment `subj` variable?\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-04T20:21:29.036812Z",
     "start_time": "2020-11-04T20:20:14.689494Z"
    }
   },
   "outputs": [],
   "source": [
    "# Load data required for screening\n",
    "raw = preproc_data_for_ica(subj, name_templates)\n",
    "fname_ica_unscreened = name_templates[\"ica_unscreened\"].format(subj)\n",
    "ica = mne.preprocessing.read_ica(fname_ica_unscreened)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-04T20:21:36.923093Z",
     "start_time": "2020-11-04T20:21:36.911046Z"
    }
   },
   "outputs": [],
   "source": [
    "# Which ICA components were automatically marked to be excluded?\n",
    "ica.exclude"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-04T20:22:01.356030Z",
     "start_time": "2020-11-04T20:21:39.093857Z"
    }
   },
   "outputs": [],
   "source": [
    "# Screen the data, marking components to be included or rejected\n",
    "# I.e., double checking the automatic marking\n",
    "ica.plot_sources(raw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T08:42:03.018499Z",
     "start_time": "2020-07-02T08:42:03.014496Z"
    }
   },
   "outputs": [],
   "source": [
    "# Which ICA components are now to be excluded?\n",
    "ica.exclude"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T08:42:13.734612Z",
     "start_time": "2020-07-02T08:42:13.556043Z"
    }
   },
   "outputs": [],
   "source": [
    "# Save which components to be excluded\n",
    "fname = name_templates[\"excluded_comps\"].format(subj)\n",
    "if op.exists(fname) and not overwrite:\n",
    "    raise IOError(\"fname '{}' exists! Please double check.\".format(fname))\n",
    "else:\n",
    "    with open(fname, \"w\") as fout:\n",
    "        fout.write(\"\\n\".join([str(comp_idx) for comp_idx in ica.exclude]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T08:42:36.871594Z",
     "start_time": "2020-07-02T08:42:36.847603Z"
    }
   },
   "outputs": [],
   "source": [
    "# Save data as screened\n",
    "ica.save(fname_ica_screened)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T08:43:08.321307Z",
     "start_time": "2020-07-02T08:42:39.206771Z"
    }
   },
   "outputs": [],
   "source": [
    "# Apply ICA to unprocessed raw data\n",
    "ica = mne.preprocessing.read_ica(fname_ica_screened)\n",
    "unproc_raw = mne.io.read_raw_fif(name_templates[\"rawconcat\"].format(subj), preload=True)\n",
    "ica.apply(unproc_raw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inspect cleaned vs. non-cleaned data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T08:45:34.251240Z",
     "start_time": "2020-07-02T08:45:08.331917Z"
    }
   },
   "outputs": [],
   "source": [
    "clean_raw = unproc_raw.copy()\n",
    "unclean_raw = mne.io.read_raw_fif(\n",
    "    name_templates[\"rawconcat\"].format(subj), preload=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T08:45:56.544475Z",
     "start_time": "2020-07-02T08:45:56.529537Z"
    }
   },
   "outputs": [],
   "source": [
    "# temporarily deleting annotations for plotting\n",
    "clean_raw.annotations.delete(range(len(clean_raw.annotations.description)))\n",
    "unclean_raw.annotations.delete(range(len(unclean_raw.annotations.description)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T08:46:01.154011Z",
     "start_time": "2020-07-02T08:45:57.362223Z"
    }
   },
   "outputs": [],
   "source": [
    "clean_raw.plot(n_channels=len(clean_raw.ch_names), bad_color=(1, 0, 0), duration=20.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unclean_raw.plot(\n",
    "    n_channels=len(unclean_raw.ch_names), bad_color=(1, 0, 0), duration=20.0\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finish with filtering + interpolation + re-reference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T08:48:17.532281Z",
     "start_time": "2020-07-02T08:47:57.172386Z"
    }
   },
   "outputs": [],
   "source": [
    "# Highpass filtering\n",
    "unproc_raw = unproc_raw.filter(l_freq=LOW_CUTOFF, h_freq=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lowpass filtering\n",
    "unproc_raw = unproc_raw.filter(l_freq=None, h_freq=HIGH_CUTOFF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Interpolation\n",
    "unproc_raw = unproc_raw.interpolate_bads()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Re-referencing\n",
    "unproc_raw = unproc_raw.set_eeg_reference(ref_channels=\"average\", projection=False, ch_type=\"eeg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save as cleaned data\n",
    "unproc_raw.save(name_templates[\"rawclean\"].format(subj))"
   ]
  }
 ],
 "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
}