Human Connectome Project: MEG
[WARNING: Tutorial under construction, not ready for public use]
Authors: Francois Tadel
This tutorial introduces how to download MEG recordings from the Human Connectome Project (HCP) ConnectomeDB database and process them into Brainstorm. The original processing pipeline was described in this article and this reference manual.
Note that the operations used here are not detailed, the goal of this tutorial is not to introduce Brainstorm to new users. For in-depth explanations of the interface and theoretical foundations, please refer to the introduction tutorials.
License
These data were generated and made available by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657), which is funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research and by the McDonnell Center for Systems Neuroscience at Washington University.
For additional information on how to acknowledge HCP and cite HCP publications if you have used data provided by the WU-Minn HCP consortium, see http://www.humanconnectome.org/citations.
As a reminder, users of these datasets must comply with the Data Use Terms that were agreed upon when receiving these data.
Download and installation
We will use only one subject available in the HCP-MEG2 distribution: subject #175237.
First, make sure you have at least 20Gb of free space on your hard drive.
Login in or create an account on the ConnectomeDB website.
- Unzip it in a folder that is not in any of the Brainstorm folders (program folder or database folder).
Start Brainstorm (Matlab scripts or stand-alone version). For help, see the Installation page.
Select the menu File > Create new protocol. Name it "TutorialHcp" and select the options:
"No, use individual anatomy",
"No, use one channel file per condition".
megconnectome pipeline v3.0
http://www.humanconnectome.org/documentation/HCP-pipelines/meg-pipeline.html
Scripting
The following script from the Brainstorm distribution reproduces the analysis presented in this tutorial page: brainstorm3/toolbox/script/tutorial_hcp.m
1 function tutorial_hcp(tutorial_dir)
2 % TUTORIAL_HCP: Script that reproduces the results of the online tutorial "Human Connectome Project: Resting-state MEG".
3 %
4 % CORRESPONDING ONLINE TUTORIALS:
5 % https://neuroimage.usc.edu/brainstorm/Tutorials/HCP-MEG
6 %
7 % INPUTS:
8 % tutorial_dir: Directory where the HCP files have been unzipped
9
10 % @=============================================================================
11 % This function is part of the Brainstorm software:
12 % https://neuroimage.usc.edu/brainstorm
13 %
14 % Copyright (c) University of Southern California & McGill University
15 % This software is distributed under the terms of the GNU General Public License
16 % as published by the Free Software Foundation. Further details on the GPLv3
17 % license can be found at http://www.gnu.org/copyleft/gpl.html.
18 %
19 % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
20 % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
21 % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
22 % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
23 % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
24 %
25 % For more information type "brainstorm license" at command prompt.
26 % =============================================================================@
27 %
28 % Author: Francois Tadel, 2017
29
30
31 %% ===== FILES TO IMPORT =====
32 % You have to specify the folder in which the tutorial dataset is unzipped
33 if (nargin == 0) || isempty(tutorial_dir) || ~file_exist(tutorial_dir)
34 error('The first argument must be the full path to the tutorial dataset folder.');
35 end
36 % Subject name
37 SubjectName = '175237';
38 % Build the path of the files to import
39 AnatDir = fullfile(tutorial_dir, SubjectName, 'MEG', 'anatomy');
40 Run1File = fullfile(tutorial_dir, SubjectName, 'unprocessed', 'MEG', '3-Restin', '4D', 'c,rfDC');
41 NoiseFile = fullfile(tutorial_dir, SubjectName, 'unprocessed', 'MEG', '1-Rnoise', '4D', 'c,rfDC');
42 % Check if the folder contains the required files
43 if ~file_exist(AnatDir) || ~file_exist(Run1File) || ~file_exist(NoiseFile)
44 error(['The folder ' tutorial_dir ' does not contain subject #175237 from the HCP-MEG distribution.']);
45 end
46
47
48 %% ===== CREATE PROTOCOL =====
49 % The protocol name has to be a valid folder name (no spaces, no weird characters...)
50 ProtocolName = 'TutorialHcp';
51 % Start brainstorm without the GUI
52 if ~brainstorm('status')
53 brainstorm nogui
54 end
55 % Delete existing protocol
56 gui_brainstorm('DeleteProtocol', ProtocolName);
57 % Create new protocol
58 gui_brainstorm('CreateProtocol', ProtocolName, 0, 0);
59 % Start a new report
60 bst_report('Start');
61
62
63 %% ===== IMPORT DATA =====
64 % Process: Import anatomy folder
65 bst_process('CallProcess', 'process_import_anatomy', [], [], ...
66 'subjectname', SubjectName, ...
67 'mrifile', {AnatDir, 'HCPv3'}, ...
68 'nvertices', 15000);
69
70 % Process: Create link to raw files
71 sFilesRun1 = bst_process('CallProcess', 'process_import_data_raw', [], [], ...
72 'subjectname', SubjectName, ...
73 'datafile', {Run1File, '4D'}, ...
74 'channelalign', 1);
75 sFilesNoise = bst_process('CallProcess', 'process_import_data_raw', [], [], ...
76 'subjectname', SubjectName, ...
77 'datafile', {NoiseFile, '4D'}, ...
78 'channelalign', 1);
79 sFilesRaw = [sFilesRun1, sFilesNoise];
80
81
82 %% ===== PRE-PROCESSING =====
83 % Process: Notch filter: 60Hz 120Hz 180Hz 240Hz 300Hz
84 sFilesNotch = bst_process('CallProcess', 'process_notch', sFilesRaw, [], ...
85 'freqlist', [60, 120, 180, 240, 300], ...
86 'sensortypes', 'MEG, EEG', ...
87 'read_all', 1);
88
89 % Process: High-pass:0.3Hz
90 sFilesBand = bst_process('CallProcess', 'process_bandpass', sFilesNotch, [], ...
91 'sensortypes', 'MEG, EEG', ...
92 'highpass', 0.3, ...
93 'lowpass', 0, ...
94 'attenuation', 'strict', ... % 60dB
95 'mirror', 0, ...
96 'useold', 0, ...
97 'read_all', 1);
98
99 % Process: Power spectrum density (Welch)
100 sFilesPsdAfter = bst_process('CallProcess', 'process_psd', sFilesBand, [], ...
101 'timewindow', [0 100], ...
102 'win_length', 4, ...
103 'win_overlap', 50, ...
104 'sensortypes', 'MEG, EEG', ...
105 'edit', struct(...
106 'Comment', 'Power', ...
107 'TimeBands', [], ...
108 'Freqs', [], ...
109 'ClusterFuncTime', 'none', ...
110 'Measure', 'power', ...
111 'Output', 'all', ...
112 'SaveKernel', 0));
113
114 % Mark bad channels
115 bst_process('CallProcess', 'process_channel_setbad', sFilesBand, [], ...
116 'sensortypes', 'A227, A244, A246, A248');
117
118 % Process: Snapshot: Frequency spectrum
119 bst_process('CallProcess', 'process_snapshot', sFilesPsdAfter, [], ...
120 'target', 10, ... % Frequency spectrum
121 'modality', 1); % MEG (All)
122
123 % Process: Delete folders
124 bst_process('CallProcess', 'process_delete', [sFilesRaw, sFilesNotch], [], ...
125 'target', 2); % Delete folders
126
127
128 %% ===== ARTIFACT CLEANING =====
129 % Process: Select data files in: */*
130 sFilesBand = bst_process('CallProcess', 'process_select_files_data', [], [], ...
131 'subjectname', 'All');
132
133 % Process: Select file names with tag: 3-Restin
134 sFilesRest = bst_process('CallProcess', 'process_select_tag', sFilesBand, [], ...
135 'tag', '3-Restin', ...
136 'search', 1, ... % Search the file names
137 'select', 1); % Select only the files with the tag
138
139 % Process: Detect heartbeats
140 bst_process('CallProcess', 'process_evt_detect_ecg', sFilesRest, [], ...
141 'channelname', 'ECG+, -ECG-', ...
142 'timewindow', [], ...
143 'eventname', 'cardiac');
144
145 % Process: SSP ECG: cardiac
146 bst_process('CallProcess', 'process_ssp_ecg', sFilesRest, [], ...
147 'eventname', 'cardiac', ...
148 'sensortypes', 'MEG', ...
149 'usessp', 1, ...
150 'select', 1);
151
152 % Process: Snapshot: Sensors/MRI registration
153 bst_process('CallProcess', 'process_snapshot', sFilesRest, [], ...
154 'target', 1, ... % Sensors/MRI registration
155 'modality', 1, ... % MEG (All)
156 'orient', 1); % left
157
158 % Process: Snapshot: SSP projectors
159 bst_process('CallProcess', 'process_snapshot', sFilesRest, [], ...
160 'target', 2, ... % SSP projectors
161 'modality', 1); % MEG (All)
162
163
164 %% ===== SOURCE ESTIMATION =====
165 % Process: Select file names with tag: task-rest
166 sFilesNoise = bst_process('CallProcess', 'process_select_tag', sFilesBand, [], ...
167 'tag', '1-Rnoise', ...
168 'search', 1, ... % Search the file names
169 'select', 1); % Select only the files with the tag
170
171 % Process: Compute covariance (noise or data)
172 bst_process('CallProcess', 'process_noisecov', sFilesNoise, [], ...
173 'baseline', [], ...
174 'sensortypes', 'MEG', ...
175 'target', 1, ... % Noise covariance (covariance over baseline time window)
176 'dcoffset', 1, ... % Block by block, to avoid effects of slow shifts in data
177 'identity', 0, ...
178 'copycond', 1, ...
179 'copysubj', 0, ...
180 'replacefile', 1); % Replace
181
182 % Process: Compute head model
183 bst_process('CallProcess', 'process_headmodel', sFilesRest, [], ...
184 'sourcespace', 1, ... % Cortex surface
185 'meg', 3); % Overlapping spheres
186
187 % Process: Compute sources [2018]
188 sSrcRest = bst_process('CallProcess', 'process_inverse_2018', sFilesRest, [], ...
189 'output', 2, ... % Kernel only: one per file
190 'inverse', struct(...
191 'Comment', 'dSPM: MEG', ...
192 'InverseMethod', 'minnorm', ...
193 'InverseMeasure', 'dspm2018', ...
194 'SourceOrient', {{'fixed'}}, ...
195 'Loose', 0.2, ...
196 'UseDepth', 1, ...
197 'WeightExp', 0.5, ...
198 'WeightLimit', 10, ...
199 'NoiseMethod', 'reg', ...
200 'NoiseReg', 0.1, ...
201 'SnrMethod', 'fixed', ...
202 'SnrRms', 1e-06, ...
203 'SnrFixed', 3, ...
204 'ComputeKernel', 1, ...
205 'DataTypes', {{'MEG'}}));
206
207
208 %% ===== POWER MAPS =====
209 % Process: Power spectrum density (Welch)
210 sSrcPsd = bst_process('CallProcess', 'process_psd', sSrcRest, [], ...
211 'timewindow', [0, 100], ...
212 'win_length', 4, ...
213 'win_overlap', 50, ...
214 'clusters', {}, ...
215 'scoutfunc', 1, ... % Mean
216 'edit', struct(...
217 'Comment', 'Power,FreqBands', ...
218 'TimeBands', [], ...
219 'Freqs', {{'delta', '2, 4', 'mean'; 'theta', '5, 7', 'mean'; 'alpha', '8, 12', 'mean'; 'beta', '15, 29', 'mean'; 'gamma1', '30, 59', 'mean'; 'gamma2', '60, 90', 'mean'}}, ...
220 'ClusterFuncTime', 'none', ...
221 'Measure', 'power', ...
222 'Output', 'all', ...
223 'SaveKernel', 0));
224
225 % Process: Spectrum normalization
226 sSrcPsdNorm = bst_process('CallProcess', 'process_tf_norm', sSrcPsd, [], ...
227 'normalize', 'relative', ... % Relative power (divide by total power)
228 'overwrite', 0);
229
230 % Process: Spatial smoothing (3.00)
231 sSrcPsdNorm = bst_process('CallProcess', 'process_ssmooth_surfstat', sSrcPsdNorm, [], ...
232 'fwhm', 3, ...
233 'overwrite', 1);
234
235 % Screen capture of final result
236 hFig = view_surface_data([], sSrcPsdNorm.FileName);
237 set(hFig, 'Position', [200 200 200 200]);
238 hFigContact = view_contactsheet(hFig, 'freq', 'fig');
239 bst_report('Snapshot', hFigContact, sSrcPsdNorm.FileName, 'Power');
240 close([hFig, hFigContact]);
241
242 % Save and display report
243 ReportFile = bst_report('Save', []);
244 bst_report('Open', ReportFile);
245
246
247
248