Brainstorm
  • Comments
  • Menu
    • Attachments
    • Versions
    • Raw Text
    • Print View
  • Login

Software

  • Introduction

  • Gallery

  • Download

  • Installation

Users

  • Tutorials

  • Forum

  • Courses

  • Community

  • Publications

Development

  • What's new

  • What's next

  • About us

  • Contact us

  • Contribute

Revision 32 as of 2017-01-31 17:38:19
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28

Human Connectome Project: Resting-state MEG

[WARNING: Tutorial under construction, not ready for public use]

Authors: Francois Tadel, Guiomar Niso, Elizabeth Bock, Sylvain Baillet

This tutorial explains 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. Here we will focus only on reproducing the results on resting MEG recordings presented in the OMEGA tutorial.

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.

Contents

  1. License
  2. Presentation of the experiment
  3. Download and installation
  4. Import the anatomy
  5. Access the recordings
  6. Pre-processing
  7. Auxiliary channels
  8. Artifact cleaning
  9. Source estimation
  10. Power maps
  11. References
  12. Scripting

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.

Presentation of the experiment

Experiment

  • 1 subject x 6 minute resting session
  • Supine position
  • Eyes open with a red fixation cross in a darkened room

MEG acquisition

  • Acquisition at 2035Hz, with whole head MAGNES 3600 (4D Neuroimaging) system housed in a magnetically shielded room with 248 magnetometer channels

  • Recorded at the Saint Louis University (SLU) medical campus
  • Anti-aliasing low-pass filter at 400Hz
  • Recorded channels (at least 287):
    • 23 MEG reference sensors
    • 248 MEG axial gradiometers
    • 1 ECG bipolar (#278, #283)
    • 1 vertical EOG bipolar (#279, #284)
    • 1 horizontal EOG bipolar (#280, #285)
    • 4 EMG: left foot (#281), left hand (#282), right foot (#286), right hand (#287)
  • Noise recordings (empty room) are acquired just before the subject recordings.

Head shape and fiducial points

  • 3D digitization using a Polhemus FASTRAK-III (head localization coils and head shape).
  • But the 4D file "hs_file" was removed from the distributed package (see not below).

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.

  • In section "WU-Minn HCP Data - 900 Subjects + 7T", select "Explore subjects > MEG subjects"

    connectomedb1.gif

  • In tab "Subject information", navigate until you find subject 175237.

    connectomedb2.gif

  • Click on it, then click on button "Download images".

    connectomedb3.gif

  • Navigate in the available recordings and click on "Queue for download" for the following packages:
    • 3T MRI > Unprocessed > Structural unprocessed

    • MEG > Unprocessed > Noise

    • MEG > Unprocessed > Resting state

    • MEG > Processed > Anatomy

  • Click on "Download packages" > "Download now" > Follow the instructions.

    connectomedb4.gif

  • Unzip all the downloaded files in the same folder. Note that this folder should not be 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".

Import the anatomy

Create a new subject and import the anatomy, partly processed with the HCP megconnectome v3 scripts.

  • Select the menu File > New subject > Subject name: 175237

    newsubj.gif

  • In the anatomy view, right-click on the subject folder > Import anatomy folder:
    File format: HCP MEG/anatomy (pipeline v3)
    File name: HCP/175237/MEG/anatomy

    import_anat.gif

  • The head surface (generated by Brainstorm) and the cortex surface (imported from the anatomy folder) are shown at the end of the process. You can also display the MRI of the subject.

    anat_surfaces.gif anat_mri.gif

  • The position of the anatomical landmarks (NAS, LPA, RPA) may not match exactly the subject's head, they are just set to standard positions in BTi coordinates (equivalent to SCS coordinates in Brainstorm). Same for the AC/PC points, which are defined in MNI coordinates (the linear MNI transformation is computed with the SPM function spm_maff8).

  • All the information allowing the identification of the subject has been removed from the dataset: the digitized head shape and the positions of the anatomical landmarks are not distributed, and the facial features have been blurred in the MRI. It is therefore impossible to coregister accurately the MRI and the position of the MEG sensors.
  • The FieldTrip team computed this registration with their tools and distribute only the final transformation matrix (MEG/anatomy/*_MEG_anatomy_transform.txt). This matrix is imported automatically here and gives good results, but unfortunately we cannot reproduce or double-check the quality of the registration.

Access the recordings

Link the resting MEG recordings to the Brainstorm database and run some basic quality control.

  • Switch to the functional view of the protocol (second button above the database explorer).
  • Right-click on the subject > Review raw file:
    File format: MEG/EEG: 4D-Neuroimaging/BTi
    File name: /175237/unprocessed/MEG/3-Restin/4D/c,rfDC

    meg_link.gif

  • File name: /175237/unprocessed/MEG/3-Restin/4D/c,rfDC
  • There are no events to read for these resting recordings, do not select any technical track:

    meg_noevt.gif

  • The registration MEG-MRI looks good (grey=head from the MRI, yellow=inside of the MEG helmet):

    meg_coreg.gif

  • Repeat the same operation with file containing the MEG room noise:
    Review raw file: /175237/unprocessed/MEG/1-Rnoise/4D/c,rfDC

    meg_noise.gif

  • The registration looks wrong, but this is normal: there is no subject in the MEG.
  • In Process1, drag and drop the two links to resting and noise recordings.
  • Run process: Frequency > Power spectrum density (Welch): All file, 4s, 50% overlap, Individual.

    psd_process.gif

  • Check the quality of the PSD for all the files, as documented in this tutorial. The sensor A244 looks extremely noisy compared with all the others, and others are suspicious (A227, A246, A248). We will mark them all as bad.

    psd_display.gif

  • Double-click on the link to open the rest recordings, select the montage "4D 218-248".
  • In the Record tab, select the button [Display mode for time series].
  • Select the bad channels (227, 244, 246, 248), right-click > Channels > Mark selected as bad.

    bad_a244.gif

Pre-processing

Apply frequency filters to the recordings.

  • In Process1, keep the same files selected and click on [Run].
  • Select process: Pre-process > Notch filter: 60 120 180 240 300 Hz, Process the entire file at once

  • Add process: Pre-process > Band-pass filter: High-pass filter at 0.3Hz, 60dB, Process entire file.

  • Add process: Frequency > Power spectrum density (Welch): Same options

    filter_process.gif

  • Run, then check the PSD after filtering.
  • Delete the folders corresponding to the original recordings (1-Rnoise, 3-Restin) and the notch filtered data (_notch). Only keep the fully processed files (_notch_high).

    psd_after.gif

Advanced

Auxiliary channels

The ECG, VEOG and HEOG are available in the recordings in the form of bipolar montages. For each of them, we have access to two channels that need to be subtracted before being visualized or used for artifact detection. If you display the EOG channels, you will observe something like this:

  • eog_bipolar.gif

To display these channels properly you can create the following montage, with the help of the Montage editor tutorial.

  • eog_montage.gif

Artifact cleaning

We will now run the automatic procedure for cleaning the heartbeats, as described in the introduction tutorials (detection, SSP). The results we want to illustrate here are robust enough, the recordings do not need to be processed any further. If you want to improve the quality of the data with more manual cleaning (blinks, saccades, movements, bad segments), please refer to the introduction tutorials.

  • Double-click on the resting recordings to open the MEG signals.
  • In the Record tab, select menu: Artifacts > Detect heartbeats:
    Channel name: "ECG+, -ECG-" (bipolar montage), All file, cardiac

    ecg_detect.gif

  • In the Record tab, select menu: Artifacts > SSP: Heartbeats: cardiac, MEG, Use existing SSP
    Run, evaluate the components, select the first one after making sure it removes the cardiac peaks.

    ssp_ecg.gif

Source estimation

  • In Process1, select the noise recordings (1-Rnoise).

  • Run process: Sources > Compute covariance: All file, Noise covariance, Copy to other folders.
    noiscov_process.gif

  • In Process1, select the subject rest recordings (3-Restin).

  • Run process: Sources > Compute head model: Cortex surface, MEG=Overlapping spheres.

    headmodel_process.gif

  • Run process: Sources > Compute sources [2016]: Kernel only, one per file, dSPM, constrained.

    inverse_process.gif

Power maps

We have now access to continuous recordings in source space, we can estimate the average power over this resting period, for various frequency bands. We propose here to compute the PSD from the first 100s instead of the full 379s rest recordings, it is faster and leads to very similar results. For more stable estimates, you can compute the the PSD over all the data available instead.

  • In Process1, select the sources for the subject rest recordings (3-Restin)

  • Select process: Frequency > Power spectrum density (Welch): [0,100s], Window=4s, 50% overlap, Group in frequency bands (use the default frequency bands), Save individual PSD values.

  • Add process: Standardize > Spectrum normalization: Relative power (divide by total power). This normalization is done independently for each source, the output for each frequency band is a percentage of how much this frequency band contributes to the total power of the source signal.

  • Add process: Sources > Spatial smoothing: FWHM=3mm, Overwrite.

    psd_sources_process.gif

    psd_sources_normalize.gif

  • Double-click on the average to display it. To create the figures below, right-click on the figure > Snapshot > Frequency contact sheet.

    psd_sources_left.gif psd_sources_top.gif

  • These figures are very similar to what was presented in the in the OMEGA tutorial. Please refer to it for the interpretation of the results.

  • You can compute the same results in the head volume. You just have to use a volume head model on a group grid instead of the one we used here. This is explained in the tutorial Group analysis. This seems to confirm that the power in the delta and gamma bands are mostly coming with non-brain regions (eyes and neck muscles).

    SCREEN CAPTURES

  • You can produce similar plots at the sensor level: Compute the PSD by frequency band of the rest recordings for each subject, then compute the average across subjects (not as reliable as averaging in source space, as explained here).

    SCREEN CAPTURES

References

  • ConnectomeDB

  • 500 Subjects + MEG2 Data Release Reference

  • HCP-MEG NeuroImage article

  • MEG Pipeline Reference (megconnectome v3)

  • MEG1 release reference manual

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





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


  • MoinMoin Powered
  • Python Powered
  • GPL licensed
  • Valid HTML 4.01