Getting started with Cerca OPM data
Cerca Magnetics is a joint-venture spin-out company from the University of Nottingham and Magnetic Shields Limited that develops full wearable OPM-MEG systems. The data from their current systems is stored in the .cMEG format, but a converter is supplied with each system to turn the data into the .fif format for post-processing. The .fif format is already supported by FieldTrip, which means that no special functions are needed for reading the data. A separate function can be used to read the data straight from .cMEG into FieldTrip – this is currently being developed.
There are some differences which are relevant for processing Cerca data - and OPM data as a whole. One difference is how events or triggers are represented and detected. Another difference is in the procedure used to record the position of the OPM sensors relative to the head, which is relevant for topographic plotting and for the co-registration between MEG and MRI. Furthermore, Cerca OPM systems use either dual axis or triaxial systems, and so displaying data at a sensor-level requires some further considerations.
The example data provided on the our download server includes a self-contained version of the data in the .fif format organized according to the BIDS standard. This includes sensor locations that are already co-registered to a 3D digitization of the participant. An additional folder includes the raw files in the original .cMEG if you wish to develop a co-registration pipeline separately from what is provided. All files are detailed in the README.txt.
Read and visualize the continuous data
% this starts with the fif data organized according to BIDS
filename = 'ses-001/meg/sub-004_ses-001_task-braille_run-001_meg.fif';
cfg = [];
cfg.dataset = filename;
data = ft_preprocessing(cfg);
cfg = [];
cfg.preproc.demean = 'yes';
cfg.ylim = [-1 1] * 1e-11;
cfg.blocksize = 5; % seconds
ft_databrowser(cfg, data);
Reading and processing triggers
The Cerca OPM system uses triggers that are represented as voltages in an analogue channel.
In this case the experiment consisted of a braille pattern stimulation on the left and right index fingers. Each stimulus is coded with an analogue trigger represented in the “Bra” channel (for Braille), and the stimulated finger is represented by the “Lef” or “Rig” channels (Left or Right). Here we will look at both fingers. For ease of visualisation, we will also only look at the Z (radial) channels.
%% read the triggers and the data segments of interest
hdr = ft_read_header(filename);
disp(hdr.label)
% determine the segments of interest around the triggers (aka the trials)
cfg = [];
cfg.dataset = filename;
cfg.trialdef.eventtype = 'Bra';
cfg.trialdef.prestim = 0.1;
cfg.trialdef.poststim = 0.3;
cfg = ft_definetrial(cfg);
% read the segments of interest (aka the trials)
cfg.channel = '*Z'; % Only look at the Z channels
cfg.demean = 'yes';
cfg.baselinewindow = [-inf 0];
data = ft_preprocessing(cfg);
% remove large variance trials
cfg = [];
cfg.method = 'summary';
data_clean = ft_rejectvisual(cfg, data);
% average over trials and plot the ERF
cfg = [];
timelock = ft_timelockanalysis(cfg, data_clean);
plot(timelock.time, timelock.avg)
Reading and filtering
As there is quite some drift in the signal, it is better to apply a filter to the continuous data prior to segmenting. That results in the following alternate preprocessing script, similar to the tutorial on processing continuous EEG or MEG data. Here, we apply a highpass filter of 1 Hz and a lowpass filter of 40 Hz.
%% read and filter the continuous data
cfg = [];
cfg.dataset = filename;
cfg.channel = '*Z';
cfg.demean = 'yes';
cfg.baselinewindow = [-inf inf];
cfg.hpfilter = 'yes';
cfg.hpfreq = 2;
cfg.lpfilter = 'yes';
cfg.lpfreq = 40;
data_continuous = ft_preprocessing(cfg);
% determine the segments of interest around the triggers (aka the trials)
cfg = [];
cfg.dataset = filename;
cfg.trialdef.eventtype = 'Bra';
cfg.trialdef.prestim = 0.1;
cfg.trialdef.poststim = 0.3;
cfg = ft_definetrial(cfg);
% do the actual segmentation of the data
data_segmented = ft_redefinetrial(cfg, data_continuous);
% remove large variance trials
cfg = [];
cfg.method = 'summary';
data_clean = ft_rejectvisual(cfg, data_segmented);
% average over trials and plot the ERF
cfg = [];
timelock = ft_timelockanalysis(cfg, data_clean);
plot(timelock.time, timelock.avg)
Topographic plotting
When we look at this recording, we see channel names like “A1X” and “F4Y”. These names correspond to the hardware sensors, i.e. the OPMs, as they are named in the DAQ recording, with the sensitive axis indicated by the last letter (either X, Y, or Z).
The data structure returned by ft_preprocessing or the header from ft_read_header include the “grad” field, which describes the sensor positions (see also this FAQ). The recording software itself does not know the relative locations of sensors to each other in the helmet, but it requires you to select a _HelmConfig.tsv file that lists which slot each sensor is in, as well as their positions and orientations for each sensitive axis. In this specific recording the OPM sensors were placed in the generic purple adult large helmet supplied by Cerca.
For topographic plotting, the channel names in data.label should be mapped to the corresponding locations on the head. To keep the script clean and reproducible, we will not go in the data structure and manually rename the channels, but use a montage as explained here.
The _HelmConfig.tsv file specifies the channel names (based on the name of the hardware sensors) and their placements in the helmet. It also contains 2D projection coordinates of the helmet slot positions in Layx and Layy which can be used for 2D plotting. Using the _HelmConfig.tsv file, we can create a layout for the topographic plotting.
% Rename the channel labels using a montage
montage = [];
montage.labelold = data.label';
% Read the channel mapping
HelmConfig = tdfread(‘Braille_example_data\11766_054_braille1\20240129_101236_cMEG_Data\ 20240129_101236_HelmConfig.tsv');
last_line_idx = strmatch('Helmet:',HelmConfig.Name);
field_names = fieldnames(HelmConfig);
for n = 1:size(field_names,1)
HelmConfig.(field_names{n})(last_line_idx,:) = [];
if max(strcmpi(field_names{n},{'Px';'Py';'Pz';'Ox';'Oy';'Oz';'Layx';'Layz'}))
try
HelmConfig.(field_names{n}) = str2num(HelmConfig.(field_names{n}));
end
end
end
HC_sensnames = strrep(cellstr(HelmConfig.Sensor),' ','');
HC_slotnames = strrep(cellstr(HelmConfig.Name),' ','');
sensor_slot_mapping = [];
for n = size(data.label,1):-1:1
sensor_slot_mapping(n) = find(strcmpi(HC_sensnames,data.label{n}));
end
montage.labelnew = HC_slotnames(sensor_slot_mapping);
montage.tra = eye(size(data.label,1)); % the data is simply copied one-to-one
% Create layout from the HelmConfig.tsv file
layout = [];
layout.pos(:,1) = HelmConfig.Layx(sensor_slot_mapping);
layout.pos(:,2) = HelmConfig.Layy(sensor_slot_mapping)];
layout.label = HC_slotnames(sensor_slot_mapping);
layout.width = repmat(0.02,size(layout.pos,1),1);
layout.height = repmat(0.02,size(layout.pos,1),1);
% rename the channels in the raw data (this could have been done earlier in the processing pipeline)
cfg = [];
cfg.montage = montage;
data_helmetpos = ft_preprocessing(cfg, data);
% and rename the channels in the timelocked ERF data (not needed if the raw data channels had been renamed earlier on)
timelock_helmetpos = ft_preprocessing(cfg, timelock);
% plot the ERFs on the corresponding position in the helmet
cfg = [];
cfg.layout = layout;
cfg.showoutline = 'yes';
ft_multiplotER(cfg, timelock_helmetpos)
Co-registration
There is a dedicated tutorial elsewhere that deals with various ways to co-register the OPM sensors with the head and the anatomical MRI. In this demonstration, the .fif file already contains a digitization that the sensors have been co-registered to (using the dev2head transform hdr.orig.dev_head_t). However, the raw digitization files and CAD file for the helmet have also been provided on the download server (see README.txt). The co-registration of the sensors to the digitization was performed as described in Figure 2b in Hill et al. 2020.