From 2f8b5a8101a3b03f55884bead7674d5d79b2427f Mon Sep 17 00:00:00 2001 From: aliciajcampbell Date: Tue, 20 Aug 2024 12:08:44 +0200 Subject: [PATCH] update --- .DS_Store | Bin 18436 -> 18436 bytes ...LEISURE_rsEEG_aperiodic+theta+memory.ipynb | 0 .../Anijarv_OKTOS_aoEEG_erp_analysis.ipynb | 0 .../Anijarv_OKTOS_rsEEG_classic_bp.ipynb | 0 ...tEEG_aperiodic+alpha+ERP_prepoststim.ipynb | 867 ------------------ ..._sartEEG_aperiodic+theta+ERP_ind_chs.ipynb | 417 --------- ...sartEEG_aperiodic+theta+ERP_poststim.ipynb | 867 ------------------ ...althy_Ageing_Cross_Sectional_LEISURE.ipynb | 0 .../Hermens_LABS_rsEEG_alpha_reactivity.ipynb | 0 .../LABS_EEG_processing_pipeline.ipynb | 0 .../Mills_LABS_rsEEG_classic_bp.ipynb | 0 .../Mitchell_OKTOS_rsEEG_MSE+LZC.ipynb | 0 ...periodic+classic_bp+alpha_reactivity.ipynb | 0 .../Pace_METACOG_erp_analysis.ipynb | 0 .../mixedeffectsmodel_OKTOS_temp.ipynb | 0 15 files changed, 2151 deletions(-) rename "studies/Anij\303\244rv/Anijarv_LABS-LEISURE_rsEEG_aperiodic+theta+memory.ipynb" => studies/Anijarv_LABS-LEISURE_rsEEG_aperiodic+theta+memory.ipynb (100%) rename "studies/Anij\303\244rv/Anijarv_OKTOS_aoEEG_erp_analysis.ipynb" => studies/Anijarv_OKTOS_aoEEG_erp_analysis.ipynb (100%) rename "studies/Anij\303\244rv/Anijarv_OKTOS_rsEEG_classic_bp.ipynb" => studies/Anijarv_OKTOS_rsEEG_classic_bp.ipynb (100%) delete mode 100644 studies/Campbell/SART/Campbell_LEISURE_sartEEG_aperiodic+alpha+ERP_prepoststim.ipynb delete mode 100644 studies/Campbell/SART/Campbell_LEISURE_sartEEG_aperiodic+theta+ERP_ind_chs.ipynb delete mode 100644 studies/Campbell/SART/Campbell_LEISURE_sartEEG_aperiodic+theta+ERP_poststim.ipynb rename studies/{Campbell/Resting-state => }/Campbell_Resting_EEG_Sustained_Attention_Healthy_Ageing_Cross_Sectional_LEISURE.ipynb (100%) rename "studies/Anij\303\244rv/Hermens_LABS_rsEEG_alpha_reactivity.ipynb" => studies/Hermens_LABS_rsEEG_alpha_reactivity.ipynb (100%) rename "studies/Anij\303\244rv/LABS_EEG_processing_pipeline.ipynb" => studies/LABS_EEG_processing_pipeline.ipynb (100%) rename "studies/Anij\303\244rv/Mills_LABS_rsEEG_classic_bp.ipynb" => studies/Mills_LABS_rsEEG_classic_bp.ipynb (100%) rename "studies/Anij\303\244rv/Mitchell_OKTOS_rsEEG_MSE+LZC.ipynb" => studies/Mitchell_OKTOS_rsEEG_MSE+LZC.ipynb (100%) rename "studies/Anij\303\244rv/Pace_LEISURE_rsEEG_aperiodic+classic_bp+alpha_reactivity.ipynb" => studies/Pace_LEISURE_rsEEG_aperiodic+classic_bp+alpha_reactivity.ipynb (100%) rename "studies/Anij\303\244rv/Pace_METACOG_erp_analysis.ipynb" => studies/Pace_METACOG_erp_analysis.ipynb (100%) rename "studies/Anij\303\244rv/mixedeffectsmodel_OKTOS_temp.ipynb" => studies/mixedeffectsmodel_OKTOS_temp.ipynb (100%) diff --git a/.DS_Store b/.DS_Store index 75c6f0fbfd10379f27041c33aa45390ac5a2b0f3..3fe35c7df633e02e4834918457cc531d386e816a 100644 GIT binary patch delta 1685 zcmeH_TWl0%6vxkhx|I2((2v^7>=v?x7NP4FN@16+wuDPtE&{7`d(~KHcRpmoc6QmB zc9p7YwEDs&HA93K67WG2FAWJ8H6g|a9{}OSY>bJC@?gA0V|>yl&ULvtmGZgzc_u+UT~$ z8x!MZCT9(sC&UzQm^DQ)nRm>bWsi!2EuV}n)-{E~?Hjs!w)O2C+;_0Nd__HH4Xn|f zYv?HP*I3HZr$w2ss9)aPqG{tZY4iAu={(8XRn=2lS841KL$*AaGt48VF!IIQ3TB1-@ zhom~i>pi_vG-#MkZq4u6{c@Fw%6U5t!&{{Cm$anzkbCX!VDpX$d6&0C`;KluaxZo@ zmVy+e?KDJ_lp%|Xbc#;X8M;8P(jr}@_vr@Rq+4{GKBmv;OZtkwrf=vc`V|B~!7{W! zLlA4R9-WBbUTj4#?nfW`F@ys+gd`rs1axFzK)^x~CoqSn@B+@_Wt@xRJTBo)yoE)) zgB!Ssk5IfKz5B6d|hLOT3#-yWZ&(R~uA&&xV%-|V3i|6n><{NMtXRv^m;7XG(;1#^#HYLJ= dYYy30Vu|Z@rGF8qzr&OIKT!WEpaK=3egih5oD={6 delta 1669 zcmeH_Piz!b9LIm(>9)QZY4^2VrnK(DN(EuvLjUd-L2SiBOM%98%eDn2vz?JX+0K@o zX^9%M8tM&&ID&r?iSf_DgGQ}-Ffr=MYP@LnpouY2BOcVKH#H{uW@kZ3I2dE%#W~D- zzu$ZBm*0Dz&(awzox#hG*LZ_UemN9+v_no7^seqyIvMZl?Yl1-zqc!u+Lnx`dpf-j zl^$=RJKOeiZayy@x3lvTC^ezlx~rI~U9FF2w1eJnkzHOi74fvq&FZ9>DcIxU9?QF} zskyXjb=$i3Oo9fLz7`^)tWVcyuJY16{AY**YVtm4E8VY;(ZhBIT?rm&A@7lk8o zc`IWZb5>($R8`BoAzZJw=&g*mv&n^=IKCj+G;;xYquH=~py^{ZV~A zW2Hm>!mgqv|3O~cpkJ?YgJm4qVo)QYB~|XInPbAPp53Z-F>b08c_U}#gL1uEit$l3 zyJ!_1QLuyLHf=j=L{gocDa_5u3pqi0r#`^=3cJHJO0KofaNMzZ%f8sxPW@KKrkz=% zXpNn)t>VQa!;Br_;Y!{=GGq$3P+g4b`xu+@Z!&0_vX1|X2gY>$2hG% z;P8(bW9|0ZuVNXm;ccwoeSCn=@HH;r8~lV{m3k$taHUacQm#?r%C$k$ z3aZU#t7COg65i?V)$8YT`9h9JIesNQW4rfgBY)Mi?UJ66%c~K*v6N7>m>1SJmh_f2 zt?S}ijMqvHSlJJ++aSX_sX!2JPQ+qdk?IBEmaQ?3*GR9_%I2Pg3{efz{hG!@(lE8$ z+OkI1Zqja6eX%)J_72`6EBU_s39dvH8UF^jSLjW8hgPU8!CjzlC9dDW(13`FsQPguwpVcNk+Ak+6d4?pz@{X!8CW=qSqK#2VnJdXD&c_&wg>SrPU8`HIE(Xm w8qeTaJcsA;l8-QWz)N`cUjCqXYq_s~WmW0l)af7Cl=xpK`mc40TuG-t0Q1*pG5`Po diff --git "a/studies/Anij\303\244rv/Anijarv_LABS-LEISURE_rsEEG_aperiodic+theta+memory.ipynb" b/studies/Anijarv_LABS-LEISURE_rsEEG_aperiodic+theta+memory.ipynb similarity index 100% rename from "studies/Anij\303\244rv/Anijarv_LABS-LEISURE_rsEEG_aperiodic+theta+memory.ipynb" rename to studies/Anijarv_LABS-LEISURE_rsEEG_aperiodic+theta+memory.ipynb diff --git "a/studies/Anij\303\244rv/Anijarv_OKTOS_aoEEG_erp_analysis.ipynb" b/studies/Anijarv_OKTOS_aoEEG_erp_analysis.ipynb similarity index 100% rename from "studies/Anij\303\244rv/Anijarv_OKTOS_aoEEG_erp_analysis.ipynb" rename to studies/Anijarv_OKTOS_aoEEG_erp_analysis.ipynb diff --git "a/studies/Anij\303\244rv/Anijarv_OKTOS_rsEEG_classic_bp.ipynb" b/studies/Anijarv_OKTOS_rsEEG_classic_bp.ipynb similarity index 100% rename from "studies/Anij\303\244rv/Anijarv_OKTOS_rsEEG_classic_bp.ipynb" rename to studies/Anijarv_OKTOS_rsEEG_classic_bp.ipynb diff --git a/studies/Campbell/SART/Campbell_LEISURE_sartEEG_aperiodic+alpha+ERP_prepoststim.ipynb b/studies/Campbell/SART/Campbell_LEISURE_sartEEG_aperiodic+alpha+ERP_prepoststim.ipynb deleted file mode 100644 index ea60f17..0000000 --- a/studies/Campbell/SART/Campbell_LEISURE_sartEEG_aperiodic+alpha+ERP_prepoststim.ipynb +++ /dev/null @@ -1,867 +0,0 @@ -{ - "cells": [ - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## [Working Title]\n", - "\n", - "Created by Toomas Erik Anijärv in 25.05.2023\n", - "\n", - "This notebook is a representation of EEG processing done for the publication with one of the participants as an example.\n", - "\n", - "You are free to use this or any other code from this repository for your own projects and publications. Citation or reference to the repository is not required, but would be much appreciated (see more on README.md)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Import packages\n", - "import mne, os\n", - "import numpy as np\n", - "import pandas as pd\n", - "import matplotlib.pyplot as plt\n", - "from matplotlib.ticker import MultipleLocator\n", - "import seaborn as sns\n", - "from autoreject import (get_rejection_threshold, AutoReject)\n", - "from fooof import FOOOF\n", - "from fooof.plts.spectra import plot_spectra, plot_spectra_shading\n", - "\n", - "# Set the default directory\n", - "os.chdir('/Users/aliciacampbell/Documents/GitHub/EEG-pyline')\n", - "mne.set_log_level('error')\n", - "\n", - "# Import functions\n", - "import basic.arrange_data as arrange\n", - "import signal_processing.pre_process as prep\n", - "import signal_processing.spectral_analysis as spectr\n", - "import signal_processing.erp_analysis as erpan" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Folder where to get the raw EEG files\n", - "raw_folder = 'Data/Raw/'\n", - "\n", - "# Folder where to export the clean epochs files\n", - "clean_folder = 'Data/Clean/'\n", - "\n", - "# Folder where to save the results and plots\n", - "results_folder = 'Results/'\n", - "\n", - "# Sub-folder for the experiment (i.e. timepoint or group)\n", - "exp_folder = 'LEISURE/T1/SART/'\n", - "exp_condition = 'SART_T1'" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### PRE-PROCESSING & TASK PERFORMANCE" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# EOG + mastoid channels and stimulus channel\n", - "eog_channels = ['EXG1', 'EXG2', 'EXG3', 'EXG4', 'EXG5', 'EXG6', 'EXG7', 'EXG8']\n", - "stimulus_channel = 'Status'\n", - "\n", - "# Parameters for filter design\n", - "filter_design = dict(l_freq=1, h_freq=30, filter_length='auto', method='fir',\n", - " l_trans_bandwidth='auto', h_trans_bandwidth='auto',\n", - " phase='zero', fir_window='hamming', fir_design='firwin')\n", - "\n", - "# Epoch time window from event/stimuli\n", - "tminmax = [-0.5, 0.8]\n", - "\n", - "# Baseline correction time window\n", - "baseline_correction = None\n", - "\n", - "# Event names with IDs for GO and NO-GO trials\n", - "event_dict = {'GO trial': 4,\n", - " 'NO-GO trial': 8,\n", - " 'Button press': 16}\n", - "\n", - "# Button press ID\n", - "button_id = 16" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get directories of raw EEG files and set export directory for clean files\n", - "dir_inprogress = os.path.join(raw_folder,exp_folder)\n", - "export_dir = os.path.join(clean_folder,exp_folder)\n", - "file_dirs, subject_names = arrange.read_files(dir_inprogress,'.bdf')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Loop through all the subjects' directories (EEG files directories)\n", - "df_success = pd.DataFrame()\n", - "for i in range(len(file_dirs)):\n", - " print('\\n{} ({} / {})'.format(subject_names[i], i+1, len(file_dirs)))\n", - " # Read in the raw EEG data\n", - " raw = mne.io.read_raw_bdf(file_dirs[i], infer_types=True, eog=eog_channels,\n", - " stim_channel=stimulus_channel)\n", - " try:\n", - " raw = raw.drop_channels([\"Erg1\",\n", - " \"B1\",\"B2\",\"B3\",\"B4\",\"B5\",\"B6\",\"B7\",\"B8\",\"B9\",\"B10\",\"B11\",\"B12\",\n", - " \"B13\",\"B14\",\"B15\",\"B16\",\"B17\",\"B18\",\"B19\",\"B20\",\"B21\",\"B22\",\"B23\",\"B24\",\n", - " \"B25\",\"B26\",\"B27\",\"B28\",\"B29\",\"B30\",\"B31\",\"B32\",\n", - " \"C1\",\"C2\",\"C3-1\",\"C4-1\",\"C5\",\"C6\",\"C7\",\"C8\",\"C9\",\"C10\",\"C11\",\"C12\",\n", - " \"C13\",\"C14\",\"C15\",\"C16\",\"C17\",\"C18\",\"C19\",\"C20\",\"C21\",\"C22\",\"C23\",\"C24\",\n", - " \"C25\",\"C26\",\"C27\",\"C28\",\"C29\",\"C30\",\"C31\",\"C32\",\n", - " \"D1\",\"D2\",\"D3\",\"D4\",\"D5\",\"D6\",\"D7\",\"D8\",\"D9\",\"D10\",\"D11\",\"D12\",\n", - " \"D13\",\"D14\",\"D15\",\"D16\",\"D17\",\"D18\",\"D19\",\"D20\",\"D21\",\"D22\",\"D23\",\"D24\",\n", - " \"D25\",\"D26\",\"D27\",\"D28\",\"D29\",\"D30\",\"D31\",\"D32\"]).rename_channels({'C3-0':'C3','C4-0':'C4'})\n", - " print('Wrong montage. Dropping many electrodes and renaming some to match the right one.')\n", - " except:\n", - " try:\n", - " raw = raw.drop_channels([\"Erg1\"])\n", - " print('Correct montage, but dropping respiratory data.')\n", - " except:\n", - " raw = raw\n", - " print('Correct montage.')\n", - "\n", - " # Set the right montage (Biosemi32) and set reference as average across all channels\n", - " raw = raw.set_montage(mne.channels.make_standard_montage('biosemi32')).load_data()\\\n", - " .set_eeg_reference(ref_channels='average', verbose=False)\n", - " \n", - " # Filter the signal with bandpass filter and remove EOG artefacts with SSP\n", - " filt = prep.filter_raw_data(raw, filter_design, line_remove=None, eog_channels=eog_channels,\n", - " plot_filt=False, savefig=False, verbose=False)\n", - " \n", - " # Find events from the filtered EEG data and name them\n", - " events = mne.find_events(filt, stim_channel=stimulus_channel, consecutive=False, output='onset')\n", - " \n", - " # Create an array and dataframe of successful GO (followed with button press) and NO-GO trials (no button press)\n", - " success_events_total = []\n", - " success_go = []\n", - " unsuccess_nogo = []\n", - " df_success_temp = pd.DataFrame(index=[subject_names[i]],\n", - " data={'Total GO' : np.sum(events[:, 2] == event_dict['GO trial']),\n", - " 'Total NO-GO' : np.sum(events[:, 2] == event_dict['NO-GO trial']),\n", - " 'Correct GO': 0, 'Correct NO-GO': 0,\n", - " 'Incorrect GO': 0, 'Incorrect NO-GO': 0})\n", - " # Go through all events to check if they were successful or not\n", - " for m in range(len(events)):\n", - " # If event is a GO, check for button press\n", - " if events[m][2] == event_dict['GO trial']:\n", - " if events[m+1][2] == button_id:\n", - " # If there is a button press -> success\n", - " success_events_total.append(events[m])\n", - " success_events_total.append(events[m+1])\n", - " success_go.append(events[m])\n", - " success_go.append(events[m+1])\n", - " df_success_temp['Correct GO'] += 1\n", - " else:\n", - " # If there is no button press -> fail\n", - " df_success_temp['Incorrect GO'] += 1\n", - " # If event is a NO-GO, check for no button press\n", - " if events[m][2] == event_dict['NO-GO trial']:\n", - " if events[m+1][2] != button_id:\n", - " # If there is no button press -> success\n", - " success_events_total.append(events[m])\n", - " df_success_temp['Correct NO-GO'] += 1\n", - " else:\n", - " # If there is a button press -> fail\n", - " unsuccess_nogo.append(events[m])\n", - " unsuccess_nogo.append(events[m+1])\n", - " df_success_temp['Incorrect NO-GO'] += 1\n", - " success_events_total = np.asarray(success_events_total)\n", - " success_go = np.asarray(success_go)\n", - " unsuccess_nogo = np.asarray(unsuccess_nogo)\n", - "\n", - " # Calculate response times to button press for correct GO and incorrect NO-GO\n", - " if len(success_go)!=0:\n", - " rt_go = np.diff(success_go[:, 0])[0::2]/raw.info['sfreq']\n", - " else:\n", - " rt_go = 0\n", - " if len(unsuccess_nogo)!=0:\n", - " rt_nogo = np.diff(unsuccess_nogo[:, 0])[0::2]/raw.info['sfreq']\n", - " else:\n", - " rt_nogo = 0\n", - "\n", - " # Calculate descriptives for these response times\n", - " df_success_temp['Average RT (Correct GO)'] = np.mean(rt_go)\n", - " df_success_temp['Average RT (Incorrect NO-GO)'] = np.mean(rt_nogo)\n", - " df_success_temp['SD RT (Correct GO)'] = np.std(rt_go)\n", - " df_success_temp['SD RT (Incorrect NO-GO)'] = np.std(rt_nogo)\n", - " df_success_temp['Median RT (Correct GO)'] = np.median(rt_go)\n", - " df_success_temp['Median RT (Incorrect NO-GO)'] = np.median(rt_nogo)\n", - " df_success_temp['Minimum RT (Correct GO)'] = np.min(rt_go)\n", - " df_success_temp['Minimum RT (Incorrect NO-GO)'] = np.min(rt_nogo)\n", - " df_success_temp['Maximum RT (Correct GO)'] = np.max(rt_go)\n", - " df_success_temp['Maximum RT (Incorrect NO-GO)'] = np.max(rt_nogo)\n", - " df_success_temp['RTs (Correct GO)'] = str(rt_go)\n", - " df_success_temp['RTs (Incorrect NO-GO)'] = str(rt_nogo)\n", - "\n", - " # Merge the participant dataframe with the master dataframe\n", - " df_success = pd.concat([df_success, df_success_temp])\n", - "\n", - " # Plot all the events\n", - " %matplotlib inline\n", - " fig, axs = plt.subplots(1, 1, figsize=(10, 5))\n", - " fig = mne.viz.plot_events(success_events_total, sfreq=filt.info['sfreq'],\n", - " first_samp=filt.first_samp, event_id=event_dict,\n", - " axes=axs, show=False)\n", - " fig.subplots_adjust(right=0.7)\n", - " axs.set_title('Successful events ({})'.format(subject_names[i]))\n", - " plt.show()\n", - "\n", - " # Create epochs time-locked to successful GO and NO-GO events (without including the button press events)\n", - " picks = mne.pick_types(filt.info, eeg=True, stim=False)\n", - " epochs = mne.Epochs(filt, success_events_total[success_events_total[:, 2] != 16], #event_id=event_dict,\n", - " tmin=tminmax[0], tmax=tminmax[1], baseline=baseline_correction,\n", - " picks=picks, preload=True)\n", - " \n", - " # Plot the epochs' GFP plot before artefact rejection\n", - " epochs.plot_image(title=\"GFP without AR ({})\".format(subject_names[i]))\n", - "\n", - " # Use AutoReject to repair and remove epochs which are artefactual\n", - " reject_criteria = get_rejection_threshold(epochs)\n", - " print('Dropping epochs with rejection threshold:',reject_criteria)\n", - " epochs.drop_bad(reject=reject_criteria)\n", - "\n", - " ar = AutoReject(thresh_method='random_search', random_state=1)\n", - " ar.fit(epochs)\n", - " epochs_ar, reject_log = ar.transform(epochs, return_log=True)\n", - " reject_log.plot('horizontal')\n", - "\n", - " # Plot the epochs' GFP after artefact rejection\n", - " epochs_ar.average().plot()\n", - " epochs_ar.plot_image(title=\"GFP with AR ({})\".format(subject_names[i]))\n", - "\n", - " # Display the final epochs object meta-data\n", - " display(epochs_ar)\n", - "\n", - " # Save the cleaned EEG file as .fif file\n", - " try:\n", - " os.makedirs(export_dir)\n", - " except FileExistsError:\n", - " pass\n", - " try:\n", - " mne.Epochs.save(epochs_ar, fname='{}/{}_clean-epo.fif'.format(export_dir,\n", - " subject_names[i]),\n", - " overwrite=True)\n", - " except FileExistsError:\n", - " pass\n", - " \n", - "# Save the dataframe for task performance\n", - "df_success.to_excel('{}/{}/{}_task_performance.xlsx'.format(results_folder,exp_folder,exp_condition))" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### SPECTRAL ANALYSIS: APERIODIC + ALPHA ACTIVITY" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Brain regions and their channels /// do for Fz, Cz, Pz\n", - "# region = {'Global' : ['Fp1', 'AF3', 'F3', 'FC1', 'Fp2', 'AF4', 'F4', 'FC2', 'Fz',\n", - "# 'F7', 'FC5', 'T7', 'C3', 'F8', 'FC6', 'T8', 'C4', 'Cz',\n", - "# 'CP5', 'P3', 'P7', 'CP6', 'P4', 'P8', 'CP1', 'CP2', 'Pz',\n", - "# 'PO3', 'PO4', 'O1', 'O2', 'Oz']}\n", - "\n", - "region = {'Parieto-occipital' : ['P3', 'P7', 'P4', 'P8', 'Pz',\n", - " 'PO3', 'PO4', 'O1', 'O2', 'Oz']}\n", - "\n", - "# region = {'Fronto-central' : ['F3', 'FC1', 'F4', 'FC2', 'Fz', 'C3', 'C4', 'Cz', 'CP5', 'CP6', 'CP1', 'CP2']}\n", - "\n", - "# region = {'Central' : ['C3', 'C4', 'Cz', 'CP5', 'CP6', 'CP1', 'CP2']}\n", - "# region = {'Midline' : ['Fz', 'Cz', 'Pz']}\n", - "\n", - "# Power spectra estimation parameters\n", - "psd_params = {'pre' : dict(method='welch', fminmax=[1, 30], window='hamming', window_duration=0.5,\n", - " window_overlap=0, zero_padding=3, tminmax=[-0.5, 0]),\n", - " 'post' : dict(method='welch', fminmax=[1, 30], window='hamming', window_duration=0.8,\n", - " window_overlap=0, zero_padding=3, tminmax=[0, 0.8])}\n", - "\n", - "# FOOOF (specparam) model parameters\n", - "fooof_params = dict(peak_width_limits=[1, 12], max_n_peaks=float('inf'), min_peak_height=0.225,\n", - " peak_threshold=2.0, aperiodic_mode='fixed')\n", - "\n", - "# Band power of interest\n", - "bands = {'Alpha' : [7, 14]}\n", - "\n", - "# Flattened spectra amplitude scale (linear, log)\n", - "flat_spectr_scale = 'linear'\n", - "\n", - "# Plot more information on the model fit plots or not; and save these plots or not\n", - "plot_rich = True\n", - "savefig = True\n", - "\n", - "# Event names (i.e. different stimuli) within the epochs\n", - "# event_list = ['GO trial', 'NO-GO trial']\n", - "event_list = {'4' : 'GO trial',\n", - " '8' : 'NO-GO trial'}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get directories of clean EEG files and set export directory\n", - "dir_inprogress = os.path.join(clean_folder, exp_folder)\n", - "file_dirs, subject_names = arrange.read_files(dir_inprogress, \"_clean-epo.fif\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Pre-create results folders and dataframe\n", - "arrange.create_results_folders(exp_folder=exp_folder, results_folder=results_folder, fooof=True)\n", - "df_ch = pd.DataFrame()\n", - "# Go through all the files (subjects) in the folder\n", - "for i in range(len(file_dirs)):\n", - " # Read the clean data from the disk\n", - " epochs = mne.read_epochs(fname='{}/{}_clean-epo.fif'.format(dir_inprogress, subject_names[i]),\n", - " verbose=False)\n", - "\n", - " # Loop through all different events\n", - " df_ch_ev_pre = pd.DataFrame()\n", - " df_ch_ev_post = pd.DataFrame()\n", - " df_ch_ev_post_diff = pd.DataFrame()\n", - " for ev, evname in event_list.items():\n", - " print('{} for {} ({}/{})'.format(ev, subject_names[i], i+1, len(file_dirs)))\n", - " reg = list(region.keys())[0]\n", - " chs = list(region.values())[0]\n", - "\n", - " # Choose only epochs from the current event\n", - " epochs_ev = epochs[ev]\n", - "\n", - " ### PRE-EVENT PSD ESTIMATION\n", - "\n", - " # Calculate Welch's power spectrum density (FFT) for the mean post-event\n", - " [psds_pre, freqs_pre] = spectr.calculate_psd(epochs_ev, subject_names[i], **psd_params['pre'], verbose=False, plot=False)\n", - "\n", - " # Average all epochs together for each channel and also for each region\n", - " psds_pre = psds_pre.mean(axis=(0))\n", - " df_psds_ch_pre = arrange.array_to_df(subject_names[i], epochs_ev, psds_pre).\\\n", - " reset_index().drop(columns='Subject')\n", - "\n", - " # Choose only channel of interest data\n", - " # psds_temp_pre = df_psds_ch_pre[ch].to_numpy() # edit if using single channel, maybe only\n", - " psds_temp_pre = df_psds_ch_pre[chs].mean(axis=1).to_numpy()\n", - "\n", - " ### POST-EVENT PSD ESTIMATION\n", - "\n", - " # Calculate Welch's power spectrum density (FFT) for the mean post-event\n", - " [psds_post, freqs_post] = spectr.calculate_psd(epochs_ev, subject_names[i], **psd_params['post'], verbose=False, plot=False)\n", - "\n", - " # Average all epochs together for each channel and also for each region\n", - " psds_post = psds_post.mean(axis=(0))\n", - " df_psds_ch_post = arrange.array_to_df(subject_names[i], epochs_ev, psds_post).\\\n", - " reset_index().drop(columns='Subject')\n", - "\n", - " # Choose only channel of interest data\n", - " psds_temp_post = df_psds_ch_post[chs].mean(axis=1).to_numpy() # edit if using single channel, maybe only\n", - "\n", - " ### ERP & POST-minus-ERP PSD ESTIMATIONS\n", - "\n", - " # Average the event epochs in time domain\n", - " evoked_ev = epochs_ev.average(picks=chs)\n", - " # fig = plt.figure()\n", - " # fig = evoked_ev.plot()\n", - " # plt.show()\n", - "\n", - " # Calculate Welch's power spectrum density (FFT) for the ERP\n", - " [psds_post_erp, freqs_post] = spectr.calculate_psd(evoked_ev, subject_names[i], **psd_params['post'], verbose=False, plot=False)\n", - "\n", - " # Calculate the post-minus-ERP PSD subtracting ERP PSD from post-event PSD\n", - " psds_post_diff_ch = psds_temp_post - psds_post_erp.mean(axis=0) # edit if using single channel, maybe only\n", - "\n", - " ### SPECPARAM\n", - "\n", - " # Fit the spectrums with FOOOF\n", - " fm_pre = FOOOF(**fooof_params, verbose=True)\n", - " fm_pre.fit(freqs_pre, psds_temp_pre, psd_params['pre']['fminmax'])\n", - " fm_post = FOOOF(**fooof_params, verbose=True)\n", - " fm_post.fit(freqs_post, psds_temp_post, psd_params['post']['fminmax'])\n", - " fm_post_diff = FOOOF(**fooof_params, verbose=True)\n", - " fm_post_diff.fit(freqs_post, psds_post_diff_ch, psd_params['post']['fminmax'])\n", - "\n", - " def convert_flat_spectr_amplitude(fm, flat_spectr_scale='linear'):\n", - " \"\"\"Log-linear conversion based on the chosen amplitude scale\"\"\"\n", - " if flat_spectr_scale == 'linear':\n", - " flatten_spectrum = 10 ** fm._spectrum_flat\n", - " flat_spectr_ylabel = 'Amplitude (uV\\u00b2/Hz)'\n", - " elif flat_spectr_scale == 'log':\n", - " flatten_spectrum = fm._spectrum_flat\n", - " flat_spectr_ylabel = 'Log-normalised amplitude'\n", - " return flatten_spectrum, flat_spectr_ylabel\n", - " \n", - " # Log-linear conversion based on the chosen amplitude scale\n", - " flatten_spectrum_pre, _ = convert_flat_spectr_amplitude(fm_pre, flat_spectr_scale)\n", - " flatten_spectrum_post, _ = convert_flat_spectr_amplitude(fm_post, flat_spectr_scale)\n", - " flatten_spectrum_post_diff, _ = convert_flat_spectr_amplitude(fm_post_diff, flat_spectr_scale)\n", - "\n", - " # Find individual alpha band parameters\n", - " cf, pw, bw, abs_bp, rel_bp = spectr.find_ind_band(flatten_spectrum_pre, freqs_pre, bands['Alpha'], bw_size=6)\n", - " cf, pw, bw, abs_bp, rel_bp = spectr.find_ind_band(flatten_spectrum_post, freqs_post, bands['Alpha'], bw_size=6)\n", - " cf, pw, bw, abs_bp, rel_bp = spectr.find_ind_band(flatten_spectrum_post_diff, freqs_post, bands['Alpha'], bw_size=6)\n", - "\n", - " ### PLOTTING\n", - "\n", - " def plot_fooof_spectra(fm, flat_spectr_scale, psd_params, bands=None, abs_bp=None, rel_bp=None):\n", - " # Set plot styles\n", - " data_kwargs = {'color' : 'black', 'linewidth' : 1.4, 'label' : 'Original'}\n", - " model_kwargs = {'color' : 'red', 'linewidth' : 1.4, 'alpha' : 0.75, 'label' : 'Full model'}\n", - " aperiodic_kwargs = {'color' : 'blue', 'linewidth' : 1.4, 'alpha' : 0.75,\n", - " 'linestyle' : 'dashed', 'label' : 'Aperiodic model'}\n", - " flat_kwargs = {'color' : 'black', 'linewidth' : 1.4}\n", - " hvline_kwargs = {'color' : 'blue', 'linewidth' : 1.0, 'linestyle' : 'dashed', 'alpha' : 0.75}\n", - "\n", - " # Log-linear conversion based on the chosen amplitude scale\n", - " flatten_spectrum, flat_spectr_ylabel = convert_flat_spectr_amplitude(fm, flat_spectr_scale)\n", - " \n", - " # Plot power spectrum model + aperiodic fit for MEAN POST-EVENT PSD\n", - " fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(10, 4), dpi=100)\n", - " plot_spectra(fm.freqs, fm.power_spectrum,\n", - " ax=axs[0], **data_kwargs)\n", - " plot_spectra(fm.freqs, fm.fooofed_spectrum_,\n", - " ax=axs[0], **model_kwargs)\n", - " plot_spectra(fm.freqs, fm._ap_fit,\n", - " ax=axs[0], **aperiodic_kwargs)\n", - " axs[0].set_xlim(psd_params['fminmax'])\n", - " axs[0].grid(linewidth=0.2)\n", - " axs[0].set_xlabel('Frequency (Hz)')\n", - " axs[0].set_ylabel('Log-normalised power (log$_{10}$[µV\\u00b2/Hz])')\n", - " axs[0].set_title('Spectrum model fit')\n", - " axs[0].legend()\n", - " \n", - " # Flattened spectrum plot (i.e., minus aperiodic fit)\n", - " plot_spectra(fm.freqs, flatten_spectrum,\n", - " ax=axs[1], **flat_kwargs)\n", - " if bands!=None:\n", - " axs[1].axvspan(bands[list(bands.keys())[0]][0], bands[list(bands.keys())[0]][1], facecolor='green', alpha=0.2)\n", - " axs[1].set_xlim(psd_params['fminmax'])\n", - " axs[1].grid(linewidth=0.2)\n", - " axs[1].set_xlabel('Frequency (Hz)')\n", - " axs[1].set_ylabel(flat_spectr_ylabel)\n", - " axs[1].set_title('Flattened spectrum')\n", - "\n", - " # If true, plot all the exported variables on the plots\n", - " if plot_rich == True:\n", - " axs[0].annotate('Error: ' + str(np.round(fm.get_params('error'), 4)) +\n", - " '\\nR\\u00b2: ' + str(np.round(fm.get_params('r_squared'), 4)),\n", - " (0.1, 0.16), xycoords='figure fraction', color='red', fontsize=8.5)\n", - " axs[0].annotate('Exponent: ' + str(np.round(fm.get_params('aperiodic_params','exponent'), 4)) +\n", - " '\\nOffset: ' + str(np.round(fm.get_params('aperiodic_params','offset'), 4)),\n", - " (0.19, 0.16), xycoords='figure fraction', color='blue', fontsize=8.5)\n", - " if abs_bp!=None and rel_bp!=None:\n", - " axs[1].annotate('Absolute BP: '+str(np.round(abs_bp, 4))+'\\nRelative BP: '+str(np.round(rel_bp, 4)),\n", - " (0.69, 0.16), xycoords='figure fraction', color='green', fontsize=8.5)\n", - " \n", - " return fig, axs\n", - " \n", - " fig, axs = plot_fooof_spectra(fm_pre, flat_spectr_scale, psd_params['pre'], \n", - " bands=bands, abs_bp=abs_bp_pre, rel_bp=rel_bp_pre)\n", - " plt.suptitle('Mean pre-event PSD at {} ({})'.format(reg, subject_names[i]))\n", - " plt.tight_layout()\n", - " if savefig == True:\n", - " plt.savefig(fname='{}/{}/FOOOF/{}_{}_{}_{}_mean_pre_event_PSD.png'.format(results_folder, exp_folder,\n", - " exp_condition, subject_names[i],\n", - " reg, evname), dpi=300)\n", - " # plt.show()\n", - "\n", - " fig, axs = plot_fooof_spectra(fm_post, flat_spectr_scale, psd_params['post'], \n", - " bands=bands, abs_bp=abs_bp_post, rel_bp=rel_bp_post)\n", - " plt.suptitle('Mean post-event PSD at {} ({})'.format(reg, subject_names[i]))\n", - " plt.tight_layout()\n", - " if savefig == True:\n", - " plt.savefig(fname='{}/{}/FOOOF/{}_{}_{}_{}_mean_post_event_PSD.png'.format(results_folder, exp_folder,\n", - " exp_condition, subject_names[i],\n", - " reg, evname), dpi=300)\n", - " # plt.show()\n", - "\n", - " fig, axs = plot_fooof_spectra(fm_post_diff, flat_spectr_scale, psd_params['post'], \n", - " bands=bands, abs_bp=abs_bp_post_diff, rel_bp=rel_bp_post_diff)\n", - " plt.suptitle('Post-minus-ERP PSD at {} ({})'.format(reg, subject_names[i]))\n", - " plt.tight_layout()\n", - " if savefig == True:\n", - " plt.savefig(fname='{}/{}/FOOOF/{}_{}_{}_{}_post_minus_erp_PSD.png'.format(results_folder, exp_folder,\n", - " exp_condition, subject_names[i],\n", - " reg, evname), dpi=300)\n", - " # plt.show()\n", - "\n", - " ### EXPORTING\n", - "\n", - " def ev_to_df(df_ch_ev, i, fm, ch, ev, subject_names, \n", - " bands, abs_bp=None, rel_bp=None, stimname='stimulus'):\n", - " if abs_bp == None: abs_bp = np.nan\n", - " if rel_bp == None: rel_bp = np.nan\n", - " df_ch_ev.loc[i, 'Exponent'] = fm.get_params('aperiodic_params','exponent')\n", - " df_ch_ev.loc[i, 'Offset'] = fm.get_params('aperiodic_params','offset')\n", - " df_ch_ev.loc[i, '{} absolute power'.format(list(bands.keys())[0])] = abs_bp\n", - " df_ch_ev.loc[i, '{} relative power'.format(list(bands.keys())[0])] = rel_bp\n", - " df_ch_ev.loc[i, 'R_2'] = fm.get_params('r_squared')\n", - " df_ch_ev.loc[i, 'Error'] = fm.get_params('error')\n", - " df_ch_ev['Channel'] = ch\n", - " df_ch_ev['Event'] = ev\n", - " df_ch_ev['Type'] = stimname\n", - " df_ch_ev['Subject'] = subject_names[i]\n", - " return df_ch_ev\n", - " \n", - " # Add model parameters to dataframe for mean post-event\n", - " df_ch_ev_pre = ev_to_df(df_ch_ev_pre, i, fm_pre, reg, ev, subject_names, \n", - " bands, abs_bp_pre, rel_bp_pre, stimname='Mean pre')\n", - " df_ch_ev_post = ev_to_df(df_ch_ev_post, i, fm_post, reg, ev, subject_names, \n", - " bands, abs_bp_post, rel_bp_post, stimname='Mean post')\n", - " df_ch_ev_post_diff = ev_to_df(df_ch_ev_post_diff, i, fm_post_diff, reg, ev, subject_names, \n", - " bands, abs_bp_post_diff, rel_bp_post_diff, stimname='Mean post-ERP')\n", - "\n", - " # Concatenate to master dataframe for mean post-event\n", - " df_ch = pd.concat([df_ch, df_ch_ev_pre])\n", - " df_ch = pd.concat([df_ch, df_ch_ev_post])\n", - " df_ch = pd.concat([df_ch, df_ch_ev_post_diff])\n", - " \n", - "\n", - "# Reorder the channels and reset index\n", - "df_ch = df_ch[['Subject', 'Channel', 'Type', 'Event', 'Exponent', 'Offset',\n", - " '{} absolute power'.format(list(bands.keys())[0]),\n", - " '{} relative power'.format(list(bands.keys())[0]),\n", - " 'R_2', 'Error']]\n", - "df_ch = df_ch.reset_index(drop=True)\n", - "\n", - "# Export results for post-event data\n", - "df_ch.to_excel('{}/{}/FOOOF/{}_{}_specparam.xlsx'.format(results_folder, exp_folder, exp_condition, reg))\n", - "display(df_ch)\n", - "\n", - "# Export results with R2 <0.9 to indentify model fits to check\n", - "R2_threshold_df = df_ch[df_ch['R_2'] < 0.9]\n", - "R2_threshold_df.to_excel('{}/{}/FOOOF/{}_{}_specparam_r2_threshold.xlsx'.format(results_folder, exp_folder, exp_condition, reg))" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### ERP DETECTION & IDENTIFICATION" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tminmax = [-200, 800]\n", - "\n", - "# Time windows for different ERP components\n", - "erp_wins = {'N1' : [40, 170, -1],\n", - " 'N2' : [180, 350, -1],\n", - " 'P2' : [100, 260, 1],\n", - " 'P3' : [270, 500, 1]}\n", - "\n", - "# Channel of interest: Fz, Cz, Pz\n", - "channel_picks = 'Pz' \n", - "\n", - "# Event names (i.e. different stimuli) within the epochs\n", - "# event_list = ['GO trial', 'NO-GO trial']\n", - "event_list = {4 : 'GO trial',\n", - " 8 : 'NO-GO trial'}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get directories of clean EEG files and set export directory\n", - "dir_inprogress = os.path.join(clean_folder, exp_folder)\n", - "file_dirs, subject_names = arrange.read_files(dir_inprogress, '_clean-epo.fif')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Loop through all the subjects' directories (EEG files directories)\n", - "df_erps = pd.DataFrame()\n", - "arrange.create_results_folders(exp_folder=exp_folder,results_folder=results_folder, erps=True)\n", - "for i in range(len(file_dirs)):\n", - " erp_wins_temp = erp_wins\n", - " # Read the clean data from the disk\n", - " epochs = mne.read_epochs(fname='{}/{}_clean-epo.fif'.format(dir_inprogress, subject_names[i]), verbose=False)\n", - " \n", - " # Apply baseline correction\n", - " epochs = epochs.apply_baseline(baseline=(None, 0))\n", - "\n", - " ### create loop here for going through GO and NO-GO's separately\n", - " for ev, evname in event_list.items():\n", - " print('{} for {} ({}/{})'.format(ev, subject_names[i], i, len(file_dirs)))\n", - " # Create an averaged evoked object from epochs\n", - " evoked_signal = epochs[ev].average(picks=channel_picks)\n", - "\n", - " # remove or add if save_evoked === truuuu\n", - " evoked_signal.save('{}/{}/ERP analysis/{}_{}_{}_evoked-ave.fif'.format(results_folder, exp_folder,\n", - " subject_names[i], channel_picks,\n", - " ev), overwrite=True)\n", - "\n", - " # Find all the peaks in the evoked signal\n", - " minpeak_times, minpeak_mags, maxpeak_times, maxpeak_mags = erpan.find_all_peaks(evoked_signal, epochs, \n", - " t_range=tminmax, thresh=None,\n", - " subject_name=subject_names[i],\n", - " verbose=False, plot=False)\n", - " \n", - " # Identify which peaks are which ERPs based on the pre-defined ERP time windows\n", - " erp_peaks, not_erp_peaks = erpan.identify_erps(evoked_signal, erp_wins_temp, minpeak_times, minpeak_mags,\n", - " maxpeak_times, maxpeak_mags, t_range=tminmax, subject_name=subject_names[i],\n", - " verbose=False, plot=True, savefig=False,\n", - " results_foldername=results_folder, exp_folder=exp_folder)\n", - "\n", - " # After visual inspection, it's possible to re-define the time windows to look for the peak\n", - " while input('Do you need to do any manual time window changes? (leave empty if \"no\")') != '':\n", - " print('Changing time window parameters for {}'.format(subject_names[i]))\n", - " new_time_win = [None, None, None]\n", - "\n", - " # Ask user for which ERP they want to change or add\n", - " erp_tochange = input('What ERP time window you want to change (e.g., N1)?')\n", - "\n", - " # Ask user what should be the minimum timepoint of the time window for that ERP\n", - " new_time_win[0] = int(input('Enter MIN time of the window in interest for {} (e.g., 50)'.format(erp_tochange)))\n", - "\n", - " # Ask user what should be the maximum timepoint of the time window for that ERP\n", - " new_time_win[1] = int(input('Enter MAX time of the window in interest for {} (e.g., 100)'.format(erp_tochange)))\n", - "\n", - " # Ask user whether this ERP should be a postitive (1) or negative (-1) peak\n", - " new_time_win[2] = int(input('Enter whether to look for MIN (-1) or MAX (1) voltage for {}'.format(erp_tochange)))\n", - "\n", - " # Change the temporary ERP time window parameters to the user inputted parameters\n", - " erp_wins_temp[erp_tochange] = new_time_win\n", - " print('Changing', erp_tochange, 'with new time window:', str(new_time_win))\n", - "\n", - " # Use these new parameters to find either minimum or maximum value in that range\n", - " try:\n", - " erp_peaks = erpan.find_minmax_erp(evoked_signal, erp_peaks, erp_tochange, new_time_win,\n", - " t_range=tminmax, subject_name=subject_names[i], verbose=False, plot=True,\n", - " savefig=False, results_foldername=results_folder, exp_folder=exp_folder)\n", - " except:\n", - " print('Something went wrong with manual ERP detection, try again.')\n", - "\n", - " # Add this/these new temporary ERP to the main dataframe\n", - " df_erps_temp = erpan.erp_dict_to_df(erp_peaks, erp_wins_temp, subject_names[i])\n", - " df_erps_temp['Event'] = ev\n", - " df_erps_temp['Channel'] = channel_picks\n", - " df_erps = pd.concat([df_erps, df_erps_temp])\n", - " print('ERPs have been found and added to the dataframe for {}'.format(subject_names[i]))\n", - " display(df_erps)\n", - "\n", - "# Calculate relative peak-to-peak amplitudes between the ERPs\n", - "print('Adding relative amplitudes for N1-P2, P2-N2, N2-P3')\n", - "df_erps['N1-P2 amplitude'] = df_erps['P2 amplitude'] - df_erps['N1 amplitude']\n", - "df_erps['P2-N2 amplitude'] = df_erps['N2 amplitude'] - df_erps['P2 amplitude']\n", - "df_erps['N2-P3 amplitude'] = df_erps['P3 amplitude'] - df_erps['N2 amplitude']\n", - "\n", - "# Export all the detected ERPs to an Excel spreadsheet\n", - "display(df_erps)\n", - "df_erps.to_excel('{}/{}/ERP analysis/{}_{}_grandaverage_erps.xlsx'.format(results_folder,exp_folder,exp_condition,channel_picks))" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### DATA VISUALISATION: ERPs" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Export the figure to results folder or not\n", - "savefig = True\n", - "\n", - "# Subjects which to not plot\n", - "exclude_subjects = [] # ['OKTOS_0019', 'OKTOS_0024', 'OKTOS_0033']\n", - "\n", - "# Channel of interest\n", - "ch = 'Pz'\n", - "\n", - "# Event names (i.e. different stimuli) within the epochs\n", - "event_list = ['GO trial', 'NO-GO trial']" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "sns.set_theme(context='notebook', font_scale=1.3,\n", - " style='whitegrid', palette='muted',\n", - " font='sans-serif')\n", - "\n", - "# Get directories of clean EEG files and exclude the pre-defined subjects\n", - "dir_inprogress = os.path.join(clean_folder, exp_folder)\n", - "file_dirs, subject_names = arrange.read_files(dir_inprogress, \"_clean-epo.fif\",\n", - " exclude_subjects=exclude_subjects)\n", - "\n", - "# Loop through all the subjects' directories (EEG files directories)\n", - "evoked_signal_go = [None]*len(file_dirs)\n", - "evoked_signal_nogo = [None]*len(file_dirs)\n", - "for i in range(len(file_dirs)):\n", - " # Read the clean data from the disk\n", - " epochs = mne.read_epochs(fname='{}/{}_clean-epo.fif'.format(dir_inprogress,\n", - " subject_names[i]),\n", - " verbose=False)\n", - " \n", - " # Create an averaged evoked object from epochs for both events\n", - " evoked_signal_go[i] = epochs['GO trial'].average(picks=ch)\n", - " evoked_signal_nogo[i] = epochs['NO-GO trial'].average(picks=ch)\n", - "\n", - "# Average all the averaged evoked objects, thereby creating a grand average signals\n", - "go_master_grand_evoked_data = mne.grand_average(evoked_signal_go).data[0]*1e6\n", - "go_master_grand_evoked_times = mne.grand_average(evoked_signal_go).times*1e3\n", - "nogo_master_grand_evoked_data = mne.grand_average(evoked_signal_nogo).data[0]*1e6\n", - "nogo_master_grand_evoked_times = mne.grand_average(evoked_signal_nogo).times*1e3\n", - "\n", - "# Plot all experiments' grand average signals on a single plot\n", - "fig, ax = plt.subplots(figsize=(6, 4), layout='tight', dpi=150)\n", - "ax.plot(go_master_grand_evoked_times, go_master_grand_evoked_data, linewidth=3)\n", - "ax.plot(nogo_master_grand_evoked_times, nogo_master_grand_evoked_data, linewidth=3)\n", - "ax.legend(event_list)\n", - "ax.set_title('Grand average of all participants at {}'.format(ch))\n", - "ax.set_xlim([-200, 1000])\n", - "ax.yaxis.set_major_locator(MultipleLocator(1))\n", - "ax.set_xlabel('Time (ms)')\n", - "ax.set_ylabel('Amplitude (µV)')\n", - "ax.grid(which='major', axis='y', alpha=0.2)\n", - "ax.grid(which='major', axis='x', alpha=0.7)\n", - "if savefig == True:\n", - " plt.savefig(fname='{}/{}/GRAND_erpfig_{}.png'.format(results_folder, exp_folder, ch),\n", - " dpi=300)\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "df_ch" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "df_pre_go = df_ch.loc[(df_ch['Type']=='Mean pre') & (df_ch['Event']=='4'), :]\n", - "df_post_erp_go = df_ch.loc[(df_ch['Type']=='Mean post-ERP') & (df_ch['Event']=='4'), :]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "df_post_erp_go" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "df_post_erp_go['Exponent'].to_numpy()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "df_pre_go['Exponent']" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "df_post_erp_go['Exponent'] - df_pre_go['Exponent']" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "df_ch[df_ch['Type']=='Mean pre']" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "base", - "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.12.0" - }, - "orig_nbformat": 4 - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/studies/Campbell/SART/Campbell_LEISURE_sartEEG_aperiodic+theta+ERP_ind_chs.ipynb b/studies/Campbell/SART/Campbell_LEISURE_sartEEG_aperiodic+theta+ERP_ind_chs.ipynb deleted file mode 100644 index 6981b8a..0000000 --- a/studies/Campbell/SART/Campbell_LEISURE_sartEEG_aperiodic+theta+ERP_ind_chs.ipynb +++ /dev/null @@ -1,417 +0,0 @@ -{ - "cells": [ - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## [Working Title]\n", - "\n", - "Created by Toomas Erik Anijärv in 07.04.2024\n", - "\n", - "This notebook is a representation of EEG processing done for the publication with one of the participants as an example.\n", - "\n", - "You are free to use this or any other code from this repository for your own projects and publications. Citation or reference to the repository is not required, but would be much appreciated (see more on README.md)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Import packages\n", - "import mne, os\n", - "import numpy as np\n", - "import pandas as pd\n", - "import matplotlib.pyplot as plt\n", - "from matplotlib.ticker import MultipleLocator\n", - "import matplotlib.ticker as ticker\n", - "import seaborn as sns\n", - "from autoreject import (get_rejection_threshold, AutoReject)\n", - "from fooof import FOOOF\n", - "from fooof.plts.spectra import plot_spectra, plot_spectra_shading\n", - "\n", - "# Set the default directory\n", - "os.chdir('/Users/aliciacampbell/Documents/GitHub/EEG-pyline')\n", - "mne.set_log_level('error')\n", - "\n", - "# Import functions\n", - "import basic.arrange_data as arrange\n", - "import signal_processing.pre_process as prep\n", - "import signal_processing.spectral_analysis as spectr\n", - "import signal_processing.erp_analysis as erpan" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Folder where to get the raw EEG files\n", - "raw_folder = 'Data/Raw/'\n", - "\n", - "# Folder where to export the clean epochs files\n", - "clean_folder = 'Data/Clean/'\n", - "\n", - "# Folder where to save the results and plots\n", - "results_folder = 'Results/'\n", - "\n", - "# Sub-folder for the experiment (i.e. timepoint or group)\n", - "exp_folder = 'LEISURE/T1/SART/'\n", - "exp_condition = 'SART_T1'" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### SPECTRAL ANALYSIS: APERIODIC + THETA ACTIVITY" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def convert_flat_spectr_amplitude(fm, flat_spectr_scale='linear'):\n", - " \"\"\"Log-linear conversion based on the chosen amplitude scale\"\"\"\n", - " if flat_spectr_scale == 'linear':\n", - " flatten_spectrum = 10 ** fm._spectrum_flat\n", - " flat_spectr_ylabel = 'Amplitude (uV\\u00b2/Hz)'\n", - " elif flat_spectr_scale == 'log':\n", - " flatten_spectrum = fm._spectrum_flat\n", - " flat_spectr_ylabel = 'Log-normalised amplitude'\n", - " return flatten_spectrum, flat_spectr_ylabel\n", - "\n", - "def plot_fooof_spectra(fm, flat_spectr_scale, psd_params, bands=None, abs_bp=None, rel_bp=None, plot_rich=True):\n", - " \"\"\"Plot fooof spectra, regular with model fit and flattened spectrum.\"\"\"\n", - " # Set plot styles\n", - " data_kwargs = {'color' : 'black', 'linewidth' : 1.4, 'label' : 'Original'}\n", - " model_kwargs = {'color' : 'red', 'linewidth' : 1.4, 'alpha' : 0.75, 'label' : 'Full model'}\n", - " aperiodic_kwargs = {'color' : 'blue', 'linewidth' : 1.4, 'alpha' : 0.75,\n", - " 'linestyle' : 'dashed', 'label' : 'Aperiodic model'}\n", - " flat_kwargs = {'color' : 'black', 'linewidth' : 1.4}\n", - " hvline_kwargs = {'color' : 'blue', 'linewidth' : 1.0, 'linestyle' : 'dashed', 'alpha' : 0.75}\n", - "\n", - " # Log-linear conversion based on the chosen amplitude scale\n", - " flatten_spectrum, flat_spectr_ylabel = convert_flat_spectr_amplitude(fm, flat_spectr_scale)\n", - " \n", - " # Plot power spectrum model + aperiodic fit for MEAN POST-EVENT PSD\n", - " fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(10, 4), dpi=100)\n", - " plot_spectra(fm.freqs, fm.power_spectrum,\n", - " ax=axs[0], **data_kwargs)\n", - " plot_spectra(fm.freqs, fm.fooofed_spectrum_,\n", - " ax=axs[0], **model_kwargs)\n", - " plot_spectra(fm.freqs, fm._ap_fit,\n", - " ax=axs[0], **aperiodic_kwargs)\n", - " axs[0].set_xlim(psd_params['fminmax'])\n", - " axs[0].grid(linewidth=0.2)\n", - " axs[0].set_xlabel('Frequency (Hz)')\n", - " axs[0].set_ylabel('Log-normalised power (log$_{10}$[µV\\u00b2/Hz])')\n", - " axs[0].set_title('Spectrum model fit')\n", - " axs[0].legend()\n", - " \n", - " # Flattened spectrum plot (i.e., minus aperiodic fit)\n", - " plot_spectra(fm.freqs, flatten_spectrum,\n", - " ax=axs[1], **flat_kwargs)\n", - " if bands!=None:\n", - " axs[1].axvspan(bands[list(bands.keys())[0]][0], bands[list(bands.keys())[0]][1], facecolor='green', alpha=0.2)\n", - " axs[1].set_xlim(psd_params['fminmax'])\n", - " axs[1].grid(linewidth=0.2)\n", - " axs[1].set_xlabel('Frequency (Hz)')\n", - " axs[1].set_ylabel(flat_spectr_ylabel)\n", - " axs[1].set_title('Flattened spectrum')\n", - "\n", - " # If true, plot all the exported variables on the plots\n", - " if plot_rich == True:\n", - " axs[0].annotate('Error: ' + str(np.round(fm.get_params('error'), 4)) +\n", - " '\\nR\\u00b2: ' + str(np.round(fm.get_params('r_squared'), 4)),\n", - " (0.1, 0.16), xycoords='figure fraction', color='red', fontsize=8.5)\n", - " axs[0].annotate('Exponent: ' + str(np.round(fm.get_params('aperiodic_params','exponent'), 4)) +\n", - " '\\nOffset: ' + str(np.round(fm.get_params('aperiodic_params','offset'), 4)),\n", - " (0.19, 0.16), xycoords='figure fraction', color='blue', fontsize=8.5)\n", - " if abs_bp!=None and rel_bp!=None:\n", - " axs[1].annotate('Absolute BP: '+str(np.round(abs_bp, 4))+'\\nRelative BP: '+str(np.round(rel_bp, 4)),\n", - " (0.69, 0.16), xycoords='figure fraction', color='green', fontsize=8.5)\n", - " \n", - " return fig, axs\n", - "\n", - "def ev_to_df(df_ch_ev, i, fm, ch, ev, subject_names, \n", - " bands=None, abs_bp=None, rel_bp=None, stimname='stimulus'):\n", - " \"\"\"Helper function for saving fooof model results to dataframe.\"\"\"\n", - " if abs_bp == None: abs_bp = np.nan\n", - " if rel_bp == None: rel_bp = np.nan\n", - " df_ch_ev.loc[i, 'Exponent'] = fm.get_params('aperiodic_params','exponent')\n", - " df_ch_ev.loc[i, 'Offset'] = fm.get_params('aperiodic_params','offset')\n", - " if bands != None:\n", - " df_ch_ev.loc[i, '{} absolute power'.format(list(bands.keys())[0])] = abs_bp\n", - " df_ch_ev.loc[i, '{} relative power'.format(list(bands.keys())[0])] = rel_bp\n", - " df_ch_ev.loc[i, 'R_2'] = fm.get_params('r_squared')\n", - " df_ch_ev.loc[i, 'Error'] = fm.get_params('error')\n", - " df_ch_ev['Channel'] = ch\n", - " df_ch_ev['Event'] = ev\n", - " df_ch_ev['Type'] = stimname\n", - " df_ch_ev['Subject'] = subject_names[i]\n", - " return df_ch_ev" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Brain regions and their channels\n", - "channels = ['AF3', 'F3', 'FC1', 'AF4', 'F4', 'FC2', 'Fz',\n", - " 'F7', 'FC5', 'T7', 'C3', 'F8', 'FC6', 'T8', 'C4', 'Cz',\n", - " 'CP5', 'P3', 'P7', 'CP6', 'P4', 'P8', 'CP1', 'CP2', 'Pz',\n", - " 'PO3', 'PO4', 'O1', 'O2', 'Oz'] # 'Fp1', 'Fp2'\n", - "\n", - "# Power spectra estimation parameters\n", - "psd_params = {'pre' : dict(method='welch', fminmax=[1, 30], window='hamming', window_duration=0.5,\n", - " window_overlap=0, zero_padding=3, tminmax=[-0.5, 0]),\n", - " 'post' : dict(method='welch', fminmax=[1, 30], window='hamming', window_duration=0.8,\n", - " window_overlap=0, zero_padding=3, tminmax=[0, 0.8])}\n", - "\n", - "# FOOOF (specparam) model parameters\n", - "fooof_params = dict(peak_width_limits=[1, 12], max_n_peaks=float('inf'), min_peak_height=0.225,\n", - " peak_threshold=2.0, aperiodic_mode='fixed')\n", - "\n", - "# Set plot styles\n", - "data_kwargs = {'color' : 'black', 'linewidth' : 1.4}\n", - "model_kwargs = {'color' : 'red', 'linewidth' : 1.4, 'alpha' : 0.75}\n", - "aperiodic_kwargs = {'color' : 'blue', 'linewidth' : 1.4, 'alpha' : 0.75,\n", - " 'linestyle' : 'dashed'}\n", - "\n", - "# Flattened spectra amplitude scale (linear, log)\n", - "flat_spectr_scale = 'linear'\n", - "\n", - "# Plot more information on the model fit plots or not; and save these plots or not\n", - "plot_rich = True\n", - "savefig = True\n", - "\n", - "# Event names (i.e. different stimuli) within the epochs\n", - "# event_list = ['GO trial', 'NO-GO trial']\n", - "event_list = {'4' : 'GO trial',\n", - " '8' : 'NO-GO trial'}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get directories of clean EEG files and set export directory\n", - "dir_inprogress = os.path.join(clean_folder, exp_folder)\n", - "file_dirs, subject_names = arrange.read_files(dir_inprogress, \"_clean-epo.fif\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Pre-create results folders and dataframe\n", - "arrange.create_results_folders(exp_folder=exp_folder, results_folder=results_folder, fooof=True)\n", - "df_ch = pd.DataFrame()\n", - "# Go through all the files (subjects) in the folder\n", - "for i in range(5):#len(file_dirs)):\n", - " # Read the clean data from the disk\n", - " epochs = mne.read_epochs(fname='{}/{}_clean-epo.fif'.format(dir_inprogress, subject_names[i]),\n", - " verbose=False)\n", - "\n", - " # Loop through all different events\n", - " df_ch_ev_pre = pd.DataFrame()\n", - " df_ch_ev_post = pd.DataFrame()\n", - " df_ch_ev_post_diff = pd.DataFrame()\n", - " for ev, evname in event_list.items():\n", - " print('{} for {} ({}/{})'.format(ev, subject_names[i], i+1, len(file_dirs)))\n", - "\n", - " # Choose only epochs from the current event\n", - " epochs_ev = epochs[ev]\n", - "\n", - " # Calculate Welch's power spectrum density (FFT) for the mean pre-event and mean post-event\n", - " psds_pre, freqs_pre = spectr.calculate_psd(epochs_ev, subject_names[i], **psd_params['pre'], verbose=False, plot=False)\n", - " psds_post, freqs_post = spectr.calculate_psd(epochs_ev, subject_names[i], **psd_params['post'], verbose=False, plot=False)\n", - "\n", - " # Create all-channels PSD dataframe for mean pre-event and mean post-event\n", - " df_psds_pre = arrange.array_to_df(subject_names[i], epochs_ev, psds_pre.mean(axis=(0))).\\\n", - " reset_index().drop(columns='Subject')\n", - " df_psds_post = arrange.array_to_df(subject_names[i], epochs_ev, psds_post.mean(axis=(0))).\\\n", - " reset_index().drop(columns='Subject')\n", - " \n", - " # Go through each individual channel\n", - " fig_pre, axs_pre = plt.subplots(4, 8, figsize=(20, 10))\n", - " fig_pre.suptitle(f'{subject_names[i]} - {evname} - pre')\n", - " plt.tight_layout(pad=5.0)\n", - " fig_post, axs_post = plt.subplots(4, 8, figsize=(20, 10))\n", - " fig_post.suptitle(f'{subject_names[i]} - {evname} - post')\n", - " plt.tight_layout(pad=5.0)\n", - " fig_post_diff, axs_post_diff = plt.subplots(4, 8, figsize=(20, 10))\n", - " fig_post_diff.suptitle(f'{subject_names[i]} - {evname} - post-ERP')\n", - " plt.tight_layout(pad=5.0)\n", - " for c, ch in enumerate(channels):\n", - " # Choose only channel of interest PSD data for mean pre-event and mean post-event\n", - " psds_pre_ch = df_psds_pre[ch].to_numpy()\n", - " psds_post_ch = df_psds_post[ch].to_numpy()\n", - "\n", - " # Average the event epochs in time domain for that channel to get evoked object\n", - " evoked_ev_ch = epochs_ev.average(picks=ch)\n", - "\n", - " # Calculate Welch's power spectrum density (FFT) for the ERP (only for single channel this time)\n", - " psds_post_erp_ch = spectr.calculate_psd(evoked_ev_ch, subject_names[i], **psd_params['post'], verbose=False, plot=False)[0][0]\n", - "\n", - " # Substract the ERP PSD from the post-stimulus PSD -> difference PSD\n", - " psds_post_diff_ch = psds_post_ch - psds_post_erp_ch\n", - "\n", - " ### SPECPARAM\n", - "\n", - " # Fit the spectrums with FOOOF\n", - " fm_pre = FOOOF(**fooof_params, verbose=True)\n", - " fm_pre.fit(freqs_pre, psds_pre_ch, psd_params['pre']['fminmax'])\n", - " fm_post = FOOOF(**fooof_params, verbose=True)\n", - " fm_post.fit(freqs_post, psds_post_ch, psd_params['post']['fminmax'])\n", - " fm_post_diff = FOOOF(**fooof_params, verbose=True)\n", - " fm_post_diff.fit(freqs_post, psds_post_diff_ch, psd_params['post']['fminmax'])\n", - " \n", - " # Log-linear conversion based on the chosen amplitude scale\n", - " flatten_spectrum_pre, _ = convert_flat_spectr_amplitude(fm_pre, flat_spectr_scale)\n", - " flatten_spectrum_post, _ = convert_flat_spectr_amplitude(fm_post, flat_spectr_scale)\n", - " flatten_spectrum_post_diff, _ = convert_flat_spectr_amplitude(fm_post_diff, flat_spectr_scale)\n", - "\n", - " ### PLOTTING\n", - "\n", - " def fooof_subplot(fm, ax, data_kwargs, model_kwargs, aperiodic_kwargs):\n", - " # Get fit parameters\n", - " exp = fm.get_params('aperiodic_params','exponent')\n", - " off = fm.get_params('aperiodic_params','offset')\n", - " R2 = fm.get_params('r_squared')\n", - " error = fm.get_params('error')\n", - "\n", - " # Set title with condition\n", - " if exp <= 0 or error > 0.075 or R2 < 0.8:\n", - " title_color = 'red'\n", - " title_weight = 'bold'\n", - " else:\n", - " title_color = 'black'\n", - " title_weight = 'normal'\n", - "\n", - " # Plot power spectrum model + aperiodic fit for mean pre\n", - " plot_spectra(fm.freqs, fm.power_spectrum,\n", - " ax=ax, **data_kwargs)\n", - " plot_spectra(fm.freqs, fm.fooofed_spectrum_,\n", - " ax=ax, **model_kwargs)\n", - " plot_spectra(fm.freqs, fm._ap_fit,\n", - " ax=ax, **aperiodic_kwargs)\n", - " ax.set_xlim(psd_params['pre']['fminmax'])\n", - " ax.grid(linewidth=0.2)\n", - " ax.set_xlabel('')\n", - " ax.set_ylabel('')\n", - " ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%.1f'))\n", - " ax.set_title(f'{ch}\\nexp={np.round(exp, 3)}, off={np.round(off, 3)}\\nerror={np.round(error, 3)}, R2={np.round(R2, 3)}',\n", - " fontdict={'color': title_color, 'weight': title_weight})\n", - " \n", - "\n", - " \n", - " row = c // 8\n", - " col = c % 8\n", - " fooof_subplot(fm_pre, axs_pre[row, col], data_kwargs, model_kwargs, aperiodic_kwargs)\n", - " fooof_subplot(fm_post, axs_post[row, col], data_kwargs, model_kwargs, aperiodic_kwargs)\n", - " fooof_subplot(fm_post_diff, axs_post_diff[row, col], data_kwargs, model_kwargs, aperiodic_kwargs)\n", - "\n", - " # plt.show()\n", - " # plt.show()\n", - " plt.show()\n", - "\n", - "\n", - "\n", - "# ### PLOTTING\n", - "\n", - "# fig, axs = plot_fooof_spectra(fm_pre, flat_spectr_scale, psd_params['pre'], plot_rich=plot_rich)\n", - "# plt.suptitle('Mean pre-event PSD at {} ({})'.format(reg, subject_names[i]))\n", - "# plt.tight_layout()\n", - "# if savefig == True:\n", - "# plt.savefig(fname='{}/{}/FOOOF/{}_{}_{}_mean_pre_event_PSD.png'.format(results_folder, exp_folder,\n", - "# exp_condition, subject_names[i],\n", - "# reg), dpi=300)\n", - "# plt.show()\n", - "\n", - "# fig, axs = plot_fooof_spectra(fm_post, flat_spectr_scale, psd_params['post'], plot_rich=plot_rich)\n", - "# plt.suptitle('Mean post-event PSD at {} ({})'.format(reg, subject_names[i]))\n", - "# plt.tight_layout()\n", - "# if savefig == True:\n", - "# plt.savefig(fname='{}/{}/FOOOF/{}_{}_{}_mean_post_event_PSD.png'.format(results_folder, exp_folder,\n", - "# exp_condition, subject_names[i],\n", - "# reg), dpi=300)\n", - "# plt.show()\n", - "\n", - "# fig, axs = plot_fooof_spectra(fm_post_diff, flat_spectr_scale, psd_params['post'], plot_rich=plot_rich)\n", - "# plt.suptitle('Post-minus-ERP PSD at {} ({})'.format(reg, subject_names[i]))\n", - "# plt.tight_layout()\n", - "# if savefig == True:\n", - "# plt.savefig(fname='{}/{}/FOOOF/{}_{}_{}_post_minus_erp_PSD.png'.format(results_folder, exp_folder,\n", - "# exp_condition, subject_names[i],\n", - "# reg), dpi=300)\n", - "# plt.show()\n", - "\n", - "# ### EXPORTING\n", - " \n", - "# # Add model parameters to dataframe for mean post-event\n", - "# df_ch_ev_pre = ev_to_df(df_ch_ev_pre, i, fm_pre, reg, ev, subject_names, stimname='Mean pre')\n", - "# df_ch_ev_post = ev_to_df(df_ch_ev_post, i, fm_post, reg, ev, subject_names, stimname='Mean post')\n", - "# df_ch_ev_post_diff = ev_to_df(df_ch_ev_post_diff, i, fm_post_diff, reg, ev, subject_names, stimname='Mean post-ERP')\n", - "\n", - "# # Concatenate to master dataframe for mean post-event\n", - "# df_ch = pd.concat([df_ch, df_ch_ev_pre])\n", - "# df_ch = pd.concat([df_ch, df_ch_ev_post])\n", - "# df_ch = pd.concat([df_ch, df_ch_ev_post_diff])\n", - " \n", - "\n", - "# # Reorder the channels and reset index\n", - "# df_ch = df_ch[['Subject', 'Channel', 'Type', 'Event', 'Exponent', 'Offset', 'R_2', 'Error']]\n", - "# df_ch = df_ch.reset_index(drop=True)\n", - "\n", - "# # Export results for post-event data\n", - "# df_ch.to_excel('{}/{}/FOOOF/{}_{}_specparam.xlsx'.format(results_folder, exp_folder, exp_condition, reg))\n", - "# display(df_ch)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "psds_post_ch.shape" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "base", - "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.12.0" - }, - "orig_nbformat": 4 - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/studies/Campbell/SART/Campbell_LEISURE_sartEEG_aperiodic+theta+ERP_poststim.ipynb b/studies/Campbell/SART/Campbell_LEISURE_sartEEG_aperiodic+theta+ERP_poststim.ipynb deleted file mode 100644 index 6bab615..0000000 --- a/studies/Campbell/SART/Campbell_LEISURE_sartEEG_aperiodic+theta+ERP_poststim.ipynb +++ /dev/null @@ -1,867 +0,0 @@ -{ - "cells": [ - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## [Working Title]\n", - "\n", - "Created by Toomas Erik Anijärv in 25.05.2023\n", - "\n", - "This notebook is a representation of EEG processing done for the publication with one of the participants as an example.\n", - "\n", - "You are free to use this or any other code from this repository for your own projects and publications. Citation or reference to the repository is not required, but would be much appreciated (see more on README.md)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Import packages\n", - "import mne, os\n", - "import numpy as np\n", - "import pandas as pd\n", - "import matplotlib.pyplot as plt\n", - "from matplotlib.ticker import MultipleLocator\n", - "import seaborn as sns\n", - "from autoreject import (get_rejection_threshold, AutoReject)\n", - "from fooof import FOOOF\n", - "from fooof.plts.spectra import plot_spectra, plot_spectra_shading\n", - "\n", - "# Set the default directory\n", - "os.chdir('/Users/aliciacampbell/Documents/GitHub/EEG-pyline')\n", - "mne.set_log_level('error')\n", - "\n", - "# Import functions\n", - "import basic.arrange_data as arrange\n", - "import signal_processing.pre_process as prep\n", - "import signal_processing.spectral_analysis as spectr\n", - "import signal_processing.erp_analysis as erpan" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Folder where to get the raw EEG files\n", - "raw_folder = 'Data/Raw/'\n", - "\n", - "# Folder where to export the clean epochs files\n", - "clean_folder = 'Data/Clean/'\n", - "\n", - "# Folder where to save the results and plots\n", - "results_folder = 'Results/'\n", - "\n", - "# Sub-folder for the experiment (i.e. timepoint or group)\n", - "exp_folder = 'LEISURE/T1/SART/'\n", - "exp_condition = 'SART_T1'" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### PRE-PROCESSING & TASK PERFORMANCE" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# EOG + mastoid channels and stimulus channel\n", - "eog_channels = ['EXG1', 'EXG2', 'EXG3', 'EXG4', 'EXG5', 'EXG6', 'EXG7', 'EXG8']\n", - "stimulus_channel = 'Status'\n", - "\n", - "# Parameters for filter design\n", - "filter_design = dict(l_freq=1, h_freq=30, filter_length='auto', method='fir',\n", - " l_trans_bandwidth='auto', h_trans_bandwidth='auto',\n", - " phase='zero', fir_window='hamming', fir_design='firwin')\n", - "\n", - "# Epoch time window from event/stimuli\n", - "tminmax = [-0.5, 0.8]\n", - "\n", - "# Baseline correction time window\n", - "baseline_correction = None\n", - "\n", - "# Event names with IDs for GO and NO-GO trials\n", - "event_dict = {'GO trial': 4,\n", - " 'NO-GO trial': 8,\n", - " 'Button press': 16}\n", - "\n", - "# Button press ID\n", - "button_id = 16" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get directories of raw EEG files and set export directory for clean files\n", - "dir_inprogress = os.path.join(raw_folder,exp_folder)\n", - "export_dir = os.path.join(clean_folder,exp_folder)\n", - "file_dirs, subject_names = arrange.read_files(dir_inprogress,'.bdf')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Loop through all the subjects' directories (EEG files directories)\n", - "df_success = pd.DataFrame()\n", - "for i in range(len(file_dirs)):\n", - " print('\\n{} ({} / {})'.format(subject_names[i], i+1, len(file_dirs)))\n", - " # Read in the raw EEG data\n", - " raw = mne.io.read_raw_bdf(file_dirs[i], infer_types=True, eog=eog_channels,\n", - " stim_channel=stimulus_channel)\n", - " try:\n", - " raw = raw.drop_channels([\"Erg1\",\n", - " \"B1\",\"B2\",\"B3\",\"B4\",\"B5\",\"B6\",\"B7\",\"B8\",\"B9\",\"B10\",\"B11\",\"B12\",\n", - " \"B13\",\"B14\",\"B15\",\"B16\",\"B17\",\"B18\",\"B19\",\"B20\",\"B21\",\"B22\",\"B23\",\"B24\",\n", - " \"B25\",\"B26\",\"B27\",\"B28\",\"B29\",\"B30\",\"B31\",\"B32\",\n", - " \"C1\",\"C2\",\"C3-1\",\"C4-1\",\"C5\",\"C6\",\"C7\",\"C8\",\"C9\",\"C10\",\"C11\",\"C12\",\n", - " \"C13\",\"C14\",\"C15\",\"C16\",\"C17\",\"C18\",\"C19\",\"C20\",\"C21\",\"C22\",\"C23\",\"C24\",\n", - " \"C25\",\"C26\",\"C27\",\"C28\",\"C29\",\"C30\",\"C31\",\"C32\",\n", - " \"D1\",\"D2\",\"D3\",\"D4\",\"D5\",\"D6\",\"D7\",\"D8\",\"D9\",\"D10\",\"D11\",\"D12\",\n", - " \"D13\",\"D14\",\"D15\",\"D16\",\"D17\",\"D18\",\"D19\",\"D20\",\"D21\",\"D22\",\"D23\",\"D24\",\n", - " \"D25\",\"D26\",\"D27\",\"D28\",\"D29\",\"D30\",\"D31\",\"D32\"]).rename_channels({'C3-0':'C3','C4-0':'C4'})\n", - " print('Wrong montage. Dropping many electrodes and renaming some to match the right one.')\n", - " except:\n", - " try:\n", - " raw = raw.drop_channels([\"Erg1\"])\n", - " print('Correct montage, but dropping respiratory data.')\n", - " except:\n", - " raw = raw\n", - " print('Correct montage.')\n", - "\n", - " # Set the right montage (Biosemi32) and set reference as average across all channels\n", - " raw = raw.set_montage(mne.channels.make_standard_montage('biosemi32')).load_data()\\\n", - " .set_eeg_reference(ref_channels='average', verbose=False)\n", - " \n", - " # Filter the signal with bandpass filter and remove EOG artefacts with SSP\n", - " filt = prep.filter_raw_data(raw, filter_design, line_remove=None, eog_channels=eog_channels,\n", - " plot_filt=False, savefig=False, verbose=False)\n", - " \n", - " # Find events from the filtered EEG data and name them\n", - " events = mne.find_events(filt, stim_channel=stimulus_channel, consecutive=False, output='onset')\n", - " \n", - " # Create an array and dataframe of successful GO (followed with button press) and NO-GO trials (no button press)\n", - " success_events_total = []\n", - " success_go = []\n", - " unsuccess_nogo = []\n", - " df_success_temp = pd.DataFrame(index=[subject_names[i]],\n", - " data={'Total GO' : np.sum(events[:, 2] == event_dict['GO trial']),\n", - " 'Total NO-GO' : np.sum(events[:, 2] == event_dict['NO-GO trial']),\n", - " 'Correct GO': 0, 'Correct NO-GO': 0,\n", - " 'Incorrect GO': 0, 'Incorrect NO-GO': 0})\n", - " # Go through all events to check if they were successful or not\n", - " for m in range(len(events)):\n", - " # If event is a GO, check for button press\n", - " if events[m][2] == event_dict['GO trial']:\n", - " if events[m+1][2] == button_id:\n", - " # If there is a button press -> success\n", - " success_events_total.append(events[m])\n", - " success_events_total.append(events[m+1])\n", - " success_go.append(events[m])\n", - " success_go.append(events[m+1])\n", - " df_success_temp['Correct GO'] += 1\n", - " else:\n", - " # If there is no button press -> fail\n", - " df_success_temp['Incorrect GO'] += 1\n", - " # If event is a NO-GO, check for no button press\n", - " if events[m][2] == event_dict['NO-GO trial']:\n", - " if events[m+1][2] != button_id:\n", - " # If there is no button press -> success\n", - " success_events_total.append(events[m])\n", - " df_success_temp['Correct NO-GO'] += 1\n", - " else:\n", - " # If there is a button press -> fail\n", - " unsuccess_nogo.append(events[m])\n", - " unsuccess_nogo.append(events[m+1])\n", - " df_success_temp['Incorrect NO-GO'] += 1\n", - " success_events_total = np.asarray(success_events_total)\n", - " success_go = np.asarray(success_go)\n", - " unsuccess_nogo = np.asarray(unsuccess_nogo)\n", - "\n", - " # Calculate response times to button press for correct GO and incorrect NO-GO\n", - " if len(success_go)!=0:\n", - " rt_go = np.diff(success_go[:, 0])[0::2]/raw.info['sfreq']\n", - " else:\n", - " rt_go = 0\n", - " if len(unsuccess_nogo)!=0:\n", - " rt_nogo = np.diff(unsuccess_nogo[:, 0])[0::2]/raw.info['sfreq']\n", - " else:\n", - " rt_nogo = 0\n", - "\n", - " # Calculate descriptives for these response times\n", - " df_success_temp['Average RT (Correct GO)'] = np.mean(rt_go)\n", - " df_success_temp['Average RT (Incorrect NO-GO)'] = np.mean(rt_nogo)\n", - " df_success_temp['SD RT (Correct GO)'] = np.std(rt_go)\n", - " df_success_temp['SD RT (Incorrect NO-GO)'] = np.std(rt_nogo)\n", - " df_success_temp['Median RT (Correct GO)'] = np.median(rt_go)\n", - " df_success_temp['Median RT (Incorrect NO-GO)'] = np.median(rt_nogo)\n", - " df_success_temp['Minimum RT (Correct GO)'] = np.min(rt_go)\n", - " df_success_temp['Minimum RT (Incorrect NO-GO)'] = np.min(rt_nogo)\n", - " df_success_temp['Maximum RT (Correct GO)'] = np.max(rt_go)\n", - " df_success_temp['Maximum RT (Incorrect NO-GO)'] = np.max(rt_nogo)\n", - " df_success_temp['RTs (Correct GO)'] = str(rt_go)\n", - " df_success_temp['RTs (Incorrect NO-GO)'] = str(rt_nogo)\n", - "\n", - " # Merge the participant dataframe with the master dataframe\n", - " df_success = pd.concat([df_success, df_success_temp])\n", - "\n", - " # Plot all the events\n", - " %matplotlib inline\n", - " fig, axs = plt.subplots(1, 1, figsize=(10, 5))\n", - " fig = mne.viz.plot_events(success_events_total, sfreq=filt.info['sfreq'],\n", - " first_samp=filt.first_samp, event_id=event_dict,\n", - " axes=axs, show=False)\n", - " fig.subplots_adjust(right=0.7)\n", - " axs.set_title('Successful events ({})'.format(subject_names[i]))\n", - " plt.show()\n", - "\n", - " # Create epochs time-locked to successful GO and NO-GO events (without including the button press events)\n", - " picks = mne.pick_types(filt.info, eeg=True, stim=False)\n", - " epochs = mne.Epochs(filt, success_events_total[success_events_total[:, 2] != 16], #event_id=event_dict,\n", - " tmin=tminmax[0], tmax=tminmax[1], baseline=baseline_correction,\n", - " picks=picks, preload=True)\n", - " \n", - " # Plot the epochs' GFP plot before artefact rejection\n", - " epochs.plot_image(title=\"GFP without AR ({})\".format(subject_names[i]))\n", - "\n", - " # Use AutoReject to repair and remove epochs which are artefactual\n", - " reject_criteria = get_rejection_threshold(epochs)\n", - " print('Dropping epochs with rejection threshold:',reject_criteria)\n", - " epochs.drop_bad(reject=reject_criteria)\n", - "\n", - " ar = AutoReject(thresh_method='random_search', random_state=1)\n", - " ar.fit(epochs)\n", - " epochs_ar, reject_log = ar.transform(epochs, return_log=True)\n", - " reject_log.plot('horizontal')\n", - "\n", - " # Plot the epochs' GFP after artefact rejection\n", - " epochs_ar.average().plot()\n", - " epochs_ar.plot_image(title=\"GFP with AR ({})\".format(subject_names[i]))\n", - "\n", - " # Display the final epochs object meta-data\n", - " display(epochs_ar)\n", - "\n", - " # Save the cleaned EEG file as .fif file\n", - " try:\n", - " os.makedirs(export_dir)\n", - " except FileExistsError:\n", - " pass\n", - " try:\n", - " mne.Epochs.save(epochs_ar, fname='{}/{}_clean-epo.fif'.format(export_dir,\n", - " subject_names[i]),\n", - " overwrite=True)\n", - " except FileExistsError:\n", - " pass\n", - " \n", - "# Save the dataframe for task performance\n", - "df_success.to_excel('{}/{}/{}_task_performance.xlsx'.format(results_folder,exp_folder,exp_condition))" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### SPECTRAL ANALYSIS: APERIODIC + THETA ACTIVITY" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Brain regions and their channels /// do for Fz, Cz, Pz\n", - "# region = {'Global' : ['Fp1', 'AF3', 'F3', 'FC1', 'Fp2', 'AF4', 'F4', 'FC2', 'Fz',\n", - "# 'F7', 'FC5', 'T7', 'C3', 'F8', 'FC6', 'T8', 'C4', 'Cz',\n", - "# 'CP5', 'P3', 'P7', 'CP6', 'P4', 'P8', 'CP1', 'CP2', 'Pz',\n", - "# 'PO3', 'PO4', 'O1', 'O2', 'Oz']}\n", - "\n", - "region = {'Parieto-occipital' : ['P3', 'P7', 'P4', 'P8', 'Pz',\n", - " 'PO3', 'PO4', 'O1', 'O2', 'Oz']}\n", - "\n", - "# region = {'Fronto-central' : ['F3', 'FC1', 'F4', 'FC2', 'Fz', 'C3', 'C4', 'Cz', 'CP5', 'CP6', 'CP1', 'CP2']}\n", - "\n", - "# region = {'Central' : ['C3', 'C4', 'Cz', 'CP5', 'CP6', 'CP1', 'CP2']}\n", - "# region = {'Midline' : ['Fz', 'Cz', 'Pz']}\n", - "\n", - "# Power spectra estimation parameters\n", - "psd_params = {'pre' : dict(method='welch', fminmax=[1, 30], window='hamming', window_duration=0.5,\n", - " window_overlap=0, zero_padding=3, tminmax=[-0.5, 0]),\n", - " 'post' : dict(method='welch', fminmax=[1, 30], window='hamming', window_duration=0.8,\n", - " window_overlap=0, zero_padding=3, tminmax=[0, 0.8])}\n", - "\n", - "# FOOOF (specparam) model parameters\n", - "fooof_params = dict(peak_width_limits=[1, 12], max_n_peaks=float('inf'), min_peak_height=0.225,\n", - " peak_threshold=2.0, aperiodic_mode='fixed')\n", - "\n", - "# Band power of interest\n", - "bands = {'Theta' : [4, 8]}\n", - "\n", - "# Flattened spectra amplitude scale (linear, log)\n", - "flat_spectr_scale = 'linear'\n", - "\n", - "# Plot more information on the model fit plots or not; and save these plots or not\n", - "plot_rich = True\n", - "savefig = True\n", - "\n", - "# Event names (i.e. different stimuli) within the epochs\n", - "# event_list = ['GO trial', 'NO-GO trial']\n", - "event_list = {'4' : 'GO trial',\n", - " '8' : 'NO-GO trial'}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get directories of clean EEG files and set export directory\n", - "dir_inprogress = os.path.join(clean_folder, exp_folder)\n", - "file_dirs, subject_names = arrange.read_files(dir_inprogress, \"_clean-epo.fif\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Pre-create results folders and dataframe\n", - "arrange.create_results_folders(exp_folder=exp_folder, results_folder=results_folder, fooof=True)\n", - "df_ch = pd.DataFrame()\n", - "# Go through all the files (subjects) in the folder\n", - "for i in range(len(file_dirs)):\n", - " # Read the clean data from the disk\n", - " epochs = mne.read_epochs(fname='{}/{}_clean-epo.fif'.format(dir_inprogress, subject_names[i]),\n", - " verbose=False)\n", - "\n", - " # Loop through all different events\n", - " df_ch_ev_pre = pd.DataFrame()\n", - " df_ch_ev_post = pd.DataFrame()\n", - " df_ch_ev_post_diff = pd.DataFrame()\n", - " for ev, evname in event_list.items():\n", - " print('{} for {} ({}/{})'.format(ev, subject_names[i], i+1, len(file_dirs)))\n", - " reg = list(region.keys())[0]\n", - " chs = list(region.values())[0]\n", - "\n", - " # Choose only epochs from the current event\n", - " epochs_ev = epochs[ev]\n", - "\n", - " ### PRE-EVENT PSD ESTIMATION\n", - "\n", - " # Calculate Welch's power spectrum density (FFT) for the mean post-event\n", - " [psds_pre, freqs_pre] = spectr.calculate_psd(epochs_ev, subject_names[i], **psd_params['pre'], verbose=False, plot=False)\n", - "\n", - " # Average all epochs together for each channel and also for each region\n", - " psds_pre = psds_pre.mean(axis=(0))\n", - " df_psds_ch_pre = arrange.array_to_df(subject_names[i], epochs_ev, psds_pre).\\\n", - " reset_index().drop(columns='Subject')\n", - "\n", - " # Choose only channel of interest data\n", - " # psds_temp_pre = df_psds_ch_pre[ch].to_numpy() # edit if using single channel, maybe only\n", - " psds_temp_pre = df_psds_ch_pre[chs].mean(axis=1).to_numpy()\n", - "\n", - " ### POST-EVENT PSD ESTIMATION\n", - "\n", - " # Calculate Welch's power spectrum density (FFT) for the mean post-event\n", - " [psds_post, freqs_post] = spectr.calculate_psd(epochs_ev, subject_names[i], **psd_params['post'], verbose=False, plot=False)\n", - "\n", - " # Average all epochs together for each channel and also for each region\n", - " psds_post = psds_post.mean(axis=(0))\n", - " df_psds_ch_post = arrange.array_to_df(subject_names[i], epochs_ev, psds_post).\\\n", - " reset_index().drop(columns='Subject')\n", - "\n", - " # Choose only channel of interest data\n", - " psds_temp_post = df_psds_ch_post[chs].mean(axis=1).to_numpy() # edit if using single channel, maybe only\n", - "\n", - " ### ERP & POST-minus-ERP PSD ESTIMATIONS\n", - "\n", - " # Average the event epochs in time domain\n", - " evoked_ev = epochs_ev.average(picks=chs)\n", - " # fig = plt.figure()\n", - " # fig = evoked_ev.plot()\n", - " # plt.show()\n", - "\n", - " # Calculate Welch's power spectrum density (FFT) for the ERP\n", - " [psds_post_erp, freqs_post] = spectr.calculate_psd(evoked_ev, subject_names[i], **psd_params['post'], verbose=False, plot=False)\n", - "\n", - " # Calculate the post-minus-ERP PSD subtracting ERP PSD from post-event PSD\n", - " psds_post_diff_ch = psds_temp_post - psds_post_erp.mean(axis=0) # edit if using single channel, maybe only\n", - "\n", - " ### SPECPARAM\n", - "\n", - " # Fit the spectrums with FOOOF\n", - " fm_pre = FOOOF(**fooof_params, verbose=True)\n", - " fm_pre.fit(freqs_pre, psds_temp_pre, psd_params['pre']['fminmax'])\n", - " fm_post = FOOOF(**fooof_params, verbose=True)\n", - " fm_post.fit(freqs_post, psds_temp_post, psd_params['post']['fminmax'])\n", - " fm_post_diff = FOOOF(**fooof_params, verbose=True)\n", - " fm_post_diff.fit(freqs_post, psds_post_diff_ch, psd_params['post']['fminmax'])\n", - "\n", - " def convert_flat_spectr_amplitude(fm, flat_spectr_scale='linear'):\n", - " \"\"\"Log-linear conversion based on the chosen amplitude scale\"\"\"\n", - " if flat_spectr_scale == 'linear':\n", - " flatten_spectrum = 10 ** fm._spectrum_flat\n", - " flat_spectr_ylabel = 'Amplitude (uV\\u00b2/Hz)'\n", - " elif flat_spectr_scale == 'log':\n", - " flatten_spectrum = fm._spectrum_flat\n", - " flat_spectr_ylabel = 'Log-normalised amplitude'\n", - " return flatten_spectrum, flat_spectr_ylabel\n", - " \n", - " # Log-linear conversion based on the chosen amplitude scale\n", - " flatten_spectrum_pre, _ = convert_flat_spectr_amplitude(fm_pre, flat_spectr_scale)\n", - " flatten_spectrum_post, _ = convert_flat_spectr_amplitude(fm_post, flat_spectr_scale)\n", - " flatten_spectrum_post_diff, _ = convert_flat_spectr_amplitude(fm_post_diff, flat_spectr_scale)\n", - "\n", - " # Find individual theta band parameters\n", - " abs_bp_pre, rel_bp_pre = spectr.find_bp(flatten_spectrum_pre, freqs_pre, bands['Theta'])\n", - " abs_bp_post, rel_bp_post = spectr.find_bp(flatten_spectrum_post, freqs_post, bands['Theta'])\n", - " abs_bp_post_diff, rel_bp_post_diff = spectr.find_bp(flatten_spectrum_post_diff, freqs_post, bands['Theta'])\n", - "\n", - " ### PLOTTING\n", - "\n", - " def plot_fooof_spectra(fm, flat_spectr_scale, psd_params, bands=None, abs_bp=None, rel_bp=None):\n", - " # Set plot styles\n", - " data_kwargs = {'color' : 'black', 'linewidth' : 1.4, 'label' : 'Original'}\n", - " model_kwargs = {'color' : 'red', 'linewidth' : 1.4, 'alpha' : 0.75, 'label' : 'Full model'}\n", - " aperiodic_kwargs = {'color' : 'blue', 'linewidth' : 1.4, 'alpha' : 0.75,\n", - " 'linestyle' : 'dashed', 'label' : 'Aperiodic model'}\n", - " flat_kwargs = {'color' : 'black', 'linewidth' : 1.4}\n", - " hvline_kwargs = {'color' : 'blue', 'linewidth' : 1.0, 'linestyle' : 'dashed', 'alpha' : 0.75}\n", - "\n", - " # Log-linear conversion based on the chosen amplitude scale\n", - " flatten_spectrum, flat_spectr_ylabel = convert_flat_spectr_amplitude(fm, flat_spectr_scale)\n", - " \n", - " # Plot power spectrum model + aperiodic fit for MEAN POST-EVENT PSD\n", - " fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(10, 4), dpi=100)\n", - " plot_spectra(fm.freqs, fm.power_spectrum,\n", - " ax=axs[0], **data_kwargs)\n", - " plot_spectra(fm.freqs, fm.fooofed_spectrum_,\n", - " ax=axs[0], **model_kwargs)\n", - " plot_spectra(fm.freqs, fm._ap_fit,\n", - " ax=axs[0], **aperiodic_kwargs)\n", - " axs[0].set_xlim(psd_params['fminmax'])\n", - " axs[0].grid(linewidth=0.2)\n", - " axs[0].set_xlabel('Frequency (Hz)')\n", - " axs[0].set_ylabel('Log-normalised power (log$_{10}$[µV\\u00b2/Hz])')\n", - " axs[0].set_title('Spectrum model fit')\n", - " axs[0].legend()\n", - " \n", - " # Flattened spectrum plot (i.e., minus aperiodic fit)\n", - " plot_spectra(fm.freqs, flatten_spectrum,\n", - " ax=axs[1], **flat_kwargs)\n", - " if bands!=None:\n", - " axs[1].axvspan(bands[list(bands.keys())[0]][0], bands[list(bands.keys())[0]][1], facecolor='green', alpha=0.2)\n", - " axs[1].set_xlim(psd_params['fminmax'])\n", - " axs[1].grid(linewidth=0.2)\n", - " axs[1].set_xlabel('Frequency (Hz)')\n", - " axs[1].set_ylabel(flat_spectr_ylabel)\n", - " axs[1].set_title('Flattened spectrum')\n", - "\n", - " # If true, plot all the exported variables on the plots\n", - " if plot_rich == True:\n", - " axs[0].annotate('Error: ' + str(np.round(fm.get_params('error'), 4)) +\n", - " '\\nR\\u00b2: ' + str(np.round(fm.get_params('r_squared'), 4)),\n", - " (0.1, 0.16), xycoords='figure fraction', color='red', fontsize=8.5)\n", - " axs[0].annotate('Exponent: ' + str(np.round(fm.get_params('aperiodic_params','exponent'), 4)) +\n", - " '\\nOffset: ' + str(np.round(fm.get_params('aperiodic_params','offset'), 4)),\n", - " (0.19, 0.16), xycoords='figure fraction', color='blue', fontsize=8.5)\n", - " if abs_bp!=None and rel_bp!=None:\n", - " axs[1].annotate('Absolute BP: '+str(np.round(abs_bp, 4))+'\\nRelative BP: '+str(np.round(rel_bp, 4)),\n", - " (0.69, 0.16), xycoords='figure fraction', color='green', fontsize=8.5)\n", - " \n", - " return fig, axs\n", - " \n", - " fig, axs = plot_fooof_spectra(fm_pre, flat_spectr_scale, psd_params['pre'], \n", - " bands=bands, abs_bp=abs_bp_pre, rel_bp=rel_bp_pre)\n", - " plt.suptitle('Mean pre-event PSD at {} ({})'.format(reg, subject_names[i]))\n", - " plt.tight_layout()\n", - " if savefig == True:\n", - " plt.savefig(fname='{}/{}/FOOOF/{}_{}_{}_{}_mean_pre_event_PSD.png'.format(results_folder, exp_folder,\n", - " exp_condition, subject_names[i],\n", - " reg, evname), dpi=300)\n", - " # plt.show()\n", - "\n", - " fig, axs = plot_fooof_spectra(fm_post, flat_spectr_scale, psd_params['post'], \n", - " bands=bands, abs_bp=abs_bp_post, rel_bp=rel_bp_post)\n", - " plt.suptitle('Mean post-event PSD at {} ({})'.format(reg, subject_names[i]))\n", - " plt.tight_layout()\n", - " if savefig == True:\n", - " plt.savefig(fname='{}/{}/FOOOF/{}_{}_{}_{}_mean_post_event_PSD.png'.format(results_folder, exp_folder,\n", - " exp_condition, subject_names[i],\n", - " reg, evname), dpi=300)\n", - " # plt.show()\n", - "\n", - " fig, axs = plot_fooof_spectra(fm_post_diff, flat_spectr_scale, psd_params['post'], \n", - " bands=bands, abs_bp=abs_bp_post_diff, rel_bp=rel_bp_post_diff)\n", - " plt.suptitle('Post-minus-ERP PSD at {} ({})'.format(reg, subject_names[i]))\n", - " plt.tight_layout()\n", - " if savefig == True:\n", - " plt.savefig(fname='{}/{}/FOOOF/{}_{}_{}_{}_post_minus_erp_PSD.png'.format(results_folder, exp_folder,\n", - " exp_condition, subject_names[i],\n", - " reg, evname), dpi=300)\n", - " # plt.show()\n", - "\n", - " ### EXPORTING\n", - "\n", - " def ev_to_df(df_ch_ev, i, fm, ch, ev, subject_names, \n", - " bands, abs_bp=None, rel_bp=None, stimname='stimulus'):\n", - " if abs_bp == None: abs_bp = np.nan\n", - " if rel_bp == None: rel_bp = np.nan\n", - " df_ch_ev.loc[i, 'Exponent'] = fm.get_params('aperiodic_params','exponent')\n", - " df_ch_ev.loc[i, 'Offset'] = fm.get_params('aperiodic_params','offset')\n", - " df_ch_ev.loc[i, '{} absolute power'.format(list(bands.keys())[0])] = abs_bp\n", - " df_ch_ev.loc[i, '{} relative power'.format(list(bands.keys())[0])] = rel_bp\n", - " df_ch_ev.loc[i, 'R_2'] = fm.get_params('r_squared')\n", - " df_ch_ev.loc[i, 'Error'] = fm.get_params('error')\n", - " df_ch_ev['Channel'] = ch\n", - " df_ch_ev['Event'] = ev\n", - " df_ch_ev['Type'] = stimname\n", - " df_ch_ev['Subject'] = subject_names[i]\n", - " return df_ch_ev\n", - " \n", - " # Add model parameters to dataframe for mean post-event\n", - " df_ch_ev_pre = ev_to_df(df_ch_ev_pre, i, fm_pre, reg, ev, subject_names, \n", - " bands, abs_bp_pre, rel_bp_pre, stimname='Mean pre')\n", - " df_ch_ev_post = ev_to_df(df_ch_ev_post, i, fm_post, reg, ev, subject_names, \n", - " bands, abs_bp_post, rel_bp_post, stimname='Mean post')\n", - " df_ch_ev_post_diff = ev_to_df(df_ch_ev_post_diff, i, fm_post_diff, reg, ev, subject_names, \n", - " bands, abs_bp_post_diff, rel_bp_post_diff, stimname='Mean post-ERP')\n", - "\n", - " # Concatenate to master dataframe for mean post-event\n", - " df_ch = pd.concat([df_ch, df_ch_ev_pre])\n", - " df_ch = pd.concat([df_ch, df_ch_ev_post])\n", - " df_ch = pd.concat([df_ch, df_ch_ev_post_diff])\n", - " \n", - "\n", - "# Reorder the channels and reset index\n", - "df_ch = df_ch[['Subject', 'Channel', 'Type', 'Event', 'Exponent', 'Offset',\n", - " '{} absolute power'.format(list(bands.keys())[0]),\n", - " '{} relative power'.format(list(bands.keys())[0]),\n", - " 'R_2', 'Error']]\n", - "df_ch = df_ch.reset_index(drop=True)\n", - "\n", - "# Export results for post-event data\n", - "df_ch.to_excel('{}/{}/FOOOF/{}_{}_specparam.xlsx'.format(results_folder, exp_folder, exp_condition, reg))\n", - "display(df_ch)\n", - "\n", - "# Export results with R2 <0.9 to indentify model fits to check\n", - "R2_threshold_df = df_ch[df_ch['R_2'] < 0.9]\n", - "R2_threshold_df.to_excel('{}/{}/FOOOF/{}_{}_specparam_r2_threshold.xlsx'.format(results_folder, exp_folder, exp_condition, reg))" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### ERP DETECTION & IDENTIFICATION" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tminmax = [-200, 800]\n", - "\n", - "# Time windows for different ERP components\n", - "erp_wins = {'N1' : [40, 170, -1],\n", - " 'N2' : [180, 350, -1],\n", - " 'P2' : [100, 260, 1],\n", - " 'P3' : [270, 500, 1]}\n", - "\n", - "# Channel of interest: Fz, Cz, Pz\n", - "channel_picks = 'Pz' \n", - "\n", - "# Event names (i.e. different stimuli) within the epochs\n", - "# event_list = ['GO trial', 'NO-GO trial']\n", - "event_list = {4 : 'GO trial',\n", - " 8 : 'NO-GO trial'}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get directories of clean EEG files and set export directory\n", - "dir_inprogress = os.path.join(clean_folder, exp_folder)\n", - "file_dirs, subject_names = arrange.read_files(dir_inprogress, '_clean-epo.fif')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Loop through all the subjects' directories (EEG files directories)\n", - "df_erps = pd.DataFrame()\n", - "arrange.create_results_folders(exp_folder=exp_folder,results_folder=results_folder, erps=True)\n", - "for i in range(len(file_dirs)):\n", - " erp_wins_temp = erp_wins\n", - " # Read the clean data from the disk\n", - " epochs = mne.read_epochs(fname='{}/{}_clean-epo.fif'.format(dir_inprogress, subject_names[i]), verbose=False)\n", - " \n", - " # Apply baseline correction\n", - " epochs = epochs.apply_baseline(baseline=(None, 0))\n", - "\n", - " ### create loop here for going through GO and NO-GO's separately\n", - " for ev, evname in event_list.items():\n", - " print('{} for {} ({}/{})'.format(ev, subject_names[i], i, len(file_dirs)))\n", - " # Create an averaged evoked object from epochs\n", - " evoked_signal = epochs[ev].average(picks=channel_picks)\n", - "\n", - " # remove or add if save_evoked === truuuu\n", - " evoked_signal.save('{}/{}/ERP analysis/{}_{}_{}_evoked-ave.fif'.format(results_folder, exp_folder,\n", - " subject_names[i], channel_picks,\n", - " ev), overwrite=True)\n", - "\n", - " # Find all the peaks in the evoked signal\n", - " minpeak_times, minpeak_mags, maxpeak_times, maxpeak_mags = erpan.find_all_peaks(evoked_signal, epochs, \n", - " t_range=tminmax, thresh=None,\n", - " subject_name=subject_names[i],\n", - " verbose=False, plot=False)\n", - " \n", - " # Identify which peaks are which ERPs based on the pre-defined ERP time windows\n", - " erp_peaks, not_erp_peaks = erpan.identify_erps(evoked_signal, erp_wins_temp, minpeak_times, minpeak_mags,\n", - " maxpeak_times, maxpeak_mags, t_range=tminmax, subject_name=subject_names[i],\n", - " verbose=False, plot=True, savefig=False,\n", - " results_foldername=results_folder, exp_folder=exp_folder)\n", - "\n", - " # After visual inspection, it's possible to re-define the time windows to look for the peak\n", - " while input('Do you need to do any manual time window changes? (leave empty if \"no\")') != '':\n", - " print('Changing time window parameters for {}'.format(subject_names[i]))\n", - " new_time_win = [None, None, None]\n", - "\n", - " # Ask user for which ERP they want to change or add\n", - " erp_tochange = input('What ERP time window you want to change (e.g., N1)?')\n", - "\n", - " # Ask user what should be the minimum timepoint of the time window for that ERP\n", - " new_time_win[0] = int(input('Enter MIN time of the window in interest for {} (e.g., 50)'.format(erp_tochange)))\n", - "\n", - " # Ask user what should be the maximum timepoint of the time window for that ERP\n", - " new_time_win[1] = int(input('Enter MAX time of the window in interest for {} (e.g., 100)'.format(erp_tochange)))\n", - "\n", - " # Ask user whether this ERP should be a postitive (1) or negative (-1) peak\n", - " new_time_win[2] = int(input('Enter whether to look for MIN (-1) or MAX (1) voltage for {}'.format(erp_tochange)))\n", - "\n", - " # Change the temporary ERP time window parameters to the user inputted parameters\n", - " erp_wins_temp[erp_tochange] = new_time_win\n", - " print('Changing', erp_tochange, 'with new time window:', str(new_time_win))\n", - "\n", - " # Use these new parameters to find either minimum or maximum value in that range\n", - " try:\n", - " erp_peaks = erpan.find_minmax_erp(evoked_signal, erp_peaks, erp_tochange, new_time_win,\n", - " t_range=tminmax, subject_name=subject_names[i], verbose=False, plot=True,\n", - " savefig=False, results_foldername=results_folder, exp_folder=exp_folder)\n", - " except:\n", - " print('Something went wrong with manual ERP detection, try again.')\n", - "\n", - " # Add this/these new temporary ERP to the main dataframe\n", - " df_erps_temp = erpan.erp_dict_to_df(erp_peaks, erp_wins_temp, subject_names[i])\n", - " df_erps_temp['Event'] = ev\n", - " df_erps_temp['Channel'] = channel_picks\n", - " df_erps = pd.concat([df_erps, df_erps_temp])\n", - " print('ERPs have been found and added to the dataframe for {}'.format(subject_names[i]))\n", - " display(df_erps)\n", - "\n", - "# Calculate relative peak-to-peak amplitudes between the ERPs\n", - "print('Adding relative amplitudes for N1-P2, P2-N2, N2-P3')\n", - "df_erps['N1-P2 amplitude'] = df_erps['P2 amplitude'] - df_erps['N1 amplitude']\n", - "df_erps['P2-N2 amplitude'] = df_erps['N2 amplitude'] - df_erps['P2 amplitude']\n", - "df_erps['N2-P3 amplitude'] = df_erps['P3 amplitude'] - df_erps['N2 amplitude']\n", - "\n", - "# Export all the detected ERPs to an Excel spreadsheet\n", - "display(df_erps)\n", - "df_erps.to_excel('{}/{}/ERP analysis/{}_{}_grandaverage_erps.xlsx'.format(results_folder,exp_folder,exp_condition,channel_picks))" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### DATA VISUALISATION: ERPs" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Export the figure to results folder or not\n", - "savefig = True\n", - "\n", - "# Subjects which to not plot\n", - "exclude_subjects = [] # ['OKTOS_0019', 'OKTOS_0024', 'OKTOS_0033']\n", - "\n", - "# Channel of interest\n", - "ch = 'Pz'\n", - "\n", - "# Event names (i.e. different stimuli) within the epochs\n", - "event_list = ['GO trial', 'NO-GO trial']" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "sns.set_theme(context='notebook', font_scale=1.3,\n", - " style='whitegrid', palette='muted',\n", - " font='sans-serif')\n", - "\n", - "# Get directories of clean EEG files and exclude the pre-defined subjects\n", - "dir_inprogress = os.path.join(clean_folder, exp_folder)\n", - "file_dirs, subject_names = arrange.read_files(dir_inprogress, \"_clean-epo.fif\",\n", - " exclude_subjects=exclude_subjects)\n", - "\n", - "# Loop through all the subjects' directories (EEG files directories)\n", - "evoked_signal_go = [None]*len(file_dirs)\n", - "evoked_signal_nogo = [None]*len(file_dirs)\n", - "for i in range(len(file_dirs)):\n", - " # Read the clean data from the disk\n", - " epochs = mne.read_epochs(fname='{}/{}_clean-epo.fif'.format(dir_inprogress,\n", - " subject_names[i]),\n", - " verbose=False)\n", - " \n", - " # Create an averaged evoked object from epochs for both events\n", - " evoked_signal_go[i] = epochs['GO trial'].average(picks=ch)\n", - " evoked_signal_nogo[i] = epochs['NO-GO trial'].average(picks=ch)\n", - "\n", - "# Average all the averaged evoked objects, thereby creating a grand average signals\n", - "go_master_grand_evoked_data = mne.grand_average(evoked_signal_go).data[0]*1e6\n", - "go_master_grand_evoked_times = mne.grand_average(evoked_signal_go).times*1e3\n", - "nogo_master_grand_evoked_data = mne.grand_average(evoked_signal_nogo).data[0]*1e6\n", - "nogo_master_grand_evoked_times = mne.grand_average(evoked_signal_nogo).times*1e3\n", - "\n", - "# Plot all experiments' grand average signals on a single plot\n", - "fig, ax = plt.subplots(figsize=(6, 4), layout='tight', dpi=150)\n", - "ax.plot(go_master_grand_evoked_times, go_master_grand_evoked_data, linewidth=3)\n", - "ax.plot(nogo_master_grand_evoked_times, nogo_master_grand_evoked_data, linewidth=3)\n", - "ax.legend(event_list)\n", - "ax.set_title('Grand average of all participants at {}'.format(ch))\n", - "ax.set_xlim([-200, 1000])\n", - "ax.yaxis.set_major_locator(MultipleLocator(1))\n", - "ax.set_xlabel('Time (ms)')\n", - "ax.set_ylabel('Amplitude (µV)')\n", - "ax.grid(which='major', axis='y', alpha=0.2)\n", - "ax.grid(which='major', axis='x', alpha=0.7)\n", - "if savefig == True:\n", - " plt.savefig(fname='{}/{}/GRAND_erpfig_{}.png'.format(results_folder, exp_folder, ch),\n", - " dpi=300)\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "df_ch" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "df_pre_go = df_ch.loc[(df_ch['Type']=='Mean pre') & (df_ch['Event']=='4'), :]\n", - "df_post_erp_go = df_ch.loc[(df_ch['Type']=='Mean post-ERP') & (df_ch['Event']=='4'), :]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "df_post_erp_go" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "df_post_erp_go['Exponent'].to_numpy()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "df_pre_go['Exponent']" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "df_post_erp_go['Exponent'] - df_pre_go['Exponent']" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "df_ch[df_ch['Type']=='Mean pre']" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "base", - "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.12.0" - }, - "orig_nbformat": 4 - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/studies/Campbell/Resting-state/Campbell_Resting_EEG_Sustained_Attention_Healthy_Ageing_Cross_Sectional_LEISURE.ipynb b/studies/Campbell_Resting_EEG_Sustained_Attention_Healthy_Ageing_Cross_Sectional_LEISURE.ipynb similarity index 100% rename from studies/Campbell/Resting-state/Campbell_Resting_EEG_Sustained_Attention_Healthy_Ageing_Cross_Sectional_LEISURE.ipynb rename to studies/Campbell_Resting_EEG_Sustained_Attention_Healthy_Ageing_Cross_Sectional_LEISURE.ipynb diff --git "a/studies/Anij\303\244rv/Hermens_LABS_rsEEG_alpha_reactivity.ipynb" b/studies/Hermens_LABS_rsEEG_alpha_reactivity.ipynb similarity index 100% rename from "studies/Anij\303\244rv/Hermens_LABS_rsEEG_alpha_reactivity.ipynb" rename to studies/Hermens_LABS_rsEEG_alpha_reactivity.ipynb diff --git "a/studies/Anij\303\244rv/LABS_EEG_processing_pipeline.ipynb" b/studies/LABS_EEG_processing_pipeline.ipynb similarity index 100% rename from "studies/Anij\303\244rv/LABS_EEG_processing_pipeline.ipynb" rename to studies/LABS_EEG_processing_pipeline.ipynb diff --git "a/studies/Anij\303\244rv/Mills_LABS_rsEEG_classic_bp.ipynb" b/studies/Mills_LABS_rsEEG_classic_bp.ipynb similarity index 100% rename from "studies/Anij\303\244rv/Mills_LABS_rsEEG_classic_bp.ipynb" rename to studies/Mills_LABS_rsEEG_classic_bp.ipynb diff --git "a/studies/Anij\303\244rv/Mitchell_OKTOS_rsEEG_MSE+LZC.ipynb" b/studies/Mitchell_OKTOS_rsEEG_MSE+LZC.ipynb similarity index 100% rename from "studies/Anij\303\244rv/Mitchell_OKTOS_rsEEG_MSE+LZC.ipynb" rename to studies/Mitchell_OKTOS_rsEEG_MSE+LZC.ipynb diff --git "a/studies/Anij\303\244rv/Pace_LEISURE_rsEEG_aperiodic+classic_bp+alpha_reactivity.ipynb" b/studies/Pace_LEISURE_rsEEG_aperiodic+classic_bp+alpha_reactivity.ipynb similarity index 100% rename from "studies/Anij\303\244rv/Pace_LEISURE_rsEEG_aperiodic+classic_bp+alpha_reactivity.ipynb" rename to studies/Pace_LEISURE_rsEEG_aperiodic+classic_bp+alpha_reactivity.ipynb diff --git "a/studies/Anij\303\244rv/Pace_METACOG_erp_analysis.ipynb" b/studies/Pace_METACOG_erp_analysis.ipynb similarity index 100% rename from "studies/Anij\303\244rv/Pace_METACOG_erp_analysis.ipynb" rename to studies/Pace_METACOG_erp_analysis.ipynb diff --git "a/studies/Anij\303\244rv/mixedeffectsmodel_OKTOS_temp.ipynb" b/studies/mixedeffectsmodel_OKTOS_temp.ipynb similarity index 100% rename from "studies/Anij\303\244rv/mixedeffectsmodel_OKTOS_temp.ipynb" rename to studies/mixedeffectsmodel_OKTOS_temp.ipynb