Human Connectome Project: MEG

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

Feedback: Comments, bug reports, suggestions, questions
Email address (if you expect an answer):