diff --git a/workshop/cuttingeegx.md b/workshop/cuttingeegx.md index 80f957563..cbd317666 100644 --- a/workshop/cuttingeegx.md +++ b/workshop/cuttingeegx.md @@ -9,7 +9,207 @@ CuttingEEG is turning 10 years old! This milestone calls for a special edition, Alongside the [CuttingEEG X](https://cuttingeegx.org) conference we will be organizing some local workshops in the morning. -- Who: Julia Chauvet, Konstantinos Tsilimparis, Britta Westner, Tilman Stephani, Robert Oostenveld -- When: 30 to 31 October 2024 +- Who: Konstantinos Tsilimparis, Robert Oostenveld. Jan-Mathijs Schoffelen +- When: 31 October 2024 - Where: Nijmegen - See for more details + +## The tutorial + +We will run a tutorial on [preprocessing and coregistration of SQUID-based and OPM-based data](/workshop/cuttingeegx/squids_vs_opms). We will explore the differences in data analysis between the two MEG systems. + +## Getting started with the hands-on session + +In this workshop you will work on your own laptop computer. It will be an in-person event with no possibilities for hybrid or online attendance. + +### Wifi access + +If you don't have a eduroam account through your institution, it is possible to get a visitor access. Please let us know if you need a visitor access at the day of the workshop. + +### MATLAB + +For the hands-on sessions we assume that you have a computer with a relatively recent version of MATLAB installed (preferably < 5 years old, >= 2019a/b). + +### FieldTrip + +To get the most recent copy of FieldTrip, you can follow this [link](https://github.com/fieldtrip/fieldtrip/releases/tag/20240916), download the zip-file, and unzip it at a convenient location on your laptop's hard drive. + +Alternatively, you could do the following in the MATLAB command window (less work and faster downloading): + +``` +% create a folder that will contain the code and the data, and change directory +mkdir('cuttingeegx'); +cd('cuttingeegx'); + +% download and unzip fieldtrip into the newly created folder +url_fieldtrip = 'https://github.com/fieldtrip/fieldtrip/archive/refs/tags/20240916.zip'; +unzip(url_fieldtrip); +``` + +Upon completion of this step, the folder structure should look something like this: + +```bash +cuttingeegx/ +|-- fieldtrip-20240916 + |-- bin + |-- compat + |-- connectivity + |-- contrib + |-- external + |-- fileio + |-- forward + |-- inverse + |-- plotting + |-- preproc + |-- private + |-- qsub + |-- realtime + |-- specest + |-- src + |-- statfun + |-- template + |-- test + |-- trialfun + |-- utilities +``` +{% include markup/red %} +Downloading and unzipping can take up to ~10 minutes, so please download and unzip FieldTrip prior to the workshop. +{% include markup/end %} + +{% include markup/red %} +If you have downloaded and unzipped by hand, it could be that there's an 'extra folder layer' in your directory structure. We recommend that you remove this extra layer, i.e. move all content one level up. +{% include markup/end %} + +### Test your installation in advance + +To have a smooth experience - and to avoid having to spend precious debugging time during the hands-on sessions - we recommend that you [test your MATLAB and FieldTrip installation in advance](/workshop/cuttingeegx/test_installation). + +## The data used in this tutorial + +Next, we proceed with downloading the relevant data. The data that are used in the hands-on sessions, are stored on the FieldTrip [download-server](https://download.fieldtriptoolbox.org/workshop/cuttingeegx). You can either ‘click around’ using web browsers and/or explorer windows to grab the data that are needed, or instead (less work, at least if it works) execute the MATLAB code below. + +Please ensure that your present working directory is the ``'cuttingeegx'`` folder, which you created in the previous step. Open a new .m file in the text editor and run the following (Note: do not use the command line window): + +``` +% create a folder (within cuttingeegx) that will contain the data +mkdir('data'); +cd('data'); + +% Download the SQUID and OPM dataset +url = 'https://download.fieldtriptoolbox.org/workshop/cuttingeegx'; +recursive_download(url, pwd) + + +function recursive_download(webLocation, localFolder) + + % RECURSIVE_DOWNLOAD downloads a complete directory from a RESTful web service + % + % Use as + % recursive_download(webLocation, localFolder) + % + % See also WEBREAD, WEBSAVE, UNTAR, UNZIP, GUNZIP + + % Copyright (C) 2023, Konstantinos Tsilimparis + % + % This file is part of FieldTrip, see http://www.fieldtriptoolbox.org + % for the documentation and details. + % + % FieldTrip is free software: you can redistribute it and/or modify + % it under the terms of the GNU General Public License as published by + % the Free Software Foundation, either version 3 of the License, or + % (at your option) any later version. + % + % FieldTrip is distributed in the hope that it will be useful, + % but WITHOUT ANY WARRANTY; without even the implied warranty of + % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + % GNU General Public License for more details. + % + % You should have received a copy of the GNU General Public License + % along with FieldTrip. If not, see . + % + % $Id$ + + % Read the HTML content of the URL + htmlContent = webread(webLocation); + pattern = ''; + matches = regexp(htmlContent, pattern, 'tokens'); + + % Iterate over the matches + for i = 2:numel(matches) % Ignore i=1, which is the parent directory link: '../' + item = matches{i}{1}; + + if endsWith(item, '/') % It is a folder + % Create the necessary directories if they do not exist + subfolder = fullfile(localFolder, item); + if ~isfolder(subfolder) + mkdir(subfolder); + end + + % Recursively download the subfolder + subWebLocation = strcat(webLocation, '/', item); + recursive_download(subWebLocation, subfolder); + + else % It is a file + % Create the necessary directories if they do not exist + if ~isfolder(localFolder) + mkdir(localFolder); + end + + % Download the file + fileUrl = strcat(webLocation, '/', item); + localFilePath = fullfile(localFolder, item); + websave(localFilePath, fileUrl); + end + end +end +``` + +At this stage, you ideally have a directory structure that looks like the following one: + +```bash +cuttingeegx/ +|-- data + |-- sub-001 + |-- anat + |-- meg + |-- opm + |-- squid +|-- fieldtrip-20240916 + |-- bin + |-- compat + |-- connectivity + |-- contrib + |-- external + |-- fileio + |-- forward + |-- inverse + |-- plotting + |-- preproc + |-- private + |-- qsub + |-- realtime + |-- specest + |-- src + |-- statfun + |-- template + |-- test + |-- trialfun + |-- utilities +``` + + +Whenever starting a fresh MATLAB session, to configure the right FieldTrip paths, execute the following: + +``` +% change into the 'cuttingeegx' folder and then do the following +restoredefaultpath +addpath('fieldtrip-20240916'); +addpath(genpath('data')); +ft_defaults; +``` + + `ft_defaults` command ensures that all of FieldTrip's required subdirectories are added to the path. + +{% include markup/red %} +Furthermore, please do NOT add FieldTrip with all subdirectories, subdirectories will be added automatically when needed, and only when needed (see this [FAQ](/faq/installation)). +{% include markup/end %} diff --git a/workshop/cuttingeegx/squids_vs_opms.md b/workshop/cuttingeegx/squids_vs_opms.md index 2f703a3a6..b8d1a9810 100644 --- a/workshop/cuttingeegx/squids_vs_opms.md +++ b/workshop/cuttingeegx/squids_vs_opms.md @@ -5,4 +5,238 @@ tags: [cuttingeegx] # SQUIDs versus OPMs -To be implemented, see https://github.com/fieldtrip/website/issues/798 \ No newline at end of file +## Introduction + +Traditional MEG systems use superconducting quantum interference devices (SQUIDs), which require expensive cryogenic cooling with liquid helium and involve placing the sensors 1.5-2 cm from the scalp. Optically pumped magnetometers (OPMs), a new generation of MEG, eliminate the need for cryogenic cooling, allowing sensors to be placed closer to the scalp. OPMs are small individual sensors that allow the researcher to choose where to place them over the head. These sensors are held in place by 3D-printed helmets, which can be customized for each subject. + +Even though, both SQUID and OPM system detect the brain's magnetic fields, their data analysis steps differ a bit. In this tutorial we explore the differences in preprocessing and coregistration between the two systems. + +OPMs are flexible in their placement allowing for new recording strategies. These new recording strategies lead to new data analysis strategies. For example, in this tutorial we use a small number of OPMs (32 sensors) and do sequential recordings in which we position them over different places over the scalp. + +OPMs are magnetometers, as their name suggests. Magnetometers are more sensitive to environmental noise than gradiometers, which most SQUID systems have. Several data analysis algorithms to remove environmental noise have been proposed (see [Seymour et al (2022)](https://www.sciencedirect.com/science/article/pii/S1053811921011058?via%3Dihub) for more details). In this tutorial, we apply homogeneous field correction (HFC). HFC works better with a large number of sensors (more than 32). + +Unlike SQUID systems, which have standard coregistration procedures, OPMs don't have a single standard. In this tutorial, we coregister the OPMs with the MRI using an optical 3D scanner which captures the participant’s facial features along with the OPM helmet (Zetter et al., 2019). + + +This tutorial does not cover follow-up analyses (like source reconstruction) which in principle should not differ from the SQUID follow-up analyses, or alternative coregistration methods which are covered in the tutorial on [coregistration of Optically Pumped Magnetometer (OPM) data](tutorial/coregistration_opm/). + +## Background + + +In this tutorial we will use recordings made with 32 OPM sensors placed in an adult-sized “smart” helmet with a total of 144 slots. Each slot allows the sensor to move along a single axial direction. That way, the sensor can slide in its slot until it touches the head surface, regardless of the head size and shape. To limit head movements we mounted the helmet on a wooden plate. + +To acquire a measurement for each of the 144 helmet slots, we divided the experiment into six sequential recordings while maintaining the participant's head in a fixed position. We kept 9 sensors around the participant’s head fixed for each recordings. Since the OPM sensors touched the participant’s head and the helmet was mounted, these 9 sensors were able to keep the participant's head fixed throughout the experiment. In each recording we moved the remaining 23 sensors to fill every helmet slot. + +### The dataset used in this tutorial +The data for this tutorial was recorded with a 32-sensor FieldLine HEDscan v3 system with a so-called smart helmet. Each OPM sensor has one channel that measures the normal component of the magnetic field. + +We perform an experiment with left median nerve stimulation on a single participant using both the SQUID and the OPM system. We expect the activity to be modelled 20 ms post-stimulation with a dipole located at the right Primary Somatosensory (S1) area (Andersen & Dalal, 2021; Boto et al., 2017; Buchner et al., 1994)). + +The dataset can be downloaded from our download server. + +## Procedure +In this tutorial for the SQUIDs we will take the following steps: +- Define trials and read the data using **[ft_definetrial](/reference/ft_definetrial)** and **[ft_preprocessing](/reference/ft_preprocessing)** +- Removing artifacts using **[ft_rejectvisual](/reference/ft_rejectvisual)** +- Compute the averaged ERFs using **[ft_timelockanalysis](/reference/ft_timelockanalysis)** +- Visualize the results for all the channels with **[ft_multiplotER](/reference/ft_multiplotER)** +- Plot the 2D sensor topography for a specified latency with **[ft_topoplotER](/reference/ft_topoplotER)** +- Coregister MRI with SQUIDs using **[ft_volumerealign](/reference/ft_volumerealign)** +- Plot the 3D sensor topography for a specified latency with **[ft_plot_topo3d](/reference/plotting/ft_plot_topo3d)** + +For the OPMs we will take the following steps: +- Define trials and read the data using **[ft_definetrial](/reference/ft_definetrial)** and **[ft_preprocessing](/reference/ft_preprocessing)** +- Removing artifacts using **[ft_denoise_hfc](/reference/ft_denoise_hfc)** and **[ft_rejectvisual](/reference/ft_rejectvisual)** +- Compute the averaged ERFs using **[ft_timelockanalysis](/reference/ft_timelockanalysis)** +- Renaming duplicate channels +- Concatenate the data over the six recordings using **[ft_appendtimelock](/reference/ft_appendtimelock)** +- Add NaNs to the missing channels +- Visualize the results for all the channels with **[ft_multiplotER](/reference/ft_multiplotER)** +- Plot the 2D sensor topography for a specified latency with **[ft_topoplotER](/reference/ft_topoplotER)** +- Coregister MRI with OPMs using an optical 3D scanner. For this several functions are used: **[ft_volumerealign](/reference/ft_volumerealign)**, **[ft_read_headshape](/reference/fileio/ft_read_headshape)**, **[ft_meshrealign](/reference/ft_meshrealign)**, **[ft_defacemesh](/reference/ft_defacemesh)**, and **[ft_transform_geometry](/reference/utilities/ft_transform_geometry)**. +- Append sensors from the six recordings using **[ft_appendsens](/reference/ft_appendsens)** +- Plot the 3D sensor topography for a specified latency with **[ft_plot_topo3d](/reference/plotting/ft_plot_topo3d)** + + +## SQUID +### Preprocessing & trial definition + +We begin by loading the SQUID data and defining trials. In the experiment, the inter-trial interval ranged from 800-1200 ms. We select a 200 ms prestimulus and 400 ms poststimulus window. We then select trials where left median nerve stimulation occurred (trigger code = 1). + +``` +cfg = []; +cfg.dataset = 'M:\Documents\cuttingeegx-workshop\data\squid\subj001ses002_3031000.01_20231208_01.ds'; +cfg.trialdef.eventtype = 'UPPT001'; +cfg.trialdef.eventvalue = [1]; % select stimulation trials only (which are ~85 % of all the trials) +cfg.trialdef.prestim = 0.2; % ISI varied from 0.8 to 1.2 sec +cfg.trialdef.poststim = 0.4; +cfg = ft_definetrial(cfg); + +cfg.demean = 'yes'; +cfg.baselinewindow = [-Inf 0]; +cfg.continuous = 'yes'; % without that I get the error: 'requested data segment extends over a discontinuous trial boundary'. I guess I need cfg.continuous since 'subj001ses002_3031000.01_20231208_01.ds' is initially a continuous recording and I do not want to lose the ITI during preprocessing. +data_stim = ft_preprocessing(cfg); + +save data_stim data_stim +``` + +### Removing artifacts + +The summary method in **[ft_rejectvisual](/reference/ft_rejectvisual)** allows to manually remove trials and/or channels that have high variance. + +``` +load data_stim data_stim + +cfg = []; +cfg.method = 'summary'; +data_stim_clean = ft_rejectvisual(cfg,data_stim); + +save data_stim_clean data_stim_clean +``` + +### Compute ERFs + +We start by removing trials that have higher variance than 8e-25 Tesla-squared, then assess the channels. Removing trials also reduces channel variance, so no need to reject any channels. We then calculate the ERFs: + +``` +load data_stim_clean data_stim_clean + +cfg = []; +avg_stim = ft_timelockanalysis(cfg, data_stim_clean); + +save avg_stim avg_stim +``` +### Visualisation + +We plot the activity from all the sensors. + +``` +load avg_stim avg_stim + +% find the time window of the activity +cfg = []; +cfg.layout = 'CTF275_helmet'; +ft_multiplotER(cfg, avg_stim); % use interactive +``` + +We are now going to use the interactive feature of **[ft_multiplotER](/reference/ft_multiplotER)** to find our activity of interest. We select the sensors that are on top of the right primary somatosensory area since that is where we expect our activity to be localised. We see a negative peak around 20 ms (more specifically 35-50 ms) post-stimulation. We can select this time window to see the topography. The dipolar pattern the right primary somatosensory area is now visible. We can also plot this dipolar pattern with **[ft_topoplotER](/reference/ft_topoplotER)**. + +``` +% plot the activity at [0.035 0.050] +cfg = []; +cfg.xlim = [0.035 0.050]; +cfg.layout = 'CTF275_helmet'; +ft_topoplotER(cfg, avg_stim); + +print -dpng figures/cuttingeegx_topo_squid.png +``` + +### Coregister the anatomical MRI to the SQUID coordinate system + +We read the MRI, SQUID channels and the Polhemus headshape in memory. I recommend converting all quantities to SI units to ensure consistent units throughout your pipeline. + +``` +mri = ft_read_mri('M:\Documents\cuttingeegx-workshop\data\dicoms\00001_1.3.12.2.1107.5.2.19.45416.2022110716263882460203497.IMA'); +ctf275 = ft_read_sens('M:\Documents\cuttingeegx-workshop\data\squid\subj001ses002_3031000.01_20231208_01.ds', 'senstype', 'meg'); +shape = ft_read_headshape('M:\Documents\cuttingeegx-workshop\data\squid\subj001ses001ses002.pos'); + +mri = ft_convert_units(mri, 'm'); +ctf275 = ft_convert_units(ctf275, 'm'); +shape = ft_convert_units(shape, 'm'); + +save mri mri +save ctf275 ctf275 +save shape shape +``` + +We coregister the MRI with the SQUIDs by converting the MRI coordinate system to match that of the SQUIDs. In other words, we ensure that the three fiducials (nasion, LPA and RPA) defining the MRI coordinate system are aligned to the same three fiducials defining the SQUID coordinate system. + +``` +load shape shape +load ctf275 ctf275 + +% assign the 'ctf' coordinate system to the MRI +cfg = []; +cfg.method = 'interactive'; +cfg.coordsys = 'ctf'; +mri_realigned1 = ft_volumerealign(cfg, mri); + +save mri_realigned1 mri_realigned1 +``` + + +To refine the coregistration, we can also coregister the Polhemus head shape with the skin surface that is extracted from the previously aligned MRI. Now not only the fiducials but also other skin surface points acquired with the Polhemus are used for the alignment. + +``` +load mri_realigned1 mri_realigned1 + +cfg = []; +cfg.method = 'headshape'; +cfg.headshape.headshape = shape; +mri_realigned2 = ft_volumerealign(cfg, mri_realigned1); + +% check +cfg = []; +cfg.method = 'headshape'; +cfg.headshape.headshape = shape; +cfg.headshape.icp = 'no'; +ft_volumerealign(cfg, mri_realigned2); % takes too long + +save mri_realigned2 mri_realigned2 +``` + + +### 3D sensor topography +Let's plot the sensor topography in 3D using **[ft_plot_topo3d](/reference/plotting/ft_plot_topo3d)**. This will help us visualize where the sensors are positioned relative to the head. To make this clearer, we'll also display the scalp and brain surfaces. + +First, we'll prepare the surface mesh of the scalp and brain: + +``` +load mri_realigned2 mri_realigned2 + +cfg = []; +cfg.output = {'brain', 'scalp'}; +mri_segmented = ft_volumesegment(cfg, mri_realigned2); + +cfg = []; +cfg.tissue = {'brain'}; +mesh_brain = ft_prepare_mesh(cfg, mri_segmented); + +cfg = []; +cfg.tissue = {'scalp'}; +mesh_scalp = ft_prepare_mesh(cfg, mri_segmented); + +save mesh_brain mesh_brain +save mesh_scalp mesh_scalp +``` +Now, we'll plot the 3D sensor topography for the time window [0.035, 0.050] seconds, along with the brain and scalp: + +``` +sampling_rate = 1200; % in Hz +prestim = 0.2; + +I1 = (prestim+0.035)*sampling_rate; +I2 = (prestim+0.050)*sampling_rate; + +selected_avg = mean(avg_stim.avg(:, I1:I2), 2); + +figure; +ft_plot_mesh(mesh_scalp, 'facealpha', 0.5, 'facecolor', 'skin', 'edgecolor', 'none', 'edgecolor', 'skin' ) +hold on +ft_plot_mesh(mesh_brain, 'facecolor', 'brain', 'edgecolor', 'none'); +hold on +ft_plot_topo3d(pos275,selected_avg, 'facealpha', 0.9) +camlight +view([360 0]) + + +print -dpng M:\Documents\cuttingeegx-workshop\code\figures\squid\cuttingeegx_topo3d_squid.png +``` + +The SQUID sensors are located 1.5-2.0 cm from the scalp, resulting in a less focal sensor topography. + +## OPM +### Preprocessing +HFC is a method for denoising MEG data based on a spatially homogeneous model of the background magnetic field across the OPM array. This method has previously been used successfully for reducing magnetic interference in OPM magnetometers (see Supplementary Material Fig. A3; Hill et al., 2022; Mellor et al., 2023; Seymour et al., 2022) that are more sensitive to the environmental noise than the SQUID gradiometers. + +### Coregistration diff --git a/workshop/cuttingeegx/test_installation.md b/workshop/cuttingeegx/test_installation.md new file mode 100644 index 000000000..668d580bb --- /dev/null +++ b/workshop/cuttingeegx/test_installation.md @@ -0,0 +1,26 @@ +--- +title: Test your MATLAB and FieldTrip installation in advance +tags: [toolkit2024] +--- + +# Test your MATLAB and FieldTrip installation in advance + +Prior to the hands-on sessions, you need to check the functionality of the computational setup you have on your end. This is needed to hit the ground running, and to avoid to spend time on debugging the MATLAB environment during the hands on session. We recommend to use a clean install of a [recent version](https://github.com/fieldtrip/fieldtrip/releases/tag/20240916) of the toolbox. Once you have set this up, you may want to execute the below code, and check the output. + + fprintf('################################################################\n'); + fprintf('computer: %s\n', computer); + ver('MATLAB'); + + [ftver, ftpath] = ft_version; + fprintf('FieldTrip path is at: %s\n', ftpath); + fprintf('FieldTrip version is: %s\n', ftver); + fprintf('ttest is: %s\n', which('ttest')) + fprintf('imdilate is: %s\n', which('imdilate')); + fprintf('dpss is: %s\n', which('dpss')); + fprintf('fminunc is: %s\n', which('fminunc')); + fprintf('ft_read_data is: %s\n', which('ft_read_data')); + fprintf('runica is: %s\n', which('runica')); % should not be found yet, or the fieldtrip version + fprintf('spm is: %s\n', which('spm')); % should not be found yet, or the fieldtrip version + fprintf('################################################################\n'); + +If you get an error that says `Undefined function or variable 'ft_version'.` please add FieldTrip to the MATLAB path and try again. See [here](/faq/installation) for more information about adding FieldTrip to the MATLAB search path (Hint: You should use `'ft_defaults'` to add the necessary FieldTrip paths).