Skip to content

Commit f02f209

Browse files
AlexanderNAAlexanderNA
committed
OPM Neuro1 tutorial
1 parent c4a9984 commit f02f209

File tree

8 files changed

+301
-0
lines changed

8 files changed

+301
-0
lines changed

‎docs/assets/figures/opm/aef_example_erf.PNG

81.3 KB
Loading

‎docs/assets/figures/opm/aef_example_psd_after.PNG

199 KB
Loading

‎docs/assets/figures/opm/aef_example_psd_after_filters.PNG

273 KB
Loading

‎docs/assets/figures/opm/aef_example_psd_before.PNG

398 KB
Loading

‎docs/assets/figures/opm/aef_example_topo.PNG

9.86 KB
Loading

‎docs/assets/figures/opm/aef_example_topo_ft.PNG

54.9 KB
Loading

‎docs/tutorials/opm/index.md

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,5 +18,6 @@ These tutorials cover the following areas in OPM analysis.
1818
3. [Coregistration](coreg/index.md)
1919
4. [Simulation](simulation/index.md)
2020
5. [Supported File Formats](formats/index.md)
21+
6. [Neuro1: Sensor Level AEF](neuro1_aef/index.md)
2122

2223
--8<-- "addons/abbreviations.md"

‎docs/tutorials/opm/neuro1_aef/index.md

Lines changed: 300 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,300 @@
1+
# Neuro-1 example: auditory evoked field
2+
3+
## Introduction
4+
5+
This code example will read in, process and analyse OPM-MEG data from a single subject recorded using a Neuro-1 QuSpin system. The aim here is to show how, with spm, a code-only approach can be used for reading in data, pre-processing, analysing and plotting data at sensor level.
6+
7+
## Download example data
8+
9+
The data used can be downloaded [here](https://www.fil.ion.ucl.ac.uk/spm/data/opm/opm-tutorial-neuro1.zip).
10+
11+
```matlab
12+
data_dir = 'D:\data\opm-tutorial-neuro1';
13+
cd(data_dir)
14+
```
15+
16+
## Read in the data
17+
18+
This data is in MEG BIDS format, so if there were more subjects the same script could be used by changing only these parameters.
19+
20+
```matlab
21+
BIDS = [];
22+
BIDS.sub = '001';
23+
BIDS.task = 'aef';
24+
BIDS.run = '001';
25+
BIDS.directory = data_dir;
26+
```
27+
Now use SPM to read in the Neuro-1 data and peripheral files.
28+
29+
Make sure you have downloaded spm. This tutorial has been tested on development version (spm/spm)
30+
[commit download link](https://github.com/spm/spm/archive/49e95f2d570528e2d8cccc95b621896a2456a3d5.zip)
31+
32+
```matlab
33+
% Run spm (change path as appropriate)
34+
addpath('D:\GitHub\spm')
35+
spm eeg
36+
37+
% Reference the meg and peripheral files
38+
S = [];
39+
S.category = 'meg';
40+
S.description = 'meg';
41+
S.type = '.lvm';
42+
S.derivative = 0;
43+
S.BIDS = BIDS;
44+
meg_dir = spm_BIDS_file(S);
45+
46+
S.description = 'positions';
47+
S.type = '.tsv';
48+
pos_dir = spm_BIDS_file(S);
49+
50+
S.description = 'channels';
51+
chan_dir = spm_BIDS_file(S);
52+
53+
S.category = 'anat';
54+
S.description = 'mri';
55+
S.type = '.nii';
56+
S.detailed = 0;
57+
anat_dir = spm_BIDS_file(S);
58+
59+
% Load the raw data
60+
S = [];
61+
S.data = meg_dir.file;
62+
S.sMRI = anat_dir.file;
63+
S.channels = chan_dir.file;
64+
S.positions = pos_dir.file;
65+
S.lead = 0; % We will calculate leadfields later
66+
S.fs = 375;
67+
D = spm_opm_create(S)
68+
```
69+
70+
Running spm_opm_create may take a few minutes while it processes the anatomical information. Files will be saved to disk which will speed up future loading.
71+
72+
## Pre-processing
73+
74+
We can now prepare the data for analysis by selecting bad channels, and applying temporal and spatial filters.
75+
One method for removing channels is to visually inspect the power spectral density (PSD). You can use spm_opm_psd to do this. By enabling S.selectbad, you can click on channels to inspect them and click and drag on the plot to select channels for removal. Outlier channel information is highlighted by a black star. Just close the figure when done and it outputs the selected (bad) indices.
76+
77+
```matlab
78+
% Identify and remove bad channels
79+
S = [];
80+
S.plot = 1;
81+
S.triallength = 3000;
82+
S.D = D;
83+
S.channels = 'MEG';
84+
S.selectbad = 1;
85+
S.plot = 1;
86+
[~,~,indices] = spm_opm_psd(S);
87+
88+
% Set bad channels
89+
sel_lab = chanlabels(S.D, S.channels);
90+
sel_indices = indchannel(D, sel_lab(indices));
91+
92+
D = badchannels(D, sel_indices, 1);
93+
D.save();
94+
95+
```
96+
97+
In this case we first see this figure:
98+
<figure markdown>
99+
<div class="center">
100+
<img src="../../../assets/figures/opm/aef_example_psd_before.PNG" style="width:160mm" />
101+
</div>
102+
</figure>
103+
104+
After selecting bad channels
105+
<figure markdown>
106+
<div class="center">
107+
<img src="../../../assets/figures/opm/aef_example_psd_after.PNG" style="width:160mm" />
108+
</div>
109+
</figure>
110+
111+
Working with only the remaining good channels, apply some filters.
112+
113+
```matlab
114+
% Bandpass filter (2-40 Hz)
115+
S = [];
116+
S.D = D;
117+
S.freq = [2 40];
118+
S.band = 'bandpass';
119+
S.prefix = 'filt_';
120+
filt_D = spm_eeg_ffilter(S);
121+
122+
% Bandstop filter (50 Hz)
123+
S = [];
124+
S.D = filt_D;
125+
S.freq = [48 52];
126+
S.band = 'stop';
127+
S.prefix = '';
128+
filt_D = spm_eeg_ffilter(S);
129+
130+
% Adaptive multipole modelling (AMM; https://doi.org/10.1002/hbm.26596)
131+
S = [];
132+
S.D = filt_D;
133+
S.prefix = '';
134+
filt_D = spm_opm_amm(S);
135+
```
136+
137+
Inspect the result
138+
139+
```matlab
140+
% Only consider good channels
141+
good_channels = filt_D.chanlabels(setdiff(indchantype(filt_D, 'MEG'),...
142+
badchannels(filt_D)));
143+
144+
% Plot the PSD, this time without selecting channels.
145+
S = [];
146+
S.plot = 1;
147+
S.triallength = 3000;
148+
S.D = filt_D;
149+
S.channels = good_channels;
150+
spm_opm_psd(S);
151+
```
152+
153+
<figure markdown>
154+
<div class="center">
155+
<img src="../../../assets/figures/opm/aef_example_psd_after_filters.PNG" style="width:160mm" />
156+
</div>
157+
</figure>
158+
159+
## Epoch the data
160+
161+
To avoid pre-processing the data each time you reload matlab or work with another dataset, it can be useful to load previously worked on data from disk. It is not necessary for this tutorial, but let's assume that is the case now.
162+
163+
```matlab
164+
% Restore the workspace
165+
try
166+
spm quit
167+
catch
168+
end
169+
restoredefaultpath
170+
clearvars
171+
close all
172+
addpath('D:\GitHub\spm')
173+
spm eeg
174+
175+
% BIDS
176+
data_dir = 'D:\data\opm-tutorial-neuro1';
177+
cd(data_dir)
178+
BIDS = [];
179+
BIDS.sub = '001';
180+
BIDS.task = 'aef';
181+
BIDS.run = '001';
182+
BIDS.directory = data_dir;
183+
184+
% Reference preprocessed data
185+
S = [];
186+
S.category = 'meg';
187+
S.description = 'meg';
188+
S.type = '.mat';
189+
S.derivative = 0;
190+
S.prefix = 'filt_';
191+
S.BIDS = BIDS;
192+
proc_dir = spm_BIDS_file(S);
193+
filt_D = spm_eeg_load(proc_dir.file);
194+
```
195+
196+
To calculate the event related field, we first need to identify events in the data. Contained within the data is a trigger channel 'A10' which recieved input from the system generating the auditory tones. Process that trigger.
197+
198+
```matlab
199+
S = [];
200+
S.D = filt_D;
201+
S.timewin = [-50 200];
202+
S.condLabels = {'tone'};
203+
S.bc = 1;
204+
S.triggerChannels = {'A10'};
205+
S.thresh = -0.16;
206+
[e_filt_D] = spm_opm_epoch_trigger(S);
207+
```
208+
209+
This should result in 603 trials.
210+
211+
## Sensor level analysis
212+
213+
```matlab
214+
% Average
215+
S = [];
216+
S.D = e_filt_D;
217+
S.prefix = 'avg_';
218+
avg_e_filt_D = spm_eeg_average(S);
219+
220+
% Plot ERF
221+
MEGind = indchantype(avg_e_filt_D,'MEGMAG');
222+
used = setdiff(MEGind,badchannels(avg_e_filt_D));
223+
pl = avg_e_filt_D(used,:,1)';
224+
figure();
225+
plot(avg_e_filt_D.time(),pl,"Color",[0.2 0.2 0.2 0.4])
226+
xlabel('Time (s)')
227+
ylabel('B (fT)')
228+
grid on
229+
ax = gca; % current axes
230+
ax.FontSize = 13;
231+
ax.TickLength = [0.02 0.02];
232+
fig= gcf;
233+
fig.Color=[1,1,1];
234+
xlim([-0.05 0.2])
235+
```
236+
237+
We can see the characteristic M100 peaking at ~110ms.
238+
239+
<figure markdown>
240+
<div class="center">
241+
<img src="../../../assets/figures/opm/aef_example_erf.PNG" style="width:160mm" />
242+
</div>
243+
</figure>
244+
245+
246+
```matlab
247+
% Make topoplot of radial channels around M100.
248+
S = [];
249+
S.D = avg_e_filt_D;
250+
S.mode = 'scalp';
251+
S.timewin = [90 125];
252+
S.channels = 'regexp_.*-Z.*';
253+
S.optimise = 1;
254+
S.prefix = 'M100_';
255+
spm_eeg_convert2images(S);
256+
257+
% Plot image
258+
image_dir = fullfile(BIDS.directory,['sub-', BIDS.sub], 'meg', ...
259+
['M100_avg_e_filt_','sub-', BIDS.sub, '_task-',BIDS.task, '_run-', BIDS.run,'_meg'] ...
260+
,['condition_', avg_e_filt_D.conditions{1}, '.nii']);
261+
spm_image('display',image_dir)
262+
```
263+
264+
The topography is sensible, showing two dipoles, left and right.
265+
266+
<figure markdown>
267+
<div class="center">
268+
<img src="../../../assets/figures/opm/aef_example_topo.PNG" style="width:160mm" />
269+
</div>
270+
</figure>
271+
272+
We can also visualise the image using Fieldtrip functions contained within SPM
273+
274+
```matlab
275+
% Convert to Fieldtrip structure and select good MEG channels only.
276+
ft_D = avg_e_filt_D.ftraw;
277+
cfg = [];
278+
cfg.channel = avg_e_filt_D.chanlabels(used);
279+
ft_D = ft_selectdata(cfg, ft_D);
280+
281+
% Use the layout produced from subject anatomy and sensor positions
282+
% https://doi.org/10.1111/ejn.70060)
283+
lay = avg_e_filt_D.lay;
284+
lay.mask{1} = lay.outline{1}; % Extrapolate to fill circle (demo purposes only!)
285+
286+
% Topoplot
287+
cfg = [];
288+
cfg.layout = lay;
289+
cfg.channel = {'*-Z'};
290+
cfg.xlim = [0.09 0.125];
291+
ft_topoplotER(cfg,ft_D);
292+
```
293+
294+
<figure markdown>
295+
<div class="center">
296+
<img src="../../../assets/figures/opm/aef_example_topo_ft.PNG" style="width:160mm" />
297+
</div>
298+
</figure>
299+
300+
--8<-- "addons/abbreviations.md"

0 commit comments

Comments
 (0)