Preprocess (collapsed).ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bb96ec5-cfea-46ca-8e31-e5bdf7a1e1da",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import re\n",
    "import mne\n",
    "from autoreject import AutoReject\n",
    "from mne_icalabel import label_components\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88f3a0d6-4a5a-4983-bca4-e6fbc703d619",
   "metadata": {},
   "outputs": [],
   "source": [
    "#mne.viz.set_browser_backend('qt')\n",
    "#mne.set_config('MNE_BROWSER_BACKEND', 'qt')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db4f984e-275c-47e8-87f1-ce08d1b231a9",
   "metadata": {},
   "source": [
    "# Artifact correction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "210583c5-978c-4160-b7ab-99a050fe9d09",
   "metadata": {},
   "outputs": [],
   "source": [
    "SUBJECTS = [\"R0782_20220411_CrosslingViz_0003\", \"R0810_20230120_CrosslingViz_0004\", \"R0811_20230120_CrosslingViz_0005\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e502c44-6aee-4962-b8d9-77624ea968d3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Load data\n",
    "\n",
    "subject = SUBJECTS[0]\n",
    "file = f'CrossLingViz_EEG/{subject}/{subject}.vhdr'\n",
    "\n",
    "raw = mne.io.read_raw_brainvision(file, preload=True)\n",
    "\n",
    "# Filtering\n",
    "\n",
    "raw_filt = raw.copy().filter(1, 30)\n",
    "raw_filt.set_channel_types(mapping={'VEOG': 'eog'});\n",
    "\n",
    "raw_filt.pick_types(eeg=True)\n",
    "#raw_filt.plot(butterfly=True)\n",
    "\n",
    "# Bad channel interpolation\n",
    "\n",
    "#raw_filt.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8464ed0d-e882-4a4a-a527-49edf6633a7e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#raw_filt.info['bads'].append('FC5') \n",
    "#raw_filt.interpolate_bads(reset_bads=True)\n",
    "\n",
    "# ICA correction\n",
    "\n",
    "ica = mne.preprocessing.ICA(n_components=25, max_iter='auto', random_state=97)\n",
    "ica.fit(raw_filt) \n",
    "ica\n",
    "\n",
    "ic_labels = label_components(raw_filt, ica, method=\"iclabel\")\n",
    "labels = ic_labels[\"labels\"]\n",
    "\n",
    "exclude_idx = [idx for idx, label in enumerate(labels) if label not in [\"brain\", \"other\"]]\n",
    "\n",
    "ica.exclude = exclude_idx\n",
    "reconst_raw=raw_filt.copy()\n",
    "ica.apply(reconst_raw)\n",
    "\n",
    "reconst_raw = reconst_raw.pick_types(eeg=True)\n",
    "\n",
    "print(\"exclude following components for subject\", subject, exclude_idx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "585c0443-056d-41ed-90ea-49dbf5a22aaa",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df750da7-95da-42dd-963c-c83029cca608",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%notebook inline\n",
    "ica.plot_components()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1804b3da-7860-4ec1-9fca-67243f2c7188",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "raw_filt.load_data()\n",
    "ica.plot_sources(raw_filt, show_scrollbars=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6c9e6b8-8ebc-416b-8ef8-b81a6a980ea9",
   "metadata": {},
   "outputs": [],
   "source": [
    "ica.exclude = [0, 1, 3, 7]\n",
    "reconst_raw=raw_filt.copy()\n",
    "ica.apply(reconst_raw)\n",
    "\n",
    "reconst_raw = reconst_raw.pick_types(eeg=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2b6b014-976f-4f4a-b734-e3fd7c6f918c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "exclude_idx\n",
    "\n",
    "#809: [0, 1, 3, 7]\n",
    "#810: [0, 1, 2, 3, 20, 24]\n",
    "#811: [0, 2, 6, 14, 15, 16]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9e72626-4a6a-4724-a75c-c1cd60d65952",
   "metadata": {},
   "outputs": [],
   "source": [
    "reconst_raw.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1828bec0-9789-4d60-abe4-8c7dbb21da62",
   "metadata": {},
   "source": [
    "# Epochs (OPTO alignment)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "052e5889-599e-4365-a194-8c60a7f1928e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "opto = raw.copy()\n",
    "opto.pick_channels(['OPTO'])\n",
    "\n",
    "# OPTO epochs\n",
    "events =  mne.events_from_annotations(raw);\n",
    "opto_epochs = mne.Epochs(opto, events[0], tmin=-0.1, tmax=1, baseline=(None, 0.), preload=True)\n",
    "\n",
    "# Get delays\n",
    "arr = []\n",
    "for i in range(len(opto_epochs)):\n",
    "    d = np.squeeze(opto_epochs[i].get_data())\n",
    "\n",
    "    bl_M = d[0:50].mean() \n",
    "    bl_V = d[0:50].var()\n",
    "\n",
    "    th = bl_M + (20e3 * bl_V)\n",
    "\n",
    "    try:\n",
    "        idx = np.where(d>=th)[0][0]\n",
    "        shft = (idx - 50) * (1 /  opto_epochs.info['sfreq'] )\n",
    "        \n",
    "    except IndexError:\n",
    "        shft = 0\n",
    "        \n",
    "    arr.append(shft)\n",
    "\n",
    "events = events[0]\n",
    "\n",
    "events_shft = []\n",
    "\n",
    "for i in range(1, len(events)):\n",
    "    events_shft[i:i+1] = mne.event.shift_time_events(events[i:i+1], ids=None, tshift=arr[i-1], sfreq=raw.info['sfreq'])\n",
    "\n",
    "events_shft = np.vstack(events_shft)\n",
    "\n",
    "# EEG epochs\n",
    "Eng = pd.read_excel(\"G:\\내 드라이브\\Research\\Decoding\\PsychopyViz\\stimuli_Engv.xlsx\")\n",
    "Kor = pd.read_excel(\"G:\\내 드라이브\\Research\\Decoding\\PsychopyViz\\stimuli_Korv.xlsx\")\n",
    "\n",
    "try:\n",
    "    event_id={'RespYes': 1, 'RespNo': 2} # if all responses = yes\n",
    "    \n",
    "    for i in range(64):\n",
    "        event_id[Eng[\"word\"][i]] = int(Eng[\"trigger\"][i])\n",
    "    \n",
    "    for i in range(64):\n",
    "        event_id[Kor[\"word\"][i]] = int(Kor[\"trigger\"][i])\n",
    "    \n",
    "    epochs_aligned = mne.Epochs(reconst_raw, events_shft, event_id, tmin=-0.3, tmax=1, baseline=(None, 0.), preload=True, event_repeated='drop')\n",
    "\n",
    "except ValueError:\n",
    "    event_id={'RespYes': 1}\n",
    "    \n",
    "    for i in range(64):\n",
    "        event_id[Eng[\"word\"][i]] = int(Eng[\"trigger\"][i])\n",
    "    \n",
    "    for i in range(64):\n",
    "        event_id[Kor[\"word\"][i]] = int(Kor[\"trigger\"][i])\n",
    "    \n",
    "    epochs_aligned = mne.Epochs(reconst_raw, events_shft, event_id, tmin=-0.3, tmax=1, baseline=(None, 0.), preload=True, event_repeated='drop')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d435055-5acc-4b7c-b94e-0ccb49ebbabb",
   "metadata": {},
   "source": [
    "# Final cleaning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41d93d71-c4a2-497f-b23f-accbdcdba5d0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Keep trials w/o responses\n",
    "\n",
    "epochs_idx=[]\n",
    "\n",
    "for i in range(len(epochs_aligned)-1):\n",
    "    if epochs_aligned.events[i+1][2]!=1 and epochs_aligned.events[i+1][2]!=2:\n",
    "        epochs_idx.append(i)\n",
    "        \n",
    "epochs_noresp = epochs_aligned[epochs_idx]\n",
    "\n",
    "# Autoreject\n",
    "\n",
    "ar = AutoReject()\n",
    "ar.fit(epochs_noresp)\n",
    "\n",
    "epochs_ar, reject_log = ar.transform(epochs_noresp, return_log=True)\n",
    "\n",
    "epochs_dropped = epochs_ar.copy()\n",
    "reject_criteria = dict(eeg=100e-6) # 100 µV\n",
    "epochs_dropped.drop_bad(reject=reject_criteria)\n",
    "epochs_dropped.plot_drop_log();\n",
    "\n",
    "print('-------------------------------------')\n",
    "print(len(epochs_dropped), \"epochs left out of\", len(epochs_noresp), \"epochs\")\n",
    "print('-------------------------------------')\n",
    "\n",
    "#R0809: dropped 116 (10.2%) 1035\n",
    "#R0810\n",
    "#R0811: dropped 31 epochs (2.6%) 1152\n",
    "\n",
    "# Rereference\n",
    "\n",
    "epochs_ref = epochs_dropped.copy(); #epochs_ar.copy()\n",
    "epochs_ref.add_reference_channels('TP9');\n",
    "epochs_ref.set_eeg_reference(['TP9', 'TP10']);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3aa93b9-bb9e-408d-90a9-93fc5bc63de5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "epochs_ref.save(\"G:\\내 드라이브\\Research\\Decoding\\EEGViz\\R0809_epo.fif\", overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad31c953-09e6-4c1a-b633-3ad1c272b515",
   "metadata": {},
   "source": [
    "# Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4df4e5a-007f-4aa3-86e7-32d9a5c0e639",
   "metadata": {},
   "outputs": [],
   "source": [
    "import mne\n",
    "import re\n",
    "epochs_ref = mne.read_epochs(\"G:\\내 드라이브\\Research\\Decoding\\EEGViz\\R0809_auto_epo.fif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42547047-8602-432c-845e-f5255489cae2",
   "metadata": {},
   "outputs": [],
   "source": [
    "Eng = pd.read_excel(\"G:\\My Drive\\Research\\Decoding\\PsychopyViz\\stimuli_Engv.xlsx\")\n",
    "Kor = pd.read_excel(\"G:\\My Drive\\Research\\Decoding\\PsychopyViz\\stimuli_Korv.xlsx\")\n",
    "\n",
    "stim = []\n",
    "for i in range(64):\n",
    "    stim.append(Eng[\"word\"][i])\n",
    "    \n",
    "for i in range(64):\n",
    "    stim.append(Kor[\"word\"][i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39f11454-24dc-4ffe-a4c8-c215a5a878b4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "r1 = re.compile(\"Kor.*goat\")\n",
    "r2 = re.compile(\"Kor.*duck\")\n",
    "r3 = re.compile(\"Kor.*swan\")\n",
    "r4 = re.compile(\"Kor.*lion\")\n",
    "\n",
    "l1 = list(filter(r1.match, stim))\n",
    "l2 = list(filter(r2.match, stim))\n",
    "l3 = list(filter(r3.match, stim))\n",
    "l4 = list(filter(r4.match, stim))\n",
    "\n",
    "evoked1 = epochs_dropped[l1].average()\n",
    "evoked2 = epochs_dropped[l2].average()\n",
    "evoked3 = epochs_dropped[l3].average()\n",
    "evoked4 = epochs_dropped[l4].average()\n",
    "evoked1.comment = 'Korean goat'\n",
    "evoked2.comment = 'Korean duck'\n",
    "evoked3.comment = 'Korean swan'\n",
    "evoked4.comment = 'Korean lion'\n",
    "\n",
    "mne.viz.plot_compare_evokeds([evoked1, evoked2, evoked3, evoked4], picks='Cz', ylim=dict(eeg=[-10, 10]));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51bde49e-2b65-42da-90b0-ef112ff57c1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "r1 = re.compile(\"Eng.*goat\")\n",
    "r2 = re.compile(\"Eng.*duck\")\n",
    "r3 = re.compile(\"Eng.*swan\")\n",
    "r4 = re.compile(\"Eng.*lion\")\n",
    "\n",
    "l1 = list(filter(r1.match, stim))\n",
    "l2 = list(filter(r2.match, stim))\n",
    "l3 = list(filter(r3.match, stim))\n",
    "l4 = list(filter(r4.match, stim))\n",
    "\n",
    "evoked1 = epochs_dropped[l1].average()\n",
    "evoked2 = epochs_dropped[l2].average()\n",
    "evoked3 = epochs_dropped[l3].average()\n",
    "evoked4 = epochs_dropped[l4].average()\n",
    "evoked1.comment = 'English goat'\n",
    "evoked2.comment = 'English duck'\n",
    "evoked3.comment = 'English swan'\n",
    "evoked4.comment = 'English lion'\n",
    "\n",
    "mne.viz.plot_compare_evokeds([evoked1, evoked2, evoked3, evoked4], picks='Cz', ylim=dict(eeg=[-10, 10]));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff569b7d-8403-4ff1-aef2-5718febec027",
   "metadata": {},
   "outputs": [],
   "source": [
    "r1 = re.compile(\"Eng.*swan\")\n",
    "l1 = list(filter(r1.match, stim))\n",
    "epochs_ref[l1].plot()\n",
    "\n",
    "#how to reject the third epoch?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2da4923d-4c34-4693-8790-5c6a5863d1c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "r1 = re.compile(\"Eng.*lion\")\n",
    "l1 = list(filter(r1.match, stim))\n",
    "epochs_ref[l1].plot_image(picks='Cz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b2a940b-504e-4633-a9e3-e2ebd04b00f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "r1 = re.compile(\"Eng.*swan\")\n",
    "l1 = list(filter(r1.match, stim))\n",
    "epochs_dropped[l1].plot_image(picks='Cz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62c7576f-d9d5-43a1-bae1-afee12b217df",
   "metadata": {},
   "outputs": [],
   "source": [
    "r1 = re.compile(\"Eng.*swan\")\n",
    "l1 = list(filter(r1.match, stim))\n",
    "epochs_ref[l1].plot_image(picks='Cz')"
   ]
  }
 ],
 "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.15"
  },
  "toc-autonumbering": true,
  "toc-showtags": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}