|
Size: 34466
Comment:
|
Size: 31322
Comment:
|
| Deletions are marked like this. | Additions are marked like this. |
| Line 1: | Line 1: |
| = MEG visual tutorial: Single subject = ''Authors: Francois Tadel, Elizabeth Bock.'' The aim of this tutorial is to reproduce the analysis and results o The aim of this tutorial is to provide high-quality recordings of a simple auditory stimulation and illustrate the best analysis paths possible with Brainstorm and FieldTrip. This page presents the workflow in the Brainstorm environment, the equivalent documentation for the FieldTrip environment will be available on the [[http://fieldtrip.fcdonders.nl/|FieldTrip website]]. Note that the operations used here are not detailed, the goal of this tutorial is not to teach Brainstorm to a new inexperienced user. For in depth explanations of the interface and the theory, please refer to the introduction tutorials. |
= MEG visual tutorial: Single subject (BIDS) = '''[WORK IN PROGRESS: THIS TUTORIAL IS NOT READY FOR PUBLIC USE]''' ''Authors: Francois Tadel, Elizabeth Bock. '' The aim of this tutorial is to reproduce in the Brainstorm environment the analysis described in the SPM tutorial "[[ftp://ftp.mrc-cbu.cam.ac.uk/personal/rik.henson/wakemandg_hensonrn/Publications/SPM12_manual_chapter.pdf|Multimodal, Multisubject data fusion]]". The data processed here consists of simultaneous MEG/EEG recordings from 16 participants performing a simple visual recognition task from presentations of famous, unfamiliar and scrambled faces. The analysis is split in two tutorial pages: the present tutorial describes the detailed analysis of one single subject; the second tutorial describes batch processing and [[Tutorials/VisualGroup|group analysis of all 16 participants]]. 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 [[http://neuroimage.usc.edu/brainstorm/Tutorials#Get_started|introduction tutorials]]. |
| Line 12: | Line 14: |
| == Download and installation == * '''Requirements''': You have already followed all the basic tutorials and you have a working copy of Brainstorm installed on your computer. * Go to the [[http://neuroimage.usc.edu/brainstorm3_register/download.php|Download]] page of this website, and download the file: '''sample_auditory.zip ''' * Unzip it in a folder that is not in any of the Brainstorm folders (program folder or database folder). This is really important that you always keep your original data files in a separate folder: the program folder can be deleted when updating the software, and the contents of the database folder is supposed to be manipulated only by the program itself. * Start Brainstorm (Matlab scripts or stand-alone version) * Select the menu File > Create new protocol. Name it "'''TutorialAuditory'''" and select the options: |
== License == This dataset was obtained from the OpenfMRI project (http://www.openfmri.org), accession #ds117. It is made available under the Creative Commons Attribution 4.0 International Public License. Please cite the following reference if you use these data: <<BR>>Wakeman DG, Henson RN, [[http://www.nature.com/articles/sdata20151|A multi-subject, multi-modal human neuroimaging dataset]], Scientific Data (2015) Any questions regarding the data, please contact: rik.henson@mrc-cbu.cam.ac.uk == Presentation of the experiment == ==== Experiment ==== * 16 subjects (the original version of the dataset included 19 subjects, but 3 were excluded from the group analysis for [[http://neuroimage.usc.edu/brainstorm/Tutorials/VisualSingleOrig#Bad_subjects|various reasons]]) * 6 acquisition runs of approximately 10mins for each subject * Presentation of series of images: familiar faces, unfamiliar faces, phase-scrambled faces * Participants had to judge the left-right symmetry of each stimulus * Total of nearly 300 trials for each of the 3 conditions ==== MEG acquisition ==== * Acquisition at 1100Hz with an Elekta-Neuromag VectorView system (simultaneous MEG+EEG). * Recorded channels (404): * 102 magnetometers * 204 planar gradiometers * 70 EEG electrodes recorded with a nose reference. * MEG data have been "cleaned" using Signal-Space Separation as implemented in MaxFilter 2.2. * A Polhemus device was used to digitize three fiducial points and a large number of other points across the scalp, which can be used to coregister the M/EEG data with the structural MRI image. * Stimulation triggers: The triggers related with the visual presentation are saved in the STI101 channel, with the following event codes (bit 3 = face, bit 4 = unfamiliar, bit 5 = scrambled): * Famous faces: 5 (00101), 6 (00110), 7 (00111) * Unfamiliar faces: 13 (01101), 14 (01110), 15 (01111) * Scrambled images: 17 (10001), 18 (10010), 19 (10011) * Delays between the trigger in STI101 and the actual presentation of stimulus: '''34.5ms''' * The data distribution includes MEG noise recordings acquired around the dates of the experiment, processed with MaxFilter 2.2 in the same way as the experimental data. ==== Subject anatomy ==== * MRI data acquired on a 3T Siemens TIM Trio: 1x1x1mm T1-weighted structural MRI. * The face was removed from the strucural images for anonymization purposes. * Processed with FreeSurfer 5.3. == Download and installation [TODO] == * First, make sure you have enough space on your hard drive, at least '''350Gb''': * Raw files: '''100Gb''' * Processed files: '''250Gb''' * The data is hosted on the OpenfMRI website: https://openfmri.org/dataset/ds000117/ * Download the following files for revision 1.0.0 (approximately 70Gb): * [[https://s3.amazonaws.com/openneuro/ds000117/ds000117_R1.0.0/compressed/ds000117_R1.0.0_derivatives_sub01-04.zip|Derivatives for subjects 01-04 (16.6 GB)]] * [[https://s3.amazonaws.com/openneuro/ds000117/ds000117_R1.0.0/compressed/ds000117_R1.0.0_derivatives_sub05-08.zip|Derivatives for subjects 05-08 (16.1 GB)]] * [[https://s3.amazonaws.com/openneuro/ds000117/ds000117_R1.0.0/compressed/ds000117_R1.0.0_derivatives_sub09-12.zip|Derivatives for subjects 09-12 (16.6 GB)]] * [[https://s3.amazonaws.com/openneuro/ds000117/ds000117_R1.0.0/compressed/ds000117_R1.0.0_derivatives_sub13-16.zip|Derivatives for subjects 13-16 (16.5 GB)]] * [[https://s3.amazonaws.com/openneuro/ds000117/ds000117_R1.0.0/compressed/ds000117_R1.0.0_metadata.zip|Metadata (235.4 KB)]] * FreeSurfer segmentation: '''NOT AVAILABLE YET [TODO]''' * Empty room recordings: '''NOT AVAILABLE YET ''''''[TODO]''' * Unzip all the .zip files in the same folder.<<BR>>'''Reminder''': Do not save the downloaded files in the Brainstorm folders (program or database folders) * '''MaxFilter/tSSS''': The MEG recordings from this study were recorded on an Elekta MEG, and need to be processed with cleaning algorithms SSS or tSSS. Brainstorm currently does not offer any free alternative to Elekta's MaxFilter, therefore in this tutorial we will import the recordings already processed with MaxFilter's tSSS, available in the "derivatives" folder of the BIDS architecture. To save disk space, will not download the raw MEG recordings. This means we will not get any of the additional files available in the BIDS structure (headshape, events, channels) but it doesn't matter because all the information is directly available in the .fif files. * Start Brainstorm (Matlab scripts or stand-alone version). For help, see the [[Installation]] page. * Select the menu File > Create new protocol. Name it "'''TutorialVisual'''" and select the options: |
| Line 23: | Line 70: |
| This dataset is formatted following the [[https://www.biorxiv.org/content/early/2017/08/08/172684|BIDS-MEG specifications]], therefore we could import all the relevant information (MRI, FreeSurfer segmentation, MEG+EEG recordings) in just one click, with the menu '''File > Load protocol''' > '''Import BIDS dataset''', as described in the online tutorial [[http://neuroimage.usc.edu/brainstorm/Tutorials/RestingOmega#BIDS_specifications|MEG resting state & OMEGA database]]. However, because your own data might not be organized following the BIDS standards, in this tutorial we preferred illustrating all the detailed steps for importing the data rather than the BIDS shortcut. This page explains how to import and process the first run of '''subject''' '''#01''' only. All the other files will have to be processed in the same way. |
|
| Line 24: | Line 75: |
| * Right-click on the TutorialAuditory folder > New subject > '''Subject01''' * Leave the default options you set for the protocol * Right-click on the subject node > Import anatomy folder: |
* Right-click on the TutorialVisual folder > New subject > '''sub-01''' * Leave the default options you defined for the protocol. * Right-click on the subject node > '''Import anatomy folder''': |
| Line 28: | Line 79: |
| * Select the folder: '''sample_auditory/anatomy''' | * Select the folder: '''derivatives/freesurfer/sub-01''' |
| Line 30: | Line 81: |
| * Set the 6 required fiducial points (indicated in MRI coordinates): * NAS: x=127, y=213, z=139 * LPA: x=52, y=113, z=96 * RPA: x=202, y=113, z=91 * AC: x=127, y=119, z=149 * PC: x=128, y=93, z=141 * IH: x=131, y=114, z=206 (anywhere on the midsagittal plane) * At the end of the process, make sure that the file "cortex_15000V" is selected (downsampled pial surface, that will be used for the source estimation). If it is not, double-click on it to select it as the default cortex surface. {{attachment:anatomy.gif||height="264",width="335"}} |
* When asked to select the anatomical fiducials, click on "'''Compute MNI transformation'''". This will register the MRI volume to an MNI template with an affine transformation, using SPM functions embedded in Brainstorm (spm_maff8.m). This will also place default fiducial points NAS/LPA/RPA in the MRI, based on MNI coordinates. We will use head points for the MEG-MRI coregistration, therefore we don't need accurate anatomical positions here. * Note that if you don't have a good digitized head shape for the subject, or if the final registration MEG-MRI doesn't look good with this head shape, you should not skip this step and mark accurate NAS/LPA/RPA fiducial points in the MRI, using the same anatomical convention that was used during the MRI acquisition. * Then click on [Save] to start the import. * At the end of the process, make sure that the file "cortex_15000V" is selected (downsampled pial surface, that will be used for the source estimation). If it is not, double-click on it to select it as the default cortex surface. Do not worry about the big holes in the head surface, parts of MRI have been remove voluntarily for anonymization purposes.<<BR>><<BR>> {{attachment:anatomy_import.gif||width="613",height="384"}} * All the anatomical atlases [[Tutorials/LabelFreeSurfer|generated by FreeSurfer]] were automatically imported: the cortical atlases (Desikan-Killiany, Mindboggle, Destrieux, Brodmann) and the sub-cortical regions (ASEG atlas). <<BR>><<BR>> {{attachment:anatomy_atlas.gif||width="550",height="211"}} |
| Line 42: | Line 91: |
| We need to attach the continuous .fif files containing the recordings to the database. |
|
| Line 43: | Line 94: |
| * Right-click on the subject folder > Review raw file * Select the file format: "'''MEG/EEG: CTF (*.ds...)'''" * Select all the .ds folders in: '''sample_auditory/data''' {{attachment:raw1.gif||height="156",width="423"}} * Refine registration now? '''YES'''<<BR>><<BR>> {{attachment:raw2.gif||height="224",width="353"}} === Multiple runs and head position === * The two AEF runs 01 and 02 were acquired successively, the position of the subject's head in the MEG helmet was estimated twice, once at the beginning of each run. The subject might have moved between the two runs. To evaluate visually the displacement between the two runs, select at the same time all the channel files you want to compare (the ones for run 01 and 02), right-click > Display sensors > MEG.<<BR>><<BR>> {{attachment:raw3.gif||height="220",width="441"}} * Typically, we would like to group the trials coming from multiple runs by experimental conditions. However, because of the subject's movements between runs, it's not possible to directly compare the sensor values between runs because they probably do not capture the brain activity coming from the same regions of the brain. * You have three options if you consider grouping information from multiple runs: * '''Method 1''': Process all the runs separately and average between runs at the source level: The more accurate option, but requires a lot more work, computation time and storage. * '''Method 2''': Ignore movements between runs: This can be acceptable for commodity if the displacements are really minimal, less accurate but much faster to process and easier to manipulate. * '''Method 3''': Co-register properly the runs using the process Standardize > Co-register MEG runs: Can be a good option for displacements under 2cm.<<BR>>Warning: This method has not be been fully evaluated on our side, to use at your own risk. Also, it does not work correctly if you have different SSP projectors calculated for multiple runs. * In this tutorial, we will illustrate only method 1: runs are not co-registered. === Epoched vs. continuous === * The CTF MEG system can save two types of files: epoched (.ds) or continuous (_AUX.ds). * Here we have an intermediate storage type: continuous recordings saved in an "epoched" file. The file is saved as small blocks of recordings of a constant time length (1 second in this case). All these time blocks are contiguous, there is no gap between them. * Brainstorm can consider this file either as a continuous or an epoched file. By default it imports the regular .ds folders as epoched, but we can change this manually, to process it as a continuous file. * Double-click on the "Link to raw file" for run 01 to view the MEG recordings. You can navigate in the file by blocks of 1s, and switch between blocks using the "Epoch" box in the Record tab. The events listed are relative to the current epoch.<<BR>><<BR>> {{attachment:raw4.gif||height="206",width="575"}} * Right-click on the "Link to raw file" for run 01 > '''Switch epoched/continuous''' * Double-click on the "Link to raw file" again. Now you can navigate in the file without interruptions. The box "Epoch" is disabled and all the events in the file are displayed at once.<<BR>><<BR>> {{attachment:raw5.gif||height="209",width="576"}} * Repeat this operation twice to convert all the files to a continuous mode. * '''Run 02''' > Switch epoched/continuous * '''Noise''' > Switch epoched/continuous == Stimulation triggers delay == === Evaluation === * Right-click on Run01/Link to raw file > '''Stim '''> Display time series (stimulus channel, UPPT001)<<BR>>Right-click on Run01/Link to raw file > '''ADC V''' > Display time series (audio signal generated, UADC001) * In the Record tab, set the duration of display window to '''0.200s'''.<<BR>>Jump to the third event in the "standard" category. * We can observe that there is a delay of about '''13ms''' between the time where the stimulus trigger is generated by the stimulation computer and the moment where the sound is actually played by the sound card of the stimulation computer ('''delay #1'''). This is matching the documentation of the experiment in the first section of this tutorial. <<BR>><<BR>> {{attachment:stim1.gif||height="282",width="590"}} === Correction === * '''Delay #1''': We can detect the triggers from the analog audio signal (ADC V/UADC001) rather than using the events already detected by the CTF software from the stim channel (Stim/UPPT001). * Drag and drop '''Run01 '''and '''Run02 '''to the Process1 box. * Add __'''twice'''__ the process "'''Events > Detect analog triggers'''".<<BR>>Once with event name="standard_fix" and reference event="standard".<<BR>>Once with event name="deviant_fix" and reference event="deviant".<<BR>>Set the other options as illustrated below:<<BR>><<BR>> {{attachment:stim2.gif}} * Open Run01 (channel ADC V) to evaluate the correction that was performed by this process. If you look at the third trigger in the "standard" category, you can measure a 14.6ms delay between the original event "standard" and the new event "standard_fix".<<BR>><<BR>> {{attachment:stim3.gif||height="151",width="570"}} * Open '''Run01''' to re-organize the event categories: * '''Delete '''the unused event categories: '''standard''', '''deviant'''. * '''Rename '''standard_fix and deviant_fix to '''standard''' and '''deviant'''. * Open '''Run02''' and do the same cleaning operations:<<BR>><<BR>> {{attachment:stim5.gif||height="149",width="502"}} * '''Important note''': We compensated for the jittered delays (delay #1), but not for the other ones (delays #2, #3 and #4). There is still a''' constant 5ms delay''' between the stimulus triggers ("standard" and "deviant") and the time where the sound actually reaches the subject's ears. == Detect and remove artifacts == |
* Right-click on the subject folder > '''Review raw file'''. * Select the file format: "'''MEG/EEG: Neuromag FIFF (*.fif)'''" * Go to the folder: '''derivatives/meg-derivatives/sub-01/ses-meg/meg''' * Select file: '''sub-01_ses-meg_task-facerecognition_run-01_proc-tsss_meg.fif''' <<BR>><<BR>> {{attachment:review_raw.gif||width="583",height="205"}} * Events:''' Ignore'''. We will read the stimulus triggers later.<<BR>><<BR>> {{attachment:review_ignore.gif||width="330",height="186"}} * Refine registration now? '''NO''' <<BR>>The head points that are available in the FIF files contain all the points that were digitized during the MEG acquisition, including the ones corresponding to the parts of the face that have been removed from the MRI. If we run the fitting algorithm, all the points around the nose will not match any close points on the head surface, leading to a wrong result. We will first remove the face points and then run the registration manually. === Channel classification === A few non-EEG channels are mixed in with the EEG channels, we need to change this before applying any operation on the EEG channels. * Right-click on the channel file > '''Edit channel file'''. Double-click on a cell to edit it. * Change the type of '''EEG062''' to '''EOG '''(electrooculogram). * Change the type of '''EEG063 '''to '''ECG '''(electrocardiogram). * Change the type of '''EEG061''' and '''EEG064''' to '''NOSIG''' (or any other non-informative type). Close the window and save the modifications. <<BR>><<BR>> {{attachment:channel_edit.gif||width="561",height="252"}} === MRI registration === At this point, the registration MEG/MRI is based only on the three anatomical landmarks NAS/LPA/RPA, which are not even accurately set (we used default MNI positions). All the MRI scans were anonymized (defaced) so all the head points below the nasion cannot be used. We will try to refine this registration using the additional head points that were digitized above the nasion. * Right-click on the channel file > '''Digitized head points > Remove points below nasion'''. <<BR>>(45 points removed, 92 head points left)<<BR>><<BR>> {{attachment:channel_remove.gif||width="329",height="194"}} * Right-click on the channel file > '''MRI registration > Refine using head points'''.<<BR>><<BR>> {{attachment:channel_refine.gif||width="294",height="173"}} * MEG/MRI registration, before (left) and after (right) this automatic registration procedure: <<BR>><<BR>> {{attachment:registration.gif||width="236",height="209"}} {{attachment:registration_final.gif||width="237",height="208"}} * Right-click on the channel file > '''MRI registration > EEG: Edit'''.<<BR>>Click on ['''Project electrodes on surface'''], then close the figure to save the modifications.<<BR>><<BR>> {{attachment:channel_project.gif||width="477",height="207"}} === Read stimulus triggers === We need to read the stimulus markers from the STI channels. The following tasks can be done in an interactive way with menus in the Record tab, as in the introduction tutorials. We will illustrate here how to do this with the pipeline editor, it will be easier to batch it for all the runs and all the subjects. * In Process1, select the "Link to raw file", click on [Run]. * Select process '''Events > Read from channel''', Channel: '''STI101''', Detection mode: '''Bit'''.<<BR>>Do not execute the process, we will add other processes to classify the markers.<<BR>><<BR>> {{attachment:raw_read_events.gif||width="499",height="290"}} * We want to create three categories of events, based on their numerical codes: * '''Famous faces''': 5 (00101), 6 (00110), 7 (00111) => Bit 3 only * '''Unfamiliar faces''': 13 (01101), 14 (01110), 15 (01111) => Bit 3 and 4 * '''Scrambled images''': 17 (10001), 18 (10010), 19 (10011) => Bit 5 only * We will start by creating the category "Unfamiliar" (combination of events "3" and "4") and remove the remove the initial events. Then we just have to rename the remaining "3" in "Famous", and all the "5" in "Scrambled". * Add process '''Events > Group by name''': "'''Unfamiliar=3,4'''", Delay=0, '''Delete original events''' * Add process '''Events > Rename event''': 3 => Famous * Add process '''Events > Rename event''': 5 => Scrambled <<BR>><<BR>> {{attachment:events_merge.gif||width="532",height="402"}} * Add process '''Events > Add time offset''' to correct for the presentation delays:<<BR>>Event names: "'''Famous, Unfamiliar, Scrambled'''", Time offset = '''34.5ms'''<<BR>><<BR>> {{attachment:events_offset.gif||width="348",height="355"}} * Finally run the script. Double-click on the recordings to make sure the labels were detected correctly. You can delete the unwanted events that are left in the recordings (1,2,9,13). Edit the colors if you don't like the default ones.<<BR>><<BR>> {{attachment:events_display.gif||width="595",height="194"}} == Pre-processing == |
| Line 89: | Line 137: |
| * One of the typical pre-processing steps consist in getting rid of the contamination due to the power lines (50 Hz or 60Hz). Let's start with the spectral evaluation of this file. * Drag '''ALL''' the "Link to raw file" to the Process1 box, or easier, just drag the node "Subject01", it will select recursively all the files in it. * Run the process "'''Frequency > Power spectrum density (Welch)'''": * Time window: '''[All file]''', Window length: '''4s''', Overlap: '''50%''', Sensor types: '''MEG''' * Note that you need at least 8Gb of RAM to run the PSD on the entire file. If you don't or if you get "Out of memory" errors, you can try running the PSD on a shorter time window. * Click on '''[Edit]''' and select option "'''Save individual PSD values''' (for each trial)".<<BR>><<BR>> {{attachment:psd1.gif||height="394",width="615"}} * Double-click on the new PSD files to display them. * Observations for '''Run''''''01''': <<BR>><<BR>> {{attachment:psd_eval01.gif||height="207",width="394"}} * Peaks related with the power lines: '''60Hz, 120Hz, 180Hz '''(240Hz and 300Hz could be observed as well depending on the window length used for the PSD) * The drop after '''600Hz''' corresponds to the low-pass filter applied at the acquisition time. * One channel indicates a higher noise than the others in high frequencies: '''MLO52''' (in red).<<BR>>We will probably mark it as bad later, when reviewing the recordings. * Observations for '''Run02''': <<BR>><<BR>> {{attachment:psd_eval02.gif||height="208",width="394"}} <<BR>> {{attachment:psd_eval02_zoom.gif||height="208",width="595"}} * Same peaks related with the power lines: '''60Hz, 120Hz, 180Hz''' * Same drop after '''600Hz'''. * Same noisy channel: '''MLO52'''. * Additionally, we observe higher level of noise in frequencies in the range of 30Hz to 100Hz on '''many occipital sensors'''. This is probably due to some tension in the neck due to an uncomfortable position. We will see later whether these channels need to be tagged as bad. === Power line contamination === * Put '''ALL''' the "Link to raw file" into the Process1 box (or directly the Subject01 folder) * Run the process: '''Pre-process > Notch filter''' * Select the frequencies: '''60, 120, 180 Hz''' * Sensor types or names: '''MEG''' * The higher harmonics are too high to bother us in this analysis, plus they are not clearly visible in all the recordings. * In output, this process creates new .ds folders in the same folder as the original files, and links the new files to the database.<<BR>><<BR>> {{attachment:psd3.gif}} * Run again the PSD process "'''Frequency > Power spectrum density (Welch)'''" on these new files, with the same parameters, to evaluate the quality of the correction. * Double-click on the new PSD files to open them.<<BR>><<BR>> {{attachment:psd5.gif||height="177",width="416"}} * Zoom in with the mouse wheel to observe what is happening around 60Hz (before / after).<<BR>><<BR>> {{attachment:psd6.gif||height="143",width="453"}} * To avoid the confusion later, delete the links to the original files: Select the folders containing the original unfiltered files and press the Delete key (or right-click > File > Delete).<<BR>><<BR>> {{attachment:psd7.gif||height="184",width="356"}} === Heartbeats and eye blinks === * Select the two AEF runs in the Process1 box. * Select successively the following processes, then click on [Run]: * '''Events > Detect heartbeats:''' Select channel '''ECG''', check "All file", event name "cardiac". * '''Events > Detect eye blinks:''' Select channel '''VEOG''', check "All file", event name "blink". * '''Events > Remove simultaneous''': Remove "'''cardiac'''", too close to "'''blink'''", delay '''250ms'''. * '''Compute SSP: Heartbeats''': Event name "cardiac", sensors="MEG", '''do not use existing SSP'''. * '''Compute SSP: Eye blinks''': Event name "blink", sensors="MEG", '''do not use existing SSP'''.<<BR>><<BR>> {{attachment:ssp_pipeline.gif||height="359",width="527"}} * Double-click on '''Run01 '''to open the MEG.<<BR>>You can change the color of events "standard" and "deviant" to make the figure more readable. * Review the '''EOG '''and '''ECG '''channels and make sure the events detected make sense.<<BR>><<BR>> {{attachment:events.gif||height="369",width="557"}} * In the Record tab, menu SSP > '''Select active projectors'''. * Blink: The first component is selected and looks good. * Cardiac: The category is disabled because no component has a value superior to 12%. * Select the first component of the cardiac category and display its topography. * It looks exactly like a cardiac topography, keep it selected and click on [Save].<<BR>><<BR>> {{attachment:ssp_result.gif||height="172",width="700"}} * Repeat the same operations for '''Run02''': * Review the events. * Select the first cardiac component. === Bad segments === * At this point, you should review the entire files, by pages of a few seconds scrolling with the F3 key, to identify all the bad channels and the noisy segments of recordings. Do this with the the EOG channel open at the same time to identify saccades or blinks that were not completely corrected with the SSP projectors. As this is a complicated task that requires some expertise, we have prepared a list of bad segments for these datasets. * Open '''Run01'''. In the Record tab, select '''File > Add events from file''': * File name: sample_auditory/data/S01_AEF_20131218_01_notch/'''events_bad_01.mat ''' * File type: Brainstorm (events*.mat) * It adds '''12 bad segments''' to the file. * Open '''Run02'''. In the Record tab, select '''File > Add events from file''': * File name: sample_auditory/data/S01_AEF_20131218_02_notch/'''events_bad_02.mat ''' * File type: Brainstorm (events*.mat) * It adds '''9 bad segments''' and '''16 saccades''' to the file. === Saccades === * Run02 contains a few saccades that generate a large amount of noise in the MEG recordings. They are not identified well by the automatic detection process based on the horizontal EOG. We have marked some of them, you have already loaded these events together with the bad segments. We are going to use again the SSP technique to remove the spatial components associated with these saccades. * Open the MEG recordings for '''Run02''' and select the right-frontal sensors (Record tab > CTF RF). * In the Record tab, menu SSP > Compute SSP: Generic<<BR>>Event name='''saccade''', Time='''[0,500]ms''', Frequency='''[1,15]Hz''', '''Use existing SSP'''<<BR>><<BR>> {{attachment:ssp_saccade_process.gif||height="442",width="329"}} * Example of saccade without correction:<<BR>><<BR>> {{attachment:ssp_saccade_before.gif||height="259",width="711"}} * With the first component of saccade SSP applied: <<BR>><<BR>> {{attachment:ssp_saccade_after.gif||height="283",width="507"}} * This first component removes really well the saccade, keep it selected and click on [Save]. === Bad channels === * During the visual exploration, some channels appeared generally noisier than the others. Example: <<BR>><<BR>> {{attachment:badchannel02.gif||height="200",width="456"}} * Right-click on '''Run01 '''> Good/bad channels > Mark some channels as bad<<BR>> > '''MRT51, MLO52''' * Right-click on '''Run02 '''> Good/bad channels > Mark some channels as bad<<BR>> > '''MRT51, ''''''MLO52, MLO42, MLO43''' |
* Keep the "Link to raw file" in Process1. * Run process '''Frequency > Power spectrum density (Welch)''' with the options illustrated below. <<BR>><<BR>> {{attachment:psd_process.gif||width="669",height="314"}} * Right-click on the PSD file > Power spectrum. <<BR>><<BR>> {{attachment:psd_plot.gif||width="408",height="235"}} * The MEG sensors look awful, because of one small segment of data located around 248s. Open the MEG recordings and scroll to 248s (just before the first Unfamiliar event). <<BR>><<BR>> {{attachment:psd_error.gif||width="527",height="183"}} * In these recordings, the continuous head tracking was activated, but it starts only at the time of the stimulation (248s) while the acquisition of the data starts 20s before (226s). The first 20s do not have head localization coils (HPI) coils activity and are not corrected by MaxFilter. After 248s, the HPI coils are on, and MaxFilter filters them out. The transition between the two states is not smooth and creates important distortions in the spectral domain. For a proper evaluation of the recordings, we should '''compute the PSD only after the HPI coils are turned on'''. * Run process '''Events > Detect cHPI activity (Elekta)'''. This detects the changes in the cHPI activity from channel STI201 and marks all the data without head localization as bad. <<BR>><<BR>> {{attachment:psd_badsegment.gif||width="503",height="234"}} * Re-run the process '''Frequency > Power spectrum density (Welch)'''. All the bad segments are excluded from the computation, therefore the PSD is now estimated only with the data '''after 248s'''. <<BR>><<BR>> {{attachment:psd_fix.gif||width="668",height="283"}} * Observations: * Three groups of sensors, from top to bottom: EEG, MEG gradiometers, MEG magnetometers. * Power lines: '''50 '''Hz and harmonics * Alpha peak around 10 Hz * Artifacts due to Elekta electronics (HPI coils): '''293'''Hz, '''307'''Hz, '''314'''Hz, '''321'''Hz, '''328'''Hz. * Peak from unknown source at '''103.4Hz''' in the MEG only. * Suspected bad EEG channels: '''EEG016''' * Close all the windows. === Remove line noise === * Keep the "Link to raw file" in Process1. * Select process '''Pre-process > Notch filter''' to remove the line noise (50-200Hz).<<BR>>Add the process '''Frequency > Power spectrum density (Welch)'''. <<BR>><<BR>> {{attachment:notch_process.gif||width="562",height="262"}} * Double-click on the PSD for the new continuous file to evaluate the quality of the correction. <<BR>><<BR>> {{attachment:notch_result.gif||width="566",height="215"}} * Close all the windows (use the [X] button at the top-right corner of the Brainstorm window). === EEG reference and bad channels === * Right-click on link to the processed file ("Raw | notch(50Hz ...") > '''EEG > Display time series'''. * Select channel '''EEG016''' and mark it as '''bad''' (using the popup menu or pressing the Delete key). <<BR>><<BR>> {{attachment:channel_bad.gif||width="587",height="214"}} * In the Record tab, menu '''Artifacts > Re-reference EEG''' > "AVERAGE". <<BR>><<BR>> {{attachment:channel_ref.gif||width="535",height="270"}} * At the end, the window "select active projectors" is open to show the new re-referencing projector. Just close this window. To get it back, use the menu Artifacts > Select active projectors. == Artifact detection == === Heartbeats: Detection === * Empty the Process1 list (right-click > Clear list). * Drag and drop the continuous processed file ("Raw | notch(50Hz...)") to the Process1 list. * Run process''' Events > Detect heartbeats''': Channel name='''EEG063''', All file, Event name=cardiac <<BR>><<BR>> {{attachment:detect_cardiac.gif||width="310",height="230"}} === Eye blinks: Detection === * In many of the other tutorials, we detect the blinks and remove them with SSP. In this experiment, we are particularly interested in the subject's response to seeing the stimulus. Therefore we will exclude from the analysis all the recordings contaminated with blinks or other eye movements. * Run process '''Events > Detect events above threshold''': <<BR>>Event name='''blink_BAD''', Channel='''EEG062''', All file, Maximum threshold='''100''', Threshold units='''uV''', Filter='''[0.30,20.00]Hz''', '''Use absolute value of signal'''.<<BR>><<BR>> {{attachment:detect_blinks.gif||width="285",height="455"}} * In many other tutorials, we use the process "Detect eye blinks", which creates simple events indicating the peak of the EOG signal during a blink. In this study, we preferred using this process "Detect events above threshold" because it creates extended events marking as bad all the segment for which the EOG value is above a given threshold. It's more manual but it makes sure we really exclude all the ocular activity. * Inspect visually the two new cateogries of events: cardiac and blink_BAD.<<BR>><<BR>> {{attachment:detect_display.gif||width="655",height="237"}} * Close all the windows (using the [X] button). === Heartbeats: Correction with SSP === * Keep the "Raw | notch" file selected in Process1. * Select process '''Artifacts > SSP: Heartbeats''' > Sensor type: '''MEG MAG''' * Add process '''Artifacts > SSP: Heartbeats''' > Sensor type: '''MEG GRAD''' - Run the execution <<BR>><<BR>> {{attachment:ssp_ecg_process.gif||width="293",height="225"}} * Double-click on the continuous file to show all the MEG sensors. <<BR>>In the Record tab, select sensors "Left-temporal". * Menu Artifacts > Select active projectors. * In category '''cardiac/MEG MAG''': Select '''component #1''' and view topography. * In category '''cardiac/MEG GRAD''': Select '''component #1 '''and view topography. <<BR>><<BR>> {{attachment:ssp_ecg_topo.gif||width="604",height="175"}} * Make sure that selecting the two components removes the cardiac artifact. Then click '''[Save]'''. === Additional bad segments === * Process1>Run: Select the process "'''Events > Detect other artifacts'''". This should be done separately for MEG and EEG to avoid confusion about which sensors are involved in the artifact. <<BR>><<BR>> {{attachment:detect_other.gif||width="486",height="253"}} * Display the MEG sensors. Review the segments that were tagged as artifact, determine if each event represents an artifact and then mark the time of the artifact as BAD. This can be done by selecting the time window around the artifact, then right-click > '''Reject time segment'''. Note that this detection process marks 1-second segments but the artifact can be shorter.<<BR>><<BR>> {{attachment:artifacts_mark_bad.png||width="640"}} * Once all the events in the two categories are reviewed and bad segments are marked, the two categories (1-7Hz and 40-240Hz) can be deleted. * Do this detection and review again for the EEG. === SQUID jumps === MEG signals recorded with Elekta-Neuromag systems frequently contain SQUID jumps ([[http://neuroimage.usc.edu/brainstorm/Tutorials/BadSegments?highlight=(squid)#Elekta-Neuromag_SQUID_jumps|more information]]). These sharp steps followed by a change of baseline value are easy to identify visually but more complicated to detect automatically. The process "Detect other artifacts" usually detects most of them in the category "1-7Hz". If you observe that some are skipped, you can try re-running it with a higher sensitivity. It is important to review all the sensors and all the time in each run to be sure these events are marked as bad segments. . {{attachment:detect_jumps.gif||width="547",height="214"}} |
| Line 165: | Line 208: |
| === Import recordings === To import epochs from '''Run01''': * Right-click on the "Link to raw file" > '''Import in database''' * Use events: "'''standard'''" and "'''deviant'''" * Epoch time: '''[-100, +500] ms''' * Apply the existing SSP (make sure that you have 2 selected projectors) * '''Remove DC''' '''offset '''based on time window: '''[-100, 0] ms''' * '''UNCHECK''' the option "Create a separate folder for each epoch type", this way all the epochs are going to be saved in the same Run01 folder, and we will able to separate the trials from Run01 and Run02.<<BR>><<BR>> {{attachment:import1.gif||height="344",width="521"}} * Note that the trials that are overlapping with a BAD segment are tagged as bad in the database explorer (marked with a red dot).<<BR>><<BR>> {{attachment:import2.gif}} Repeat the same operation for '''Run02''': * Right-click on the "Link to raw file" > '''Import in database''' * Use the same options as for the previous run.<<BR>><<BR>> {{attachment:import3.gif||height="235",width="216"}} === Average responses === * As said previously, it is usually not recommended to average recordings in sensor space across multiple acquisition runs because the subject might have moved between the sessions. Different head positions were recorded for each run, we will reconstruct the sources separately for each each run to take into account these movements. * However, in the case of event-related studies it makes sense to start our data exploration with an average across runs, just to evaluate the quality of the evoked responses. We have seen that the subject almost didn't move between the two runs, so the error would be minimal. We will compute now an approximate sensor average between runs, and we will run a more formal average in source space later. * We have 80 good "deviant" trials that we want to average together. * Select the trial groups "deviant" from both runs in Process1, run process "'''Average > Average files'''"<<BR>>Select the option "'''By trial group (subject average)'''"<<BR>><<BR>> {{attachment:process_average_data1.gif||height="456",width="470"}} * To compare properly this "deviant" average with the other condition, we need to use the same number of trials in the "standard" condition. We are going to pick 40 "standard" trials from Run01 and 40 from Run02. To make it easy, let's take the 40 first good trials. * Select the '''41 '''first "standard" trials of Run01 + the '''41 '''first "standard" trials of Run02 in Process1.<<BR>>This will sum to '''80 '''selected files, because the Process1 tab ignores the bad trials (trial #37 is bad in Run01, trial #36 is bad in Run02) * Run again process "'''Average > Average files'''" > "'''By trial group (subject average)'''"<<BR>><<BR>> {{attachment:process_average_data2.gif||height="444",width="471"}} * The average for the two conditions "standard" and "deviant" are saved in the folder ''(intra-subject)''. The channel file added to this folder is an average of the channel files from Run01 and Run02.<<BR>><<BR>> {{attachment:average_sensor_files.gif||height="205",width="220"}} === Visual exploration === * Display the two averages, "standard" and "deviant": * Right-click on average > MEG > Display time series * Right-click on average > MISC > Display time series (EEG electrodes Cz and Pz) * Right-click on average > MEG > 2D Sensor cap * In the Filter tab, add a '''low-pass filter''' at '''100Hz'''. * Right-click on the 2D topography figures > Snapshot > Time contact sheet. * Here are results for the standard (top) and deviant (bottom) beeps:<<BR>><<BR>> {{attachment:average_sensor.gif||height="399",width="703"}} * '''P50''': 50ms, bilateral auditory response in both conditions. * '''N100''': 95ms, bilateral auditory response in both conditions. * '''MMN''': 100-200ms, mismatch negativity in the deviant condition only (detection of deviant). * '''P200''': 170ms, in both conditions but much stronger in the standard condition. * '''P300''': 300-400ms, deviant condition only (decision making in preparation of the button press). * '''Standard '''(right-click on the topography figure > Snapshot > Time contact sheet) : <<BR>><<BR>> {{attachment:average_sensor_standard.gif||height="286",width="370"}} * '''Deviant''': <<BR>><<BR>> {{attachment:average_sensor_deviant.gif||height="300",width="371"}} === Difference deviant-standard === * In the Process2 tab, select the deviant average (Files A) and the standard average (Files B). * Run the process "'''Other > Difference A-B'''"<<BR>><<BR>> {{attachment:average_sensor_diff.gif||height="293",width="340"}} * The difference deviant-standard does not show anymore the early responses (P50, P100) but emphasizes the difference in the later process (MMN/P200 and P300). <<BR>><<BR>> {{attachment:average_sensor_diff2.gif||height="214",width="456"}} == Source estimation == === Head model === * Select the two imported folders at once, right-click > Compute head model<<BR>><<BR>> {{attachment:headmodel1.gif||height="221",width="413"}} * Use the '''overlapping spheres''' model and keep all of the options at their default values.<<BR>><<BR>> {{attachment:headmodel2.gif||height="234",width="209"}} {{attachment:headmodel3.gif||height="208",width="244"}} * For more information: [[Tutorials/TutHeadModel|Head model tutorial]]. === Noise covariance matrix === * We want to calculate the noise covariance from the empty room measurements and use it for the other runs. * In the '''Noise''' folder, right-click on the Link to raw file > Noise covariance > Compute from recordings.<<BR>><<BR>> {{attachment:noisecov1.gif||height="253",width="392"}} * Keep all the default options and click [OK]. <<BR>><<BR>> {{attachment:noisecov2.gif||height="285",width="291"}} * Right-click on the noise covariance file > Copy to other conditions.<<BR>><<BR>> {{attachment:noisecov3.gif||height="181",width="232"}} * You can double-click on one the copied noise covariance files to check what it looks like:<<BR>><<BR>> {{attachment:noisecov4.gif||height="232",width="201"}} * For more information: [[Tutorials/TutNoiseCov|Noise covariance tutorial]]. === Inverse model === * Select the two imported folders at once, right-click > Compute sources<<BR>><<BR>> {{attachment:inverse1.gif||height="192",width="282"}} * Select '''dSPM''' and keep all the default options.<<BR>><<BR>> {{attachment:inverse2.gif||height="205",width="199"}} * Then you are asked to confirm the list of bad channels to use in the source estimation for each run. Just leave the defaults, which are the channels that we set as bad earlier.<<BR>><<BR>> {{attachment:inverse_bad1.gif||height="202",width="326"}} * One inverse operator is created in each condition, with one link per data file.<<BR>><<BR>> {{attachment:inverse3.gif||height="240",width="202"}} * For more information: [[Tutorials/TutSourceEstimation|Source estimation tutorial]]. === Average in source space === * Now we have the source maps available for all the trials, we average them in source space. * Select the folders for '''Run01 '''and '''Run02 '''and the ['''Process sources'''] button on the left. * Run process "'''Average > Average files'''":<<BR>>Select "'''By trial group (subject average)'''"<<BR>><<BR>> {{attachment:process_average_results.gif||height="376",width="425"}} * Double-click on the source averages to display them (standard=top, deviant=bottom).<<BR>><<BR>> {{attachment:average_source.gif||height="321",width="780"}} * Note that opening the source maps can be very long because of the online filters. Check in the Filter tab, you probably still have a '''100Hz low-pass filter''' applied for the visualization. In the case of averaged source maps, the 15000 source signals are filtered on the fly when you load a source file. This can take a significant amount of time. You may consider unchecking this option if the display is too slow on your computer. * '''Standard:''' (Right-click on the 3D figures > Snapshot > Time contact sheet)<<BR>><<BR>> {{attachment:average_source_standard_left.gif||height="263",width="486"}} <<BR>> {{attachment:average_source_standard_right.gif||height="263",width="486"}} * '''Deviant:'''<<BR>><<BR>> {{attachment:average_source_deviant_left.gif||height="263",width="486"}} <<BR>> {{attachment:average_source_deviant_right.gif||height="263",width="486"}} * '''Movies''': Right-click on any figure > Snapshot > '''Movie (time): All figures''' (click to download video)<<BR>><<BR>>[[http://neuroimage.usc.edu/wikidocs/average_sources.avi|{{attachment:average_source_video.gif|http://neuroimage.usc.edu/brainstorm/wikidocs/average_sources.avi|height="258",width="484"}}]] === Difference deviant-standard === * In the Process2 tab, select the average of the deviant sources (Files A) and the average of the standard sources (Files B). * Run the process "'''Other > Difference A-B'''"<<BR>><<BR>> {{attachment:average_source_diff.gif||height="347",width="424"}} * Double-click on the difference to display it, explore it in time.<<BR>><<BR>> {{attachment:average_source_diff_left.gif||height="252",width="455"}} <<BR>> {{attachment:average_source_diff_right.gif||height="252",width="455"}} * The first observations we can make are the following: * '''P50''': No important difference. * '''N100''': Stronger response in the right auditory system for the deviant condition. * '''MMN '''(125ms): Stronger response for the deviant (left auditory, right temporal/frontal/motor). * '''P200 '''(175ms): Stronger response in the auditory system for the standard condition. * '''After 200ms''': Stronger response in the deviant condition (left auditory, left motor, right auditory, right temporal, right motor, right parietal) * Alternatively, you could calculate the difference of the average of all the "deviant" (80) and all the "standard" (388) trials, using the process "'''Other > Weighted difference'''". It is an attempt to compensate for the difference of number of trials. === Student's t-test === * Using a t-test instead of the difference of the two averages, you can reproduce similar results but with a significance level attached to each value. With this test, we can also use all the trials we have, unlike the difference of the means: The t-test behaves very well with imbalanced designs, we can keep all the standard trials. * In the Process2 tab, select the following files: * Files A: All the deviant trials, with the '''[Process sources]''' button selected. * Files B: All the standard trials, with the '''[Process sources]''' button selected. * Run the process "'''Test > Student's t-test'''", Equal variance, Absolute value of average.<<BR>><<BR>> {{attachment:ttest_source1.gif||height="406",width="699"}} * Double-click on the t-test file to open it. Set the options in the Stat tab: * p-value threshold: '''0.05''' * Multiple comparisons: '''FDR''' * Control over dimensions: '''1.Signals''' and '''2.Time''' * Explore the results in time.<<BR>><<BR>> {{attachment:ttest_source2_left.gif||height="231",width="415"}} <<BR>> {{attachment:ttest_source2_right.gif||height="231",width="415"}} == Regions of interest == === Manual tracing === * Let's place all the regions of interest starting from the easiest to identify. * Open the average source files (standard and deviant), together with the average recordings for the standard condition for an easier time navigation. * In the Surface tab, smooth the cortical surface at '''70%'''. * For each region: go to the indicated time point, adjust the amplitude threshold in the Surface tab, identify the area of interest, click on its center, grow the scout, rename it. * Grow all the regions to the same size: '''20 vertices'''. * Note that all the following screen captures are produced with a online low-pass filter at 100Hz. * '''A1L''': Left primary auditory cortex (Heschl gyrus) * The most visible region in both conditions. Active during all the main steps of the auditory processing: P50, N100, MMN, P200, P300. * '''Standard '''condition, t='''90ms''', amplitude threshold='''70%'''<<BR>><<BR>> {{attachment:scout_a1l.gif||height="162",width="650"}} * '''A1R''': Right primary auditory cortex (Heschl gyrus) * The position of this region is a lot less obvious than A1L, we don't see one focal region with a sustained activity. These binaural auditory stimulations should be generating similar bilateral responses in both left and right auditory cortices at early latencies. Possible explanations for this observation: * The earplug was not adjusted on the right side and the sound was not well delivered. * The subject's hearing from the right ear is impaired. * The response is actually stronger in the left auditory cortex for this subject. * The orientation of the source makes it difficult to capture for the MEG sensors. * We are trying to find a region that peaks at the same time as A1L (95ms and 200ms in the standard condition). It is very difficult to find anything that behaves this way in both the deviant and the standard condition, so we will pick something very approximate, knowing that we cannot really rely on this region. The auditory system is very dynamic, in squared centimeters of cortex we can observe many functionally independent regions activated at different moments. * '''Deviant '''condition, t='''34ms''', amplitude threshold='''20%''' * Results are ok for the deviant condition but not so good for the standard condition. <<BR>><<BR>> {{attachment:scout_a1r.gif||height="162",width="648"}} * '''IFGL''': Left inferior frontal gyrus (Brodmann area 44) * Involved in the auditory processing, particularly while processing irregularities. * You can use the atlas "Brodmann-thresh" available in the Scout tab for identifying this region. * '''Deviant '''condition, t='''140ms''', amplitude threshold='''40% '''<<BR>><<BR>> {{attachment:scout_ifgl.gif||height="162",width="648"}} * '''IFGR''': Right inferior frontal gyrus (Brodmann area 44''')''' * Expected to have an activity similar to the left IFG. * '''Deviant '''condition, t='''110ms''', amplitude threshold='''40% '''<<BR>><<BR>> {{attachment:scout_ifgr.gif||height="162",width="648"}} * '''M1L''': Left motor cortex * The subject taps with the right index when a deviant is presented. * The motor cortex responds at very early latencies together with the auditory cortex, in both conditions (50ms and 100ms). The subject is ready for a fast response to the task. * At 175ms, the peak in the deviant condition probably corresponds to an inhibition: the sound heard is not a deviant, there is no further motor processing required. * At 225ms, the peak in the standard condition is probably a motor preparation. At 350ms, the motor task begins, the subject moves the right hand (recorded reaction times 500ms +/- 108ms). * '''Deviant '''condition, t='''240ms''', amplitude threshold='''50% '''<<BR>><<BR>> {{attachment:scout_m1l.gif||height="162",width="648"}} * '''M1R''': Right motor cortex * Probably involved in the preparation of the motor response as well. Less recruited during the actual motor command. * '''Deviant '''condition, t='''35ms''', amplitude threshold='''25%'''''' '''<<BR>><<BR>> {{attachment:scout_m1r.gif||height="162",width="648"}} * '''PPCR''': Right posterior parietal cortex * Known to play a role as a relay in the auditory processing. * '''Deviant '''condition, t='''225ms''', amplitude threshold='''60%'''''' '''<<BR>><<BR>> {{attachment:scout_ppcr.gif||height="163",width="652"}} === Latency comparisons [TODO] === === Influence of the number of trials === * We have decided to run the source analysis on the same number of trials for both conditions. We have been working so far with an average of the standard condition calculated from 80 trials. Just out of curiosity, we can recalculate another average with all the good standard trials (388). * Here are the scouts traces for both averages (80 trials in green, 388 trials in red):<<BR>><<BR>> {{attachment:scouts_ntrials.gif||height="266",width="607"}} * As expected, the signal is cleaner in the average with more trials, but it is interesting to note that the overall shape of the traces does not change. The main effects observed are similar, the latencies are identical, multiplying the number of trials by five does not change much the interpretation. == Time-frequency [TODO] == * Select in Process1 all the trials for the deviant condition (80 files) and the 41 first trials of each run for the deviant condition (80 files), just like you did when calculating the averages. * Select the ['''Process sources'''] button. * Run the process "'''Frequency > Time-frequency (Morlet wavelets)'''" * This produced one time-frequency file per trial in input: 80 for each condition. Now we want to average the time-frequency maps of all the trials for each condition separately. * Select the ['''Process time-freq'''] button. * Run the process "'''Average > Average files'''" > "'''By trial group (subject average)'''" * Z-score + Difference of the averages * t-test == Coherence [TODO] == == Dipole scanning [TODO] == == Phase-amplitude coupling [TODO] == == Discussion == Discussion about the choice of the dataset on the FieldTrip bug tracker:<<BR>>http://bugzilla.fcdonders.nl/show_bug.cgi?id=2300 == Scripting [TODO] == ==== Process selection ==== ==== Graphic edition ==== ==== Generate Matlab script ==== The operations described in this tutorial can be reproduced from a Matlab script, available in the Brainstorm distribution: '''brainstorm3/toolbox/script/tutorial_auditory.m ''' <<EmbedContent(http://neuroimage.usc.edu/bst/get_feedback.php?Tutorials/Auditory)>> |
=== Import epochs === * Keep the "Raw | notch" file selected in Process1. * Select process: '''Import recordings > Import MEG/EEG: Events''' (do not run immediately)<<BR>>Event names "Famous, Unfamiliar, Scrambled", All file, Epoch time=[-500,1200]ms * Add process: '''Pre-process > Remove DC offset''': Baseline=[-500,-0.9]ms - Run execution.<<BR>><<BR>> {{attachment:import_epochs.gif||width="647",height="399"}} === Average by run === * In Process1, select all the imported trials. * Run process: '''Average > Average files''': '''By trial groups (folder average)''' <<BR>><<BR>> {{attachment:average_process.gif||width="478",height="504"}} === Review EEG ERP === * EEG evoked response (famous, scrambled, unfamiliar), low-pass filter at 32Hz: <<BR>><<BR>> {{attachment:average_topo.gif||width="684",height="392"}} * Open the [[Tutorials/ChannelClusters|Cluster tab]] and create a cluster with the channel '''EEG065''' (button [NEW IND]).<<BR>><<BR>> {{attachment:cluster_create.gif||width="569",height="166"}} * Select the cluster, select the three average files, right-click > '''Clusters time series''' (Overlay:Files). <<BR>><<BR>> {{attachment:cluster_display.gif||width="658",height="220"}} * Basic observations for EEG065 (right parieto-occipital electrode): * Around 170ms (N170): greater negative deflection for Famous than Scrambled faces. * After 250ms: difference between Famous and Unfamiliar faces. == Source estimation [TODO] == === MEG noise covariance: Empty room recordings [TODO] === The minimum norm model we will use next to estimate the source activity can be improved by modeling the the noise contaminating the data. The introduction tutorials explain how to [[Tutorials/NoiseCovariance|estimate the noise covariance]] in different ways for EEG and MEG. For the MEG recordings we will use the empty room measurements we have, and for the EEG we will compute it from the pre-stimulus baselines we have in all the imported epochs. * Create a new subject: '''emptyroom''' * Right-click on the new subject > '''Review raw file'''.<<BR>>Select file: '''sub-emptyroom/ses-meg/meg/090707_raw_st.fif '''(processed with MaxFilter/tSSS)<<BR>>Do not apply default transformation, Ignore event channel.<<BR>><<BR>> {{attachment:noise_review.gif||width="506",height="164"}} * Select this new file in Process1 and run process '''Pre-process > Notch filter''': '''50 100 150 200Hz'''. When using empty room measurements to compute the noise covariance, they must be processed exactly in the same way as the other recordings (same MaxFilter parameters, same frequency filters). <<BR>><<BR>> {{attachment:noise_notch.gif||width="436",height="225"}} * Right-click on the filtered noise recordings > '''Noise covariance > Compute from recordings''': <<BR>><<BR>> {{attachment:noise_compute.gif||width="649",height="240"}} * Right-click on the '''Noise covariance > Copy to other subjects '''<<BR>><<BR>> {{attachment:noise_copy.gif||width="273",height="126"}} === EEG noise covariance: Pre-stimulus baseline === * In folder sub-01/run-01_tsss_notch, select all the imported the imported trials, right-click > '''Noise covariance > Compute from recordings''', Time=[-500,-0.9]ms, '''EEG only''', '''Merge'''. {{attachment:noise_eeg.gif||width="650",height="280"}} * This computes the noise covariance only for EEG, and combines it with the existing MEG information. {{attachment:noise_display.gif||width="375",height="223"}} === BEM layers === We will compute a [[Tutorials/TutBem|BEM forward model]] to estimate the brain sources from the EEG recordings. For this, we need some layers defining the separation between the different tissues of the head (scalp, inner skull, outer skull). * Go to the anatomy view (first button above the database explorer). * Right-click on the subject folder > '''Generate BEM surfaces''': The number of vertices to use for each layer depends on your computing power and the accuracy you expect. You can try for instance with '''1082 vertices '''(scalp) and '''642 vertices '''(outer skull and inner skull). <<BR>><<BR>> {{attachment:anatomy_bem.gif||width="603",height="232"}} === Forward model: EEG and MEG === * Go back to the functional view (second button above the database explorer). * Model used: Overlapping spheres for MEG, OpenMEEG BEM for EEG ([[Tutorials/HeadModel|more information]]). * In folder sub-01/run_01_tsss_notch, right-click on the channel file > '''Compute head model'''.<<BR>>Keep all the default options . Expect this to take a while...<<BR>><<BR>> {{attachment:headmodel_compute.gif||width="643",height="313"}} === Inverse model: Minimum norm estimates === * Right-click on the new head model > '''Compute sources [2016]''': '''MEG MAG + GRAD''' (default options) <<BR>><<BR>> {{attachment:sources_compute.gif||width="626",height="406"}} * Right-click on the new head model > '''Compute sources [2016]''': '''EEG''' (default bad channels). * At the end we have two inverse operators, that are shared for all the files of the run (single trials and averages). If we wanted to look at the run-level source averages, we could normalize the source maps with a Z-score wrt baseline. In this tutorial, we will first average across runs and normalize the subject-level averages. This will be done in the next tutorial (group analysis).<<BR>><<BR>> {{attachment:sources_files.gif||width="243",height="147"}} == Time-frequency analysis == We will compute the time-frequency decomposition of each trial using Morlet wavelets, and average the power of the Morlet coefficients for each condition and each run separately. We will restrict the computation to the MEG magnetometers and the EEG channels to limit the computation time and disk usage. * In Process1, select the imported trials "Famous" for run#01. * Run process '''Frequency > Time-frequency (Morlet wavelets)''': Sensor types=MEG MAG,EEG <<BR>>Not normalized, Frequency='''Log(6:20:60)''', Measure=Power, '''Save average''' <<BR>><<BR>> {{attachment:tf_process.gif||width="628",height="479"}} * Double-click on the file to display it. In the Display tab, select the option "Hide edge effects" to exclude from the display all the values that could not be estimated in a reliable way. Let's extract only the good values from this file (-200ms to +900ms). <<BR>><<BR>> {{attachment:tf_display.gif||width="270",height="155"}} * In Process1, select the time-frequency file. * Run process '''Extract > Extract time''': Time window='''[-200, 900]ms''', '''Overwrite input files''' <<BR>><<BR>> {{attachment:tf_cut.gif||width="411",height="197"}} * Display the file again, observe that all the possibly bad values are gone. <<BR>><<BR>> {{attachment:tf_cutdisplay.gif||width="410",height="178"}} * You can display all the sensors at once (MEG MAG or EEG): right-click > 2D Layout (maps).<<BR>><<BR>> {{attachment:tf_2dlayout.gif||width="653",height="270"}} * Repeat these steps for other conditions (Scrambled and Unfamiliar) and the other runs (2-6). There is no way with this process to compute all the averages at once, as it is with the process "Average files". This will be easier to run from a script. * If we wanted to look at the run-level source averages, we could normalize these time-frequency maps. In this tutorial, we will first average across runs and normalize the subject-level averages. This will be done in the next tutorial (group analysis). == Scripting == We have now all the files we need for the group analysis ([[Tutorials/VisualGroup|next tutorial]]). We need to repeat the same operations for all the runs and all the subjects. Some of these steps are '''fully automatic''' and take a lot of time (filtering, computing the forward model), they should be executed from a script. However, we recommend you always '''review manually''' some of the pre-processing steps (selection of the bad segments and bad channels, SSP/ICA components). Do not trust blindly any fully automated cleaning procedure. For the strict reproducibility of this analysis, we provide a script that processes all the 19 subjects: '''brainstorm3/toolbox/script/tutorial_visual_single.m''' (execution time: 10-30 hours)<<BR>>Report for the first subject: [[http://neuroimage.usc.edu/bst/examples/report_TutorialVisual_sub001.html|report_TutorialVisual_sub001.html]] <<HTML(<div style="border:1px solid black; background-color:#EEEEFF; width:720px; height:500px; overflow:scroll; padding:10px; font-family: Consolas,Menlo,Monaco,Lucida Console,Liberation Mono,DejaVu Sans Mono,Bitstream Vera Sans Mono,Courier New,monospace,sans-serif; font-size: 13px; white-space: pre;">)>><<EmbedContent("http://neuroimage.usc.edu/bst/viewcode.php?f=tutorial_visual_single.m")>><<HTML(</div >)>> <<BR>>You should note that this is not the result of a fully automated procedure. The bad channels were identified manually and are defined for each run in the script. The bad segments were detected automatically, confirmed manually for each run and saved in external files, then exported as text ans copied at the end of this script. All the process calls (bst_process) were generated automatically using with the script generator (menu '''Generate .m script''' in the pipeline editor). Everything else was added manually (loops, bad channels, file copies). <<EmbedContent("http://neuroimage.usc.edu/bst/get_prevnext.php?prev=&next=Tutorials/VisualGroup")>> <<EmbedContent(http://neuroimage.usc.edu/bst/get_feedback.php?Tutorials/VisualSingle)>> |
MEG visual tutorial: Single subject (BIDS)
[WORK IN PROGRESS: THIS TUTORIAL IS NOT READY FOR PUBLIC USE]
Authors: Francois Tadel, Elizabeth Bock.
The aim of this tutorial is to reproduce in the Brainstorm environment the analysis described in the SPM tutorial "Multimodal, Multisubject data fusion". The data processed here consists of simultaneous MEG/EEG recordings from 16 participants performing a simple visual recognition task from presentations of famous, unfamiliar and scrambled faces.
The analysis is split in two tutorial pages: the present tutorial describes the detailed analysis of one single subject; the second tutorial describes batch processing and group analysis of all 16 participants.
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
License
This dataset was obtained from the OpenfMRI project (http://www.openfmri.org), accession #ds117. It is made available under the Creative Commons Attribution 4.0 International Public License. Please cite the following reference if you use these data:
Wakeman DG, Henson RN, A multi-subject, multi-modal human neuroimaging dataset, Scientific Data (2015)
Any questions regarding the data, please contact: rik.henson@mrc-cbu.cam.ac.uk
Presentation of the experiment
Experiment
16 subjects (the original version of the dataset included 19 subjects, but 3 were excluded from the group analysis for various reasons)
- 6 acquisition runs of approximately 10mins for each subject
- Presentation of series of images: familiar faces, unfamiliar faces, phase-scrambled faces
- Participants had to judge the left-right symmetry of each stimulus
- Total of nearly 300 trials for each of the 3 conditions
MEG acquisition
Acquisition at 1100Hz with an Elekta-Neuromag VectorView system (simultaneous MEG+EEG).
- Recorded channels (404):
- 102 magnetometers
- 204 planar gradiometers
- 70 EEG electrodes recorded with a nose reference.
MEG data have been "cleaned" using Signal-Space Separation as implemented in MaxFilter 2.2.
- A Polhemus device was used to digitize three fiducial points and a large number of other points across the scalp, which can be used to coregister the M/EEG data with the structural MRI image.
- Stimulation triggers: The triggers related with the visual presentation are saved in the STI101 channel, with the following event codes (bit 3 = face, bit 4 = unfamiliar, bit 5 = scrambled):
- Famous faces: 5 (00101), 6 (00110), 7 (00111)
- Unfamiliar faces: 13 (01101), 14 (01110), 15 (01111)
- Scrambled images: 17 (10001), 18 (10010), 19 (10011)
Delays between the trigger in STI101 and the actual presentation of stimulus: 34.5ms
The data distribution includes MEG noise recordings acquired around the dates of the experiment, processed with MaxFilter 2.2 in the same way as the experimental data.
Subject anatomy
- MRI data acquired on a 3T Siemens TIM Trio: 1x1x1mm T1-weighted structural MRI.
- The face was removed from the strucural images for anonymization purposes.
Processed with FreeSurfer 5.3.
Download and installation [TODO]
First, make sure you have enough space on your hard drive, at least 350Gb:
Raw files: 100Gb
Processed files: 250Gb
The data is hosted on the OpenfMRI website: https://openfmri.org/dataset/ds000117/
- Download the following files for revision 1.0.0 (approximately 70Gb):
FreeSurfer segmentation: NOT AVAILABLE YET [TODO]
Empty room recordings: NOT AVAILABLE YET [TODO]
Unzip all the .zip files in the same folder.
Reminder: Do not save the downloaded files in the Brainstorm folders (program or database folders)MaxFilter/tSSS: The MEG recordings from this study were recorded on an Elekta MEG, and need to be processed with cleaning algorithms SSS or tSSS. Brainstorm currently does not offer any free alternative to Elekta's MaxFilter, therefore in this tutorial we will import the recordings already processed with MaxFilter's tSSS, available in the "derivatives" folder of the BIDS architecture. To save disk space, will not download the raw MEG recordings. This means we will not get any of the additional files available in the BIDS structure (headshape, events, channels) but it doesn't matter because all the information is directly available in the .fif files.
Start Brainstorm (Matlab scripts or stand-alone version). For help, see the Installation page.
Select the menu File > Create new protocol. Name it "TutorialVisual" and select the options:
"No, use individual anatomy",
"No, use one channel file per condition".
Import the anatomy
This dataset is formatted following the BIDS-MEG specifications, therefore we could import all the relevant information (MRI, FreeSurfer segmentation, MEG+EEG recordings) in just one click, with the menu File > Load protocol > Import BIDS dataset, as described in the online tutorial MEG resting state & OMEGA database. However, because your own data might not be organized following the BIDS standards, in this tutorial we preferred illustrating all the detailed steps for importing the data rather than the BIDS shortcut.
This page explains how to import and process the first run of subject #01 only. All the other files will have to be processed in the same way.
- Switch to the "anatomy" view.
Right-click on the TutorialVisual folder > New subject > sub-01
- Leave the default options you defined for the protocol.
Right-click on the subject node > Import anatomy folder:
Set the file format: "FreeSurfer folder"
Select the folder: derivatives/freesurfer/sub-01
- Number of vertices of the cortex surface: 15000 (default value)
When asked to select the anatomical fiducials, click on "Compute MNI transformation". This will register the MRI volume to an MNI template with an affine transformation, using SPM functions embedded in Brainstorm (spm_maff8.m). This will also place default fiducial points NAS/LPA/RPA in the MRI, based on MNI coordinates. We will use head points for the MEG-MRI coregistration, therefore we don't need accurate anatomical positions here.
- Note that if you don't have a good digitized head shape for the subject, or if the final registration MEG-MRI doesn't look good with this head shape, you should not skip this step and mark accurate NAS/LPA/RPA fiducial points in the MRI, using the same anatomical convention that was used during the MRI acquisition.
- Then click on [Save] to start the import.
At the end of the process, make sure that the file "cortex_15000V" is selected (downsampled pial surface, that will be used for the source estimation). If it is not, double-click on it to select it as the default cortex surface. Do not worry about the big holes in the head surface, parts of MRI have been remove voluntarily for anonymization purposes.
All the anatomical atlases generated by FreeSurfer were automatically imported: the cortical atlases (Desikan-Killiany, Mindboggle, Destrieux, Brodmann) and the sub-cortical regions (ASEG atlas).
Access the recordings
Link the recordings
We need to attach the continuous .fif files containing the recordings to the database.
- Switch to the "functional data" view.
Right-click on the subject folder > Review raw file.
Select the file format: "MEG/EEG: Neuromag FIFF (*.fif)"
Go to the folder: derivatives/meg-derivatives/sub-01/ses-meg/meg
Select file: sub-01_ses-meg_task-facerecognition_run-01_proc-tsss_meg.fif
Events: Ignore. We will read the stimulus triggers later.
Refine registration now? NO
The head points that are available in the FIF files contain all the points that were digitized during the MEG acquisition, including the ones corresponding to the parts of the face that have been removed from the MRI. If we run the fitting algorithm, all the points around the nose will not match any close points on the head surface, leading to a wrong result. We will first remove the face points and then run the registration manually.
Channel classification
A few non-EEG channels are mixed in with the EEG channels, we need to change this before applying any operation on the EEG channels.
Right-click on the channel file > Edit channel file. Double-click on a cell to edit it.
Change the type of EEG062 to EOG (electrooculogram).
Change the type of EEG063 to ECG (electrocardiogram).
Change the type of EEG061 and EEG064 to NOSIG (or any other non-informative type). Close the window and save the modifications.
MRI registration
At this point, the registration MEG/MRI is based only on the three anatomical landmarks NAS/LPA/RPA, which are not even accurately set (we used default MNI positions). All the MRI scans were anonymized (defaced) so all the head points below the nasion cannot be used. We will try to refine this registration using the additional head points that were digitized above the nasion.
Right-click on the channel file > Digitized head points > Remove points below nasion.
(45 points removed, 92 head points left)
Right-click on the channel file > MRI registration > Refine using head points.
MEG/MRI registration, before (left) and after (right) this automatic registration procedure:
Right-click on the channel file > MRI registration > EEG: Edit.
Click on [Project electrodes on surface], then close the figure to save the modifications.
Read stimulus triggers
We need to read the stimulus markers from the STI channels. The following tasks can be done in an interactive way with menus in the Record tab, as in the introduction tutorials. We will illustrate here how to do this with the pipeline editor, it will be easier to batch it for all the runs and all the subjects.
- In Process1, select the "Link to raw file", click on [Run].
Select process Events > Read from channel, Channel: STI101, Detection mode: Bit.
Do not execute the process, we will add other processes to classify the markers.
- We want to create three categories of events, based on their numerical codes:
Famous faces: 5 (00101), 6 (00110), 7 (00111) => Bit 3 only
Unfamiliar faces: 13 (01101), 14 (01110), 15 (01111) => Bit 3 and 4
Scrambled images: 17 (10001), 18 (10010), 19 (10011) => Bit 5 only
- We will start by creating the category "Unfamiliar" (combination of events "3" and "4") and remove the remove the initial events. Then we just have to rename the remaining "3" in "Famous", and all the "5" in "Scrambled".
Add process Events > Group by name: "Unfamiliar=3,4", Delay=0, Delete original events
Add process Events > Rename event: 3 => Famous
Add process Events > Rename event: 5 => Scrambled
Add process Events > Add time offset to correct for the presentation delays:
Event names: "Famous, Unfamiliar, Scrambled", Time offset = 34.5ms
Finally run the script. Double-click on the recordings to make sure the labels were detected correctly. You can delete the unwanted events that are left in the recordings (1,2,9,13). Edit the colors if you don't like the default ones.
Pre-processing
Spectral evaluation
- Keep the "Link to raw file" in Process1.
Run process Frequency > Power spectrum density (Welch) with the options illustrated below.
Right-click on the PSD file > Power spectrum.
The MEG sensors look awful, because of one small segment of data located around 248s. Open the MEG recordings and scroll to 248s (just before the first Unfamiliar event).
In these recordings, the continuous head tracking was activated, but it starts only at the time of the stimulation (248s) while the acquisition of the data starts 20s before (226s). The first 20s do not have head localization coils (HPI) coils activity and are not corrected by MaxFilter. After 248s, the HPI coils are on, and MaxFilter filters them out. The transition between the two states is not smooth and creates important distortions in the spectral domain. For a proper evaluation of the recordings, we should compute the PSD only after the HPI coils are turned on.
Run process Events > Detect cHPI activity (Elekta). This detects the changes in the cHPI activity from channel STI201 and marks all the data without head localization as bad.
Re-run the process Frequency > Power spectrum density (Welch). All the bad segments are excluded from the computation, therefore the PSD is now estimated only with the data after 248s.
- Observations:
- Three groups of sensors, from top to bottom: EEG, MEG gradiometers, MEG magnetometers.
Power lines: 50 Hz and harmonics
- Alpha peak around 10 Hz
Artifacts due to Elekta electronics (HPI coils): 293Hz, 307Hz, 314Hz, 321Hz, 328Hz.
Peak from unknown source at 103.4Hz in the MEG only.
Suspected bad EEG channels: EEG016
- Close all the windows.
Remove line noise
- Keep the "Link to raw file" in Process1.
Select process Pre-process > Notch filter to remove the line noise (50-200Hz).
Add the process Frequency > Power spectrum density (Welch).
Double-click on the PSD for the new continuous file to evaluate the quality of the correction.
- Close all the windows (use the [X] button at the top-right corner of the Brainstorm window).
EEG reference and bad channels
Right-click on link to the processed file ("Raw | notch(50Hz ...") > EEG > Display time series.
Select channel EEG016 and mark it as bad (using the popup menu or pressing the Delete key).
In the Record tab, menu Artifacts > Re-reference EEG > "AVERAGE".
At the end, the window "select active projectors" is open to show the new re-referencing projector. Just close this window. To get it back, use the menu Artifacts > Select active projectors.
Artifact detection
Heartbeats: Detection
Empty the Process1 list (right-click > Clear list).
- Drag and drop the continuous processed file ("Raw | notch(50Hz...)") to the Process1 list.
Run process Events > Detect heartbeats: Channel name=EEG063, All file, Event name=cardiac
Eye blinks: Detection
- In many of the other tutorials, we detect the blinks and remove them with SSP. In this experiment, we are particularly interested in the subject's response to seeing the stimulus. Therefore we will exclude from the analysis all the recordings contaminated with blinks or other eye movements.
Run process Events > Detect events above threshold:
Event name=blink_BAD, Channel=EEG062, All file, Maximum threshold=100, Threshold units=uV, Filter=[0.30,20.00]Hz, Use absolute value of signal.
- In many other tutorials, we use the process "Detect eye blinks", which creates simple events indicating the peak of the EOG signal during a blink. In this study, we preferred using this process "Detect events above threshold" because it creates extended events marking as bad all the segment for which the EOG value is above a given threshold. It's more manual but it makes sure we really exclude all the ocular activity.
Inspect visually the two new cateogries of events: cardiac and blink_BAD.
- Close all the windows (using the [X] button).
Heartbeats: Correction with SSP
- Keep the "Raw | notch" file selected in Process1.
Select process Artifacts > SSP: Heartbeats > Sensor type: MEG MAG
Add process Artifacts > SSP: Heartbeats > Sensor type: MEG GRAD - Run the execution
Double-click on the continuous file to show all the MEG sensors.
In the Record tab, select sensors "Left-temporal".Menu Artifacts > Select active projectors.
In category cardiac/MEG MAG: Select component #1 and view topography.
In category cardiac/MEG GRAD: Select component #1 and view topography.
Make sure that selecting the two components removes the cardiac artifact. Then click [Save].
Additional bad segments
Process1>Run: Select the process "Events > Detect other artifacts". This should be done separately for MEG and EEG to avoid confusion about which sensors are involved in the artifact.
Display the MEG sensors. Review the segments that were tagged as artifact, determine if each event represents an artifact and then mark the time of the artifact as BAD. This can be done by selecting the time window around the artifact, then right-click > Reject time segment. Note that this detection process marks 1-second segments but the artifact can be shorter.
- Once all the events in the two categories are reviewed and bad segments are marked, the two categories (1-7Hz and 40-240Hz) can be deleted.
- Do this detection and review again for the EEG.
SQUID jumps
MEG signals recorded with Elekta-Neuromag systems frequently contain SQUID jumps (more information). These sharp steps followed by a change of baseline value are easy to identify visually but more complicated to detect automatically.
The process "Detect other artifacts" usually detects most of them in the category "1-7Hz". If you observe that some are skipped, you can try re-running it with a higher sensitivity. It is important to review all the sensors and all the time in each run to be sure these events are marked as bad segments.
Epoching and averaging
Import epochs
- Keep the "Raw | notch" file selected in Process1.
Select process: Import recordings > Import MEG/EEG: Events (do not run immediately)
Event names "Famous, Unfamiliar, Scrambled", All file, Epoch time=[-500,1200]msAdd process: Pre-process > Remove DC offset: Baseline=[-500,-0.9]ms - Run execution.
Average by run
- In Process1, select all the imported trials.
Run process: Average > Average files: By trial groups (folder average)
Review EEG ERP
EEG evoked response (famous, scrambled, unfamiliar), low-pass filter at 32Hz:
Open the Cluster tab and create a cluster with the channel EEG065 (button [NEW IND]).
Select the cluster, select the three average files, right-click > Clusters time series (Overlay:Files).
- Basic observations for EEG065 (right parieto-occipital electrode):
- Around 170ms (N170): greater negative deflection for Famous than Scrambled faces.
- After 250ms: difference between Famous and Unfamiliar faces.
Source estimation [TODO]
MEG noise covariance: Empty room recordings [TODO]
The minimum norm model we will use next to estimate the source activity can be improved by modeling the the noise contaminating the data. The introduction tutorials explain how to estimate the noise covariance in different ways for EEG and MEG. For the MEG recordings we will use the empty room measurements we have, and for the EEG we will compute it from the pre-stimulus baselines we have in all the imported epochs.
Create a new subject: emptyroom
Right-click on the new subject > Review raw file.
Select file: sub-emptyroom/ses-meg/meg/090707_raw_st.fif (processed with MaxFilter/tSSS)
Do not apply default transformation, Ignore event channel.
Select this new file in Process1 and run process Pre-process > Notch filter: 50 100 150 200Hz. When using empty room measurements to compute the noise covariance, they must be processed exactly in the same way as the other recordings (same MaxFilter parameters, same frequency filters).
Right-click on the filtered noise recordings > Noise covariance > Compute from recordings:
Right-click on the Noise covariance > Copy to other subjects
EEG noise covariance: Pre-stimulus baseline
In folder sub-01/run-01_tsss_notch, select all the imported the imported trials, right-click > Noise covariance > Compute from recordings, Time=[-500,-0.9]ms, EEG only, Merge.
This computes the noise covariance only for EEG, and combines it with the existing MEG information.
BEM layers
We will compute a BEM forward model to estimate the brain sources from the EEG recordings. For this, we need some layers defining the separation between the different tissues of the head (scalp, inner skull, outer skull).
- Go to the anatomy view (first button above the database explorer).
Right-click on the subject folder > Generate BEM surfaces: The number of vertices to use for each layer depends on your computing power and the accuracy you expect. You can try for instance with 1082 vertices (scalp) and 642 vertices (outer skull and inner skull).
Forward model: EEG and MEG
- Go back to the functional view (second button above the database explorer).
Model used: Overlapping spheres for MEG, OpenMEEG BEM for EEG (more information).
In folder sub-01/run_01_tsss_notch, right-click on the channel file > Compute head model.
Keep all the default options . Expect this to take a while...
Inverse model: Minimum norm estimates
Right-click on the new head model > Compute sources [2016]: MEG MAG + GRAD (default options)
Right-click on the new head model > Compute sources [2016]: EEG (default bad channels).
At the end we have two inverse operators, that are shared for all the files of the run (single trials and averages). If we wanted to look at the run-level source averages, we could normalize the source maps with a Z-score wrt baseline. In this tutorial, we will first average across runs and normalize the subject-level averages. This will be done in the next tutorial (group analysis).
Time-frequency analysis
We will compute the time-frequency decomposition of each trial using Morlet wavelets, and average the power of the Morlet coefficients for each condition and each run separately. We will restrict the computation to the MEG magnetometers and the EEG channels to limit the computation time and disk usage.
- In Process1, select the imported trials "Famous" for run#01.
Run process Frequency > Time-frequency (Morlet wavelets): Sensor types=MEG MAG,EEG
Not normalized, Frequency=Log(6:20:60), Measure=Power, Save average
Double-click on the file to display it. In the Display tab, select the option "Hide edge effects" to exclude from the display all the values that could not be estimated in a reliable way. Let's extract only the good values from this file (-200ms to +900ms).
- In Process1, select the time-frequency file.
Run process Extract > Extract time: Time window=[-200, 900]ms, Overwrite input files
Display the file again, observe that all the possibly bad values are gone.
You can display all the sensors at once (MEG MAG or EEG): right-click > 2D Layout (maps).
- Repeat these steps for other conditions (Scrambled and Unfamiliar) and the other runs (2-6). There is no way with this process to compute all the averages at once, as it is with the process "Average files". This will be easier to run from a script.
- If we wanted to look at the run-level source averages, we could normalize these time-frequency maps. In this tutorial, we will first average across runs and normalize the subject-level averages. This will be done in the next tutorial (group analysis).
Scripting
We have now all the files we need for the group analysis (next tutorial). We need to repeat the same operations for all the runs and all the subjects. Some of these steps are fully automatic and take a lot of time (filtering, computing the forward model), they should be executed from a script.
However, we recommend you always review manually some of the pre-processing steps (selection of the bad segments and bad channels, SSP/ICA components). Do not trust blindly any fully automated cleaning procedure.
For the strict reproducibility of this analysis, we provide a script that processes all the 19 subjects: brainstorm3/toolbox/script/tutorial_visual_single.m (execution time: 10-30 hours)
Report for the first subject: report_TutorialVisual_sub001.html
1 function tutorial_visual_single(bids_dir, reports_dir)
2 % TUTORIAL_VISUAL_SINGLE: Runs the Brainstorm/SPM group analysis pipeline (single subject, BIDS version).
3 %
4 % ONLINE TUTORIALS: https://neuroimage.usc.edu/brainstorm/Tutorials/VisualSingle
5 %
6 % INPUTS:
7 % - bids_dir: Path to folder ds000117 (https://openneuro.org/datasets/ds000117)
8 % |- derivatives/freesurfer/sub-XX : Segmentation folders generated with FreeSurfer
9 % |- derivatives/meg_derivatives/sub-XX/ses-meg/meg/*.fif : MEG+EEG recordings (processed with MaxFilter's SSS)
10 % |- derivatives/meg_derivatives/sub-emptyroom/ses-meg/meg/*.fif : Empty room measurements
11 % - reports_dir: If defined, exports all the reports as HTML to this folder
12
13 % @=============================================================================
14 % This function is part of the Brainstorm software:
15 % https://neuroimage.usc.edu/brainstorm
16 %
17 % Copyright (c) University of Southern California & McGill University
18 % This software is distributed under the terms of the GNU General Public License
19 % as published by the Free Software Foundation. Further details on the GPLv3
20 % license can be found at http://www.gnu.org/copyleft/gpl.html.
21 %
22 % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
23 % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
24 % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
25 % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
26 % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
27 %
28 % For more information type "brainstorm license" at command prompt.
29 % =============================================================================@
30 %
31 % Author: Francois Tadel, Elizabeth Bock, 2016-2018
32
33
34 %% ===== SCRIPT VARIABLES =====
35 % Full list of subjects to process
36 SubjectNames = {'sub-01', 'sub-02', 'sub-03', 'sub-04', 'sub-05', 'sub-06', 'sub-07', 'sub-08', ...
37 'sub-09', 'sub-10', 'sub-11', 'sub-12', 'sub-13', 'sub-14', 'sub-15', 'sub-16'};
38 % Empty-room dates for each subject (so that we can match automatically recordings with empty-room)
39 EmptyRoomSubj = 'sub-emptyroom';
40 AcquisitionDates = {'09-Apr-2009', '06-May-2009', '11-May-2009', '18-May-2009', '15-May-2009', '15-May-2009', '15-May-2009', '15-May-2009', ...
41 '15-May-2009', '15-May-2009', '01-Jun-2009', '01-Jun-2009', '01-Jun-2009', '26-Nov-2009', '08-Dec-2009', '08-Dec-2009'};
42 % Bad channels {iSubj} = {Run01, Run02, Run03, Run04, Run05, Run06}
43 BadChannels{1} = {'EEG016', 'EEG070', 'EEG050',{'EEG008','EEG050'}, [], []};
44 BadChannels{2} = {{'EEG027', 'EEG030', 'EEG038'}, 'EEG010', 'EEG010', 'EEG010', 'EEG010', 'EEG010'};
45 BadChannels{3} = {{'EEG008','EEG017'}, {'EEG008','EEG017'}, {'EEG008','EEG017'}, {'EEG008','EEG017'}, {'EEG008','EEG017','EEG001'}, {'EEG008','EEG017','EEG020'}};
46 BadChannels{4} = {{'EEG038'}, {'EEG038','EEG001','EEG016'}, {'EEG038','EEG001','EEG016'}, {'EEG038','EEG001'}, {'EEG038','EEG001','EEG016'}, {'EEG038','EEG001','EEG016'}};
47 BadChannels{5} = {'EEG001', 'EEG001', [], [], [], []};
48 BadChannels{6} = {'EEG068', [], 'EEG004', [], [], []};
49 BadChannels{7} = {[], [], {'EEG004','EEG008'}, {'EEG004','EEG008'},{'EEG004','EEG008','EEG043','EEG045','EEG047'}, {'EEG004','EEG008'}};
50 BadChannels{8} = {[], [], [], [], [], []};
51 BadChannels{9} = {[], 'EEG004', 'EEG004', [], 'EEG004', 'EEG004'};
52 BadChannels{10} = {[], [], [], [], [], []};
53 BadChannels{11} = {{'EEG010','EEG050'}, 'EEG050', 'EEG050', 'EEG050', 'EEG050', 'EEG050'};
54 BadChannels{12} = {{'EEG024','EEG057'}, {'EEG024','EEG057'}, {'EEG024','EEG057'}, {'EEG024','EEG057','EEG070'}, {'EEG024','EEG057'}, {'EEG024','EEG057','EEG070'}};
55 BadChannels{13} = {'EEG009', 'EEG009', {'EEG009','EEG057','EEG69'}, 'EEG009', {'EEG009','EEG044'}, {'EEG009','EEG044'}};
56 BadChannels{14} = {'EEG029', 'EEG029', 'EEG029', {'EEG004','EEG008','EEG016','EEG029'}, {'EEG004','EEG008','EEG016','EEG029'}, {'EEG004','EEG008','EEG016','EEG029'}};
57 BadChannels{15} = {'EEG038', 'EEG038', 'EEG038', 'EEG038', {'EEG054','EEG038'}, 'EEG038'};
58 BadChannels{16} = {'EEG008', 'EEG008', 'EEG008', 'EEG008', 'EEG008', 'EEG008'};
59 % SSP components to remove {iSubj} = {sRun01, sRun02, sRun03, sRun03, sRun04, sRun05, sRun06}, sRun0X={ECG_GRAD,ECG_MAG}
60 SspSelect{1} = {{1,1}, {1,1}, {1,1}, {1,1}, {1,1}, {1,1}};
61 SspSelect{2} = {{1,1}, {1,1}, {1,1}, {1,1}, {3,1}, {1,1}};
62 SspSelect{3} = {{[],1}, {1,1}, {1,1}, {1,1}, {1,1}, {1,1}};
63 SspSelect{4} = {{[],1}, {[],1}, {[],1}, {[],1}, {[],1}, {[],1}};
64 SspSelect{5} = {{2,1}, {1,1}, {1,1}, {[],1}, {1,1}, {1,1}};
65 SspSelect{6} = {{2,1}, {2,1}, {1,1}, {2,1}, {1,1}, {2,1}};
66 SspSelect{7} = {{1,1}, {1,1}, {1,1}, {1,1}, {1,1}, {1,1}};
67 SspSelect{8} = {{1,1}, {1,1}, {[],1}, {2,1}, {1,1}, {2,1}};
68 SspSelect{9} = {{1,1}, {1,1}, {[],1}, {[],1}, {1,1}, {1,1}};
69 SspSelect{10} = {{1,1}, {1,1}, {1,1}, {1,1}, {1,1}, {1,1}};
70 SspSelect{11} = {{[],1}, {[],1}, {[],1}, {[],1}, {[],[]}, {[],[]}};
71 SspSelect{12} = {{[1,2],[1,2]}, {1,1}, {1,1}, {1,1}, {1,1}, {1,1}};
72 SspSelect{13} = {{[],[]}, {[],[]}, {[],[]}, {[],[]}, {[],[]}, {[],[]}};
73 SspSelect{14} = {{1,1}, {1,1}, {1,1}, {1,1}, {1,1}, {1,1}};
74 SspSelect{15} = {{1,1}, {1,1}, {1,1}, {1,1}, {1,1}, {1,1}};
75 SspSelect{16} = {{1,1}, {1,1}, {1,1}, {2,1}, {1,1}, {1,1}};
76
77
78 %% ===== CREATE PROTOCOL =====
79 % Start brainstorm without the GUI
80 if ~brainstorm('status')
81 brainstorm nogui
82 end
83 % Output folder for reports
84 if (nargin < 2) || isempty(reports_dir) || ~isdir(reports_dir)
85 reports_dir = [];
86 end
87 % You have to specify the folder in which the tutorial dataset is unzipped
88 if (nargin < 1) || isempty(bids_dir) || ~file_exist(bids_dir) || ~file_exist(bst_fullfile(bids_dir, 'derivatives')) || ~file_exist(bst_fullfile(bids_dir, 'dataset_description.json'))
89 error('The first argument must be the full path to the tutorial folder.');
90 end
91 % The protocol name has to be a valid folder name (no spaces, no weird characters...)
92 ProtocolName = 'TutorialVisual';
93 % Delete existing protocol
94 gui_brainstorm('DeleteProtocol', ProtocolName);
95 % Create new protocol
96 gui_brainstorm('CreateProtocol', ProtocolName, 0, 0);
97 % Set visualization filters: 40Hz low-pass, no high-pass
98 panel_filter('SetFilters', 1, 40, 0, [], 0, [], 0, 0);
99 % Set colormap: local color scale
100 bst_colormaps('SetMaxMode', 'meg', 'local');
101 bst_colormaps('SetMaxMode', 'eeg', 'local');
102
103
104 %% ===== PRE-PROCESS AND IMPORT =====
105 for iSubj = 1:16
106 % Start a new report (one report per subject)
107 bst_report('Start');
108 disp(sprintf('\n===== IMPORT: SUBJECT #%d =====\n', iSubj));
109
110 % If subject already exists: delete it
111 [sSubject, iSubject] = bst_get('Subject', SubjectNames{iSubj});
112 if ~isempty(sSubject)
113 db_delete_subjects(iSubject);
114 end
115
116 % ===== FILES TO IMPORT =====
117 % Build the path of the files to import
118 AnatDir = fullfile(bids_dir, 'derivatives', 'freesurfer', SubjectNames{iSubj}, 'ses-mri', 'anat');
119 DataDir = fullfile(bids_dir, 'derivatives', 'meg_derivatives', SubjectNames{iSubj}, 'ses-meg', 'meg');
120 % Check if the folder contains the required files
121 if ~file_exist(AnatDir)
122 error(['The folder "' AnatDir '" does not exist.']);
123 end
124 if ~file_exist(DataDir)
125 error(['The folder "' DataDir '" does not exist.']);
126 end
127
128 % ===== ANATOMY =====
129 % Process: Import anatomy folder
130 bst_process('CallProcess', 'process_import_anatomy', [], [], ...
131 'subjectname', SubjectNames{iSubj}, ...
132 'mrifile', {AnatDir, 'FreeSurfer'}, ...
133 'nvertices', 15000);
134
135 % ===== PROCESS EACH RUN =====
136 for iRun = 1:6
137 % Files to import
138 FifFile = bst_fullfile(DataDir, sprintf('%s_ses-meg_task-facerecognition_run-%02d_proc-sss_meg.fif', SubjectNames{iSubj}, iRun));
139
140 % ===== LINK CONTINUOUS FILE =====
141 % Process: Create link to raw file
142 sFileRaw = bst_process('CallProcess', 'process_import_data_raw', [], [], ...
143 'subjectname', SubjectNames{iSubj}, ...
144 'datafile', {FifFile, 'FIF'}, ...
145 'channelreplace', 1, ...
146 'channelalign', 0);
147 % Set acquisition date
148 panel_record('SetAcquisitionDate', sFileRaw.iStudy, AcquisitionDates{iSubj});
149
150 % ===== PREPARE CHANNEL FILE =====
151 % Process: Set channels type
152 bst_process('CallProcess', 'process_channel_settype', sFileRaw, [], ...
153 'sensortypes', 'EEG061, EEG064', ...
154 'newtype', 'NOSIG');
155 bst_process('CallProcess', 'process_channel_settype', sFileRaw, [], ...
156 'sensortypes', 'EEG062', ...
157 'newtype', 'EOG');
158 bst_process('CallProcess', 'process_channel_settype', sFileRaw, [], ...
159 'sensortypes', 'EEG063', ...
160 'newtype', 'ECG');
161
162 % Process: Remove head points
163 sFileRaw = bst_process('CallProcess', 'process_headpoints_remove', sFileRaw, [], ...
164 'zlimit', 0);
165 % Process: Refine registration
166 sFileRaw = bst_process('CallProcess', 'process_headpoints_refine', sFileRaw, []);
167 % Process: Project electrodes on scalp
168 sFileRaw = bst_process('CallProcess', 'process_channel_project', sFileRaw, []);
169
170 % Process: Snapshot: Sensors/MRI registration
171 bst_process('CallProcess', 'process_snapshot', sFileRaw, [], ...
172 'target', 1, ... % Sensors/MRI registration
173 'modality', 1, ... % MEG (All)
174 'orient', 1, ... % left
175 'Comment', sprintf('MEG/MRI Registration: Subject #%d, Run #%d', iSubj, iRun));
176 bst_process('CallProcess', 'process_snapshot', sFileRaw, [], ...
177 'target', 1, ... % Sensors/MRI registration
178 'modality', 4, ... % EEG
179 'orient', 1, ... % left
180 'Comment', sprintf('EEG/MRI Registration: Subject #%d, Run #%d', iSubj, iRun));
181
182 % ===== IMPORT TRIGGERS =====
183 % Process: Read from channel
184 bst_process('CallProcess', 'process_evt_read', sFileRaw, [], ...
185 'stimchan', 'STI101', ...
186 'trackmode', 2, ... % Bit: detect the changes for each bit independently
187 'zero', 0);
188 % Process: Group by name
189 bst_process('CallProcess', 'process_evt_groupname', sFileRaw, [], ...
190 'combine', 'Unfamiliar=3,4', ...
191 'dt', 0, ...
192 'delete', 1);
193 % Process: Rename event
194 bst_process('CallProcess', 'process_evt_rename', sFileRaw, [], ...
195 'src', '3', ...
196 'dest', 'Famous');
197 % Process: Rename event
198 bst_process('CallProcess', 'process_evt_rename', sFileRaw, [], ...
199 'src', '5', ...
200 'dest', 'Scrambled');
201 % Process: Add time offset
202 bst_process('CallProcess', 'process_evt_timeoffset', sFileRaw, [], ...
203 'info', [], ...
204 'eventname', 'Famous, Unfamiliar, Scrambled', ...
205 'offset', 0.0345);
206 % Process: Delete events
207 bst_process('CallProcess', 'process_evt_delete', sFileRaw, [], ...
208 'eventname', '1,2,6,7,8,9,10,11,12,13,14,15,16');
209 % Process: Detect cHPI activity (Elekta):STI201
210 bst_process('CallProcess', 'process_evt_detect_chpi', sFileRaw, [], ...
211 'eventname', 'chpi_bad', ...
212 'channelname', 'STI201', ...
213 'method', 'off'); % Mark as bad when the HPI coils are OFF
214
215 % ===== FREQUENCY FILTERS =====
216 % Process: Notch filter: 50Hz 100Hz 150Hz 200Hz
217 sFileClean = bst_process('CallProcess', 'process_notch', sFileRaw, [], ...
218 'freqlist', [50, 100, 150, 200], ...
219 'sensortypes', 'MEG, EEG', ...
220 'read_all', 0);
221 % Process: Power spectrum density (Welch)
222 sFilesPsd = bst_process('CallProcess', 'process_psd', [sFileRaw, sFileClean], [], ...
223 'timewindow', [], ...
224 'win_length', 4, ...
225 'win_overlap', 50, ...
226 'sensortypes', 'MEG, EEG', ...
227 'edit', struct(...
228 'Comment', 'Power', ...
229 'TimeBands', [], ...
230 'Freqs', [], ...
231 'ClusterFuncTime', 'none', ...
232 'Measure', 'power', ...
233 'Output', 'all', ...
234 'SaveKernel', 0));
235 % Process: Snapshot: Frequency spectrum
236 bst_process('CallProcess', 'process_snapshot', sFilesPsd, [], ...
237 'target', 10, ... % Frequency spectrum
238 'Comment', sprintf('Power spctrum: Subject #%d, Run #%d', iSubj, iRun));
239
240 % ===== BAD CHANNELS =====
241 if ~isempty(BadChannels{iSubj}{iRun})
242 % Process: Set bad channels
243 bst_process('CallProcess', 'process_channel_setbad', sFileClean, [], ...
244 'sensortypes', BadChannels{iSubj}{iRun});
245 end
246
247 % ===== EEG REFERENCE =====
248 % Process: Re-reference EEG
249 bst_process('CallProcess', 'process_eegref', sFileClean, [], ...
250 'eegref', 'AVERAGE', ...
251 'sensortypes', 'EEG');
252
253 % ===== DETECT ARTIFACTS ======
254 % Process: Detect heartbeats
255 bst_process('CallProcess', 'process_evt_detect_ecg', sFileClean, [], ...
256 'channelname', 'EEG063', ...
257 'timewindow', [], ...
258 'eventname', 'cardiac');
259 % Different amplitude thresholds for different subjects
260 if strcmpi(SubjectNames{iSubj}, 'sub-05')
261 thresholdMAX = 50;
262 else
263 thresholdMAX = 100;
264 end
265 % Process: Detect: blink_BAD - Detects all events where the amplitude exceeds 100uV
266 bst_process('CallProcess', 'process_evt_detect_threshold', sFileClean, [], ...
267 'eventname', 'blink_BAD', ...
268 'channelname', 'EEG062', ...
269 'timewindow', [], ...
270 'thresholdMAX', thresholdMAX, ...
271 'units', 3, ... % uV (10^-6)
272 'bandpass', [0.3, 20], ...
273 'isAbsolute', 1, ...
274 'isDCremove', 0);
275
276 % ===== SSP COMPUTATION =====
277 % Process: SSP ECG: cardiac
278 bst_process('CallProcess', 'process_ssp_ecg', sFileClean, [], ...
279 'eventname', 'cardiac', ...
280 'sensortypes', 'MEG GRAD', ...
281 'usessp', 1, ...
282 'select', SspSelect{iSubj}{iRun}{1});
283 bst_process('CallProcess', 'process_ssp_ecg', sFileClean, [], ...
284 'eventname', 'cardiac', ...
285 'sensortypes', 'MEG MAG', ...
286 'usessp', 1, ...
287 'select', SspSelect{iSubj}{iRun}{2});
288 % Process: Snapshot: SSP projectors
289 bst_process('CallProcess', 'process_snapshot', sFileClean, [], ...
290 'target', 2, ...
291 'Comment', sprintf('Subject #%d, Run #%d', iSubj, iRun)); % SSP projectors
292
293 % ===== IMPORT BAD EVENTS =====
294 % Get bad segments: this is typically done manually, not from a script
295 BadSegments = GetBadSegments(iSubj, iRun);
296 % Process: Import from file
297 bst_process('CallProcess', 'process_evt_import', sFileClean, [], ...
298 'evtfile', {BadSegments, 'ARRAY-TIMES'}, ...
299 'evtname', 'BAD');
300
301 % ===== IMPORT TRIALS =====
302 % Process: Import MEG/EEG: Events
303 sFilesEpochs = bst_process('CallProcess', 'process_import_data_event', sFileClean, [], ...
304 'subjectname', SubjectNames{iSubj}, ...
305 'condition', '', ...
306 'eventname', 'Famous, Scrambled, Unfamiliar', ...
307 'timewindow', [], ...
308 'epochtime', [-0.5, 1.2], ...
309 'createcond', 0, ...
310 'ignoreshort', 1, ...
311 'usectfcomp', 1, ...
312 'usessp', 1, ...
313 'freq', [], ...
314 'baseline', [-0.5, -0.0009]);
315
316 % ===== AVERAGE: RUN =====
317 % Process: Average: By trial group (folder average)
318 sFilesAvg = bst_process('CallProcess', 'process_average', sFilesEpochs, [], ...
319 'avgtype', 5, ... % By trial group (folder average)
320 'avg_func', 1, ... % Arithmetic average: mean(x)
321 'weighted', 0, ...
322 'keepevents', 0);
323 % Process: Snapshot: Recordings time series
324 bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ...
325 'target', 5, ... % Recordings time series
326 'modality', 4, ... % EEG
327 'time', 0.11, ...
328 'Comment', sprintf('Subject #%d, Run #%d', iSubj, iRun));
329 % Process: Snapshot: Recordings topography
330 bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ...
331 'target', 6, ... % Recordings topography (one time)
332 'modality', 4, ... % EEG
333 'time', 0.11, ...
334 'Comment', sprintf('Subject #%d, Run #%d', iSubj, iRun));
335
336 % ===== COMPUTE NOISECOV: EEG =====
337 % Process: Compute covariance (noise or data)
338 bst_process('CallProcess', 'process_noisecov', sFilesEpochs, [], ...
339 'baseline', [-0.5, -0.0009], ...
340 'sensortypes', 'EEG', ...
341 'target', 1, ... % Noise covariance (covariance over baseline time window)
342 'dcoffset', 1, ... % Block by block, to avoid effects of slow shifts in data
343 'identity', 0, ...
344 'copycond', 0, ...
345 'copysubj', 0, ...
346 'replacefile', 1); % Replace
347 end
348
349 % Save report
350 ReportFile = bst_report('Save', []);
351 if ~isempty(reports_dir) && ~isempty(ReportFile)
352 bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_' SubjectNames{iSubj} '.html']));
353 end
354 end
355
356
357 %% ===== EMPTY ROOM RECORDINGS =====
358 disp(sprintf('\n===== IMPORT: EMPTY-ROOM =====\n'));
359 % Loop on all the noise sessions
360 NoiseFiles = {};
361 for ses = {'20090409', '20090506', '20090511', '20090515', '20090518', '20090601', '20091126', '20091208'}
362 NoiseFiles{end+1} = fullfile(bids_dir, 'derivatives', 'meg_derivatives', EmptyRoomSubj, ['ses-' ses{1}], 'meg', ['sub-emptyroom_ses-' ses{1} '_task-noise_proc-sss_meg.fif']);
363 end
364 % Process: Create link to raw file
365 sFilesNoise = bst_process('CallProcess', 'process_import_data_raw', [], [], ...
366 'subjectname', EmptyRoomSubj, ...
367 'datafile', {NoiseFiles, 'FIF'}, ...
368 'channelreplace', 1, ...
369 'channelalign', 0);
370 % Process: Notch filter: 50Hz 100Hz 150Hz 200Hz
371 sFileNoiseClean = bst_process('CallProcess', 'process_notch', sFilesNoise, [], ...
372 'freqlist', [50, 100, 150, 200], ...
373 'sensortypes', 'MEG, EEG', ...
374 'read_all', 0);
375 % Process: Compute noise covariance
376 bst_process('CallProcess', 'process_noisecov', sFileNoiseClean, [], ...
377 'baseline', [], ...
378 'sensortypes', 'MEG', ...
379 'target', 1, ... % Noise covariance (covariance over baseline time window)
380 'dcoffset', 1, ... % Block by block, to avoid effects of slow shifts in data
381 'identity', 0, ...
382 'copycond', 1, ...
383 'copysubj', 1, ...
384 'copymatch', 1, ...
385 'replacefile', 2); % Merge
386
387
388 %% ===== SOURCE ESTIMATION =====
389 % Start a new report (one report for the source estimation of all the subjects)
390 bst_report('Start');
391 % Loop on the subjects: This loop is separated from the previous one, because we should
392 % compute the BEM surfaces after importing all the runs, so that the registration is done
393 % using the high resolution head surface, instead of the smooth scalp BEM layer.
394 for iSubj = 1:length(SubjectNames)
395 disp(sprintf('\n===== SOURCES: SUBJECT #%d =====\n', iSubj));
396
397 % ===== BEM SURFACES =====
398 % Process: Generate BEM surfaces
399 bst_process('CallProcess', 'process_generate_bem', [], [], ...
400 'subjectname', SubjectNames{iSubj}, ...
401 'nscalp', 1082, ...
402 'nouter', 642, ...
403 'ninner', 642, ...
404 'thickness', 4, ...
405 'method', 'brainstorm');
406
407 % ===== SELECT ALL AVERAGES =====
408 % Process: Select data files in: */*
409 sFilesAvg = bst_process('CallProcess', 'process_select_files_data', [], [], ...
410 'subjectname', SubjectNames{iSubj});
411 % Process: Select file comments with tag: Avg
412 sFilesAvg = bst_process('CallProcess', 'process_select_tag', sFilesAvg, [], ...
413 'tag', 'Avg'); % Select only the files with the tag
414
415 % ===== COMPUTE HEAD MODELS =====
416 % Process: Compute head model (only for the first run of the subject)
417 bst_process('CallProcess', 'process_headmodel', sFilesAvg(1), [], ...
418 'sourcespace', 1, ... % Cortex surface
419 'meg', 3, ... % Overlapping spheres
420 'eeg', 3, ... % OpenMEEG BEM
421 'ecog', 1, ... %
422 'seeg', 1, ... %
423 'openmeeg', struct(...
424 'BemSelect', [1, 1, 1], ...
425 'BemCond', [1, 0.0125, 1], ...
426 'BemNames', {{'Scalp', 'Skull', 'Brain'}}, ...
427 'BemFiles', {{}}, ...
428 'isAdjoint', 0, ...
429 'isAdaptative', 1, ...
430 'isSplit', 0, ...
431 'SplitLength', 4000));
432 % Get all the runs for this subject (ie the list of the study indices)
433 iStudyOther = setdiff(unique([sFilesAvg.iStudy]), sFilesAvg(1).iStudy);
434 % Copy the forward model file to the other runs
435 sHeadmodel = bst_get('HeadModelForStudy', sFilesAvg(1).iStudy);
436 for iStudy = iStudyOther
437 db_add(iStudy, sHeadmodel.FileName);
438 end
439
440 % ===== COMPUTE SOURCES: MEG =====
441 % Process: Compute sources [2018]
442 sAvgSrcMeg = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ...
443 'output', 1, ... % Kernel only: shared
444 'inverse', struct(...
445 'Comment', 'MN: MEG ALL', ...
446 'InverseMethod', 'minnorm', ...
447 'InverseMeasure', 'amplitude', ...
448 'SourceOrient', {{'fixed'}}, ...
449 'Loose', 0.2, ...
450 'UseDepth', 1, ...
451 'WeightExp', 0.5, ...
452 'WeightLimit', 10, ...
453 'NoiseMethod', 'reg', ...
454 'NoiseReg', 0.1, ...
455 'SnrMethod', 'fixed', ...
456 'SnrRms', 1e-06, ...
457 'SnrFixed', 3, ...
458 'ComputeKernel', 1, ...
459 'DataTypes', {{'MEG GRAD', 'MEG MAG'}}));
460 % Process: Snapshot: Sources (one time) - Loop only to get a correct comment for the report
461 for i = 1:length(sAvgSrcMeg)
462 bst_process('CallProcess', 'process_snapshot', sAvgSrcMeg(i), [], ...
463 'target', 8, ... % Sources (one time)
464 'orient', 4, ... % bottom
465 'time', 0.11, ...
466 'threshold', 20, ...
467 'Comment', ['MEG sources: ' sFilesAvg(i).FileName]);
468 end
469
470 % ===== COMPUTE SOURCES: EEG =====
471 % Process: Compute sources [2018]
472 sAvgSrcEeg = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ...
473 'output', 1, ... % Kernel only: shared
474 'inverse', struct(...
475 'Comment', 'MN: EEG', ...
476 'InverseMethod', 'minnorm', ...
477 'InverseMeasure', 'amplitude', ...
478 'SourceOrient', {{'fixed'}}, ...
479 'Loose', 0.2, ...
480 'UseDepth', 1, ...
481 'WeightExp', 0.5, ...
482 'WeightLimit', 10, ...
483 'NoiseMethod', 'reg', ...
484 'NoiseReg', 0.1, ...
485 'SnrMethod', 'fixed', ...
486 'SnrRms', 1e-06, ...
487 'SnrFixed', 3, ...
488 'ComputeKernel', 1, ...
489 'DataTypes', {{'EEG'}}));
490 % Process: Snapshot: Sources (one time) - Loop only to get a correct comment for the report
491 for i = 1:length(sAvgSrcEeg)
492 bst_process('CallProcess', 'process_snapshot', sAvgSrcEeg(i), [], ...
493 'target', 8, ... % Sources (one time)
494 'orient', 4, ... % bottom
495 'time', 0.11, ...
496 'threshold', 10, ...
497 'Comment', ['EEG sources: ' sFilesAvg(i).FileName]);
498 end
499 end
500 % Save report
501 ReportFile = bst_report('Save', []);
502 if ~isempty(reports_dir) && ~isempty(ReportFile)
503 bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_sources.html']));
504 end
505
506
507 %% ===== TIME-FREQUENCY =====
508 % Start a new report (one report for the time-frequency of all the subjects)
509 bst_report('Start');
510 % List of conditions to process separately
511 AllConditions = {'Famous', 'Scrambled', 'Unfamiliar'};
512 % Channels to display in the screen capture, by order of preference (if the first channel is bad, use the following)
513 SelChannel = {'EEG070','EEG060','EEG065','EEG050','EEG003'};
514 % Compute one separate time-frequency average for each subject/run/condition
515 for iSubj = 1:length(SubjectNames)
516 disp(sprintf('\n===== TIME-FREQUENCY: SUBJECT #%d =====\n', iSubj));
517 for iRun = 1:6
518 % Process: Select data files in: Subject/Run
519 sTrialsAll = bst_process('CallProcess', 'process_select_files_data', [], [], ...
520 'subjectname', SubjectNames{iSubj}, ...
521 'condition', sprintf('sub-%02d_ses-meg_task-facerecognition_run-%02d_proc-sss_meg_notch', iSubj, iRun));
522 % Loop on the conditions
523 for iCond = 1:length(AllConditions)
524 % Comment describing this average
525 strComment = [SubjectNames{iSubj}, ' / ', sprintf('run_%02d', iRun), ' / ', AllConditions{iCond}];
526 disp(['BST> ' strComment]);
527 % Find the first good channel in the display list
528 if isempty(BadChannels{iSubj}{iRun})
529 iSel = 1;
530 else
531 iSel = find(~ismember(SelChannel,BadChannels{iSubj}{iRun}), 1);
532 end
533 % Process: Select file comments with tag: Avg
534 sTrialsCond = bst_process('CallProcess', 'process_select_tag', sTrialsAll, [], ...
535 'tag', [AllConditions{iCond}, '_trial'], ...
536 'search', 1, ... % Search the file names
537 'select', 1); % Select only the files with the tag
538 % Process: Time-frequency (Morlet wavelets), averaged across trials
539 sTimefreq = bst_process('CallProcess', 'process_timefreq', sTrialsCond, [], ...
540 'sensortypes', 'MEG MAG, EEG', ...
541 'edit', struct(...
542 'Comment', ['Avg: ' AllConditions{iCond} ', Power, 6-60Hz'], ...
543 'TimeBands', [], ...
544 'Freqs', [6, 6.8, 7.6, 8.6, 9.7, 11, 12.4, 14, 15.8, 17.9, 20.2, 22.8, 25.7, 29, 32.7, 37, 41.7, 47.1, 53.2, 60], ...
545 'MorletFc', 1, ...
546 'MorletFwhmTc', 3, ...
547 'ClusterFuncTime', 'none', ...
548 'Measure', 'power', ...
549 'Output', 'average', ...
550 'RemoveEvoked', 0, ...
551 'SaveKernel', 0), ...
552 'normalize', 'none'); % None: Save non-standardized time-frequency maps
553 % Process: Extract time: [-200ms,900ms]
554 sTimefreq = bst_process('CallProcess', 'process_extract_time', sTimefreq, [], ...
555 'timewindow', [-0.2, 0.9], ...
556 'overwrite', 1);
557 % Screen capture of one sensor
558 hFigTf = view_timefreq(sTimefreq.FileName, 'SingleSensor', SelChannel{iSel});
559 bst_report('Snapshot', hFigTf, strComment, 'Time-frequency', [200, 200, 400, 250]);
560 close(hFigTf);
561 end
562 end
563 end
564 % Save report
565 ReportFile = bst_report('Save', []);
566 if ~isempty(reports_dir) && ~isempty(ReportFile)
567 bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_timefreq.html']));
568 end
569 end
570
571
572
573
574 %% ===== SUPPORT FUNCTIONS =====
575 function BadSeg = GetBadSegments(iSubj, iRun)
576 BadSegments{1} = {...
577 [247.867 248.185; 598.999 598.999; 598.999 598.999; 611.999 611.999; 612.999 612.999; 613.999 613.999; 616.999 616.999; 617.999 617.999; 623.999 623.999; 715.209 715.467], ...
578 [84.791 85.166], ...
579 [79.183 80.167], ...
580 [64.309 65.185], ...
581 [90.958 91.167; 178.005 178.355; 293.282 295.919; 312.298 316.479; 353.835 357.716], ...
582 [60.292 66.802; 69.975 71.210; 105.233 107.586; 108.822 109.506; 376.225 376.325]};
583 BadSegments{2} = {...
584 [223.806 224.199; 279.772 279.895; 453.241 455.108; 692.423 692.593], ...
585 [65.298 66.194; 304.727 306.178; 399.165 400.732], ...
586 [203.141 205.085; 281.579 287.883; 420.395 421.128], ...
587 [387.118 388.229; 440.318 441.900; 554.825 558.744], ...
588 [71.000 80.999; 82.750 87.367; 149.528 149.667; 264.747 267.995; 368.415 371.973; 376.263 378.763; 398.334 401.551; 537.410 541.645], ...
589 [38.000 47.999; 47.825 50.046; 61.298 61.384; 249.653 253.379; 282.917 283.820; 286.135 287.616; 298.167 300.196; 328.254 329.511; 335.957 337.817; 478.277 480.707]};
590 BadSegments{3} = {...
591 [406.312 407.207; 727.055 728.714], ...
592 [84.894 85.156; 152.028 152.946; 297.835 298.915; 418.272 421.845; 554.084 554.794], ...
593 [73.758 74.159; 378.212 378.536; 406.065 407.099; 470.541 471.698; 488.900 491.168; 529.596 530.453], ...
594 [94.874 95.152; 317.385 321.374; 325.696 327.055; 439.220 439.829; 454.473 455.175; 486.196 486.829; 518.660 522.015; 524.400 525.249; 562.417 570.325], ...
595 [96.208 97.181; 98.942 99.096; 135.005 135.754; 143.990 144.599; 250.139 250.247; 300.459 300.559; 338.265 339.322; 545.913 546.067], ...
596 [91.415 92.156; 284.843 286.525; 297.886 298.404; 317.046 317.163; 332.698 332.791; 358.946 359.402; 428.405 428.775; 478.374 478.690; 549.866 550.128]};
597 BadSegments{4} = {...
598 [22.967 22.967; 50.036 50.098; 52.058 52.058; 156.653 156.653; 171.565 173.386; 239.544 242.105; 268.162 270.175; 268.992 268.992; 316.032 316.032; 338.283 339.000; 357.959 361.909; 370.871 370.871; 381.579 383.677; 437.731 437.731; 463.482 468.505; 476.135 479.838; 486.652 488.272; 504.860 508.999], ...
599 [309.493 311.707; 342.681 344.525; 354.019 357.321; 390.023 391.225; 393.926 395.855; 404.221 405.069; 432.522 435.932; 459.048 460.715; 471.763 478.529; 549.387 551.999; 591.087 594.143; 608.541 611.079; 624.847 626.615; 649.648 651.570], ...
600 [57.411 58.198; 88.346 88.955; 200.761 202.335; 227.016 227.688; 257.726 258.054; 356.798 359.005; 404.260 411.003], ...
601 [46.000 54.823; 61.000 70.332; 203.005 207.125; 275.875 278.121; 313.500 314.824; 337.973 338.636; 422.505 426.239], ...
602 [58.000 62.479; 78.250 85.166; 89.955 91.360; 116.322 117.888; 130.013 131.987; 149.509 150.489; 174.650 175.823; 182.030 183.334; 196.758 197.384; 204.458 204.697; 205.236 208.663; 311.028 316.383; 320.700 327.181; 332.437 335.354; 344.205 346.133; 374.208 374.865; 385.519 386.214; 441.942 444.241; 453.957 456.997; 486.039 487.004; 501.238 504.185; 512.962 514.675; 553.398 556.215], ...
603 [41.406 45.743; 58.681 59.144; 108.086 108.896; 140.633 143.750; 196.110 199.474; 210.778 210.971; 234.649 235.143; 258.081 259.632; 339.101 340.805; 390.277 390.609; 438.935 442.122; 528.221 534.031]};
604 BadSegments{5} = {...
605 [265.539 265.778; 266.334 266.495; 268.479 268.965; 367.428 367.636; 439.655 442.779; 453.497 453.853; 504.997 505.329; 519.513 519.683; 595.674 595.982; 602.000 602.463], ...
606 [121.113 121.499; 124.971 126.213; 253.735 254.075; 272.232 272.464; 272.895 273.104; 346.368 346.645; 368.812 369.052; 406.382 406.605; 452.920 453.113; 454.903 455.112; 507.655 507.840; 508.766 509.013; 584.853 585.030; 594.831 595.656; 597.261 602.249], ...
607 [37.251 37.497; 38.825 39.056; 40.615 41.849; 43.624 44.758; 53.333 53.641; 54.698 55.076; 57.668 59.196; 79.129 79.360; 81.475 81.714; 122.658 123.375; 284.787 285.296; 288.754 288.993; 345.790 346.022; 421.212 421.459; 481.428 482.207; 503.408 503.685; 504.272 504.449; 524.451 524.714; 526.913 531.499], ...
608 [87.322 88.178; 91.085 91.325; 95.121 95.491; 114.174 114.397; 129.874 130.113; 151.220 151.544; 281.689 281.959; 532.966 533.345], ...
609 [59.176 60.218; 74.854 75.317; 308.180 309.877; 380.705 381.059], ...
610 [182.382 183.245; 196.220 196.736; 276.018 276.327; 292.490 294.086; 370.755 370.847; 435.644 436.624; 467.535 468.460; 522.838 525.847]};
611 BadSegments{6} = {...
612 [141.690 142.424; 157.070 157.417; 355.138 356.025; 423.999 423.999; 424.999 424.999; 426.999 426.999; 427.999 427.999; 486.430 488.151; 493.999 493.999; 501.511 501.619; 501.549 501.549; 501.585 501.585; 502.999 502.999; 503.999 503.999; 540.999 540.999; 541.999 541.999; 555.999 555.999; 556.999 556.999; 561.999 561.999; 563.999 563.999; 564.999 564.999; 565.999 565.999; 567.999 567.999], ...
613 [64.898 65.161; 71.700 72.718; 226.185 226.740; 324.124 324.425; 329.062 329.301; 486.143 486.975], ...
614 [62.300 63.148; 266.254 266.639; 409.920 410.221], ...
615 [54.048 55.214; 330.893 331.255], ...
616 [185.616 186.495; 331.411 331.796; 386.843 387.028; 387.999 387.999; 389.999 389.999; 434.575 434.875; 519.802 519.995], ...
617 [44.720 45.167; 211.446 211.964; 368.955 369.172]};
618 BadSegments{7} = {...
619 [154.966 155.144; 551.639 551.855], ...
620 [88.774 89.167; 107.999 107.999; 109.999 109.999; 110.999 110.999; 112.999 112.999; 113.999 113.999; 114.999 114.999; 119.999 119.999; 121.999 121.999; 124.223 125.465; 137.011 138.871], ...
621 [81.627 82.136; 377.953 381.046], ...
622 [241.136 242.171; 543.849 544.196; 596.639 598.553; 600.227 601.075; 603.999 603.999; 605.999 605.999; 609.999 609.999; 611.305 612.809; 614.999 614.999; 615.803 616.937; 623.694 625.877; 653.999 653.999; 655.055 655.756; 663.999 663.999], ...
623 [68.852 69.481; 74.034 75.200; 78.426 78.600; 104.497 105.963; 253.951 254.034; 256.964 257.038; 257.915 258.048; 323.254 324.156; 365.880 368.131; 369.952 370.060; 371.728 372.999; 430.931 431.965; 535.521 542.999], ...
624 [70.271 71.205; 94.441 96.445; 98.613 99.126; 112.318 112.749; 131.686 132.265; 148.935 150.615; 161.120 161.600; 205.325 208.214; 215.035 215.863; 217.403 218.818; 286.178 287.171; 399.075 404.853]};
625 BadSegments{8} = {...
626 [238.505 238.546; 256.354 257.224; 316.116 316.167; 341.519 341.558; 356.493 356.566; 380.095 380.170; 391.906 392.048; 448.457 448.562; 469.931 470.028; 530.488 530.575; 555.180 558.136; 562.152 562.245; 588.403 588.502; 625.205 625.662; 638.019 638.438; 649.982 650.008; 650.691 651.357; 651.925 652.061; 665.445 665.472; 695.923 696.015; 706.528 706.720; 729.706 732.534], ...
627 [99.546 106.114; 113.245 114.199; 150.885 154.989; 277.486 278.075; 333.645 335.574; 339.358 340.646; 371.860 375.394; 497.459 499.156], ...
628 [49.008 50.205; 231.446 233.406; 329.590 329.659; 355.101 356.019; 360.733 360.973; 372.955 374.891; 389.283 392.068; 453.610 455.431; 464.632 465.265; 489.209 489.996; 514.777 515.295], ...
629 [178.368 179.232; 293.865 294.521; 418.252 418.314; 450.124 450.209; 480.236 480.445; 492.725 493.975; 495.170 497.624; 500.382 500.459; 504.247 504.402; 628.017 628.079; 628.827 628.905; 630.209 630.332], ...
630 [100.616 101.172; 107.691 108.539; 188.814 188.875; 193.119 193.281; 207.964 208.033; 432.254 432.385; 489.834 489.911; 517.960 518.037; 518.955 519.040; 520.938 521.062; 521.945 522.037; 523.959 524.052; 525.942 526.057; 526.945 527.015; 528.958 529.059; 531.138 531.192; 531.979 532.048; 532.951 533.028; 533.962 534.024], ...
631 [135.997 137.232; 383.565 383.657; 418.763 418.955]};
632 BadSegments{9} = {...
633 [215.107 216.187; 262.388 262.388; 287.519 287.635; 289.895 290.135; 311.999 311.999; 350.161 351.179; 526.154 527.033; 564.999 564.999; 584.935 585.059; 587.999 587.999; 601.999 601.999; 603.999 603.999; 608.999 608.999; 612.999 612.999], ...
634 [47.667 47.775; 49.172 54.681; 67.195 70.805; 78.988 80.477; 138.701 138.948; 138.727 138.727; 138.728 138.728; 138.728 138.728; 138.734 138.734; 138.734 138.734; 138.881 138.881; 138.881 138.881; 138.907 138.907; 140.993 141.093; 141.037 141.037; 141.045 141.045; 155.795 155.864; 155.822 155.822; 155.849 155.849; 168.999 168.999; 206.999 206.999; 219.999 219.999; 225.999 225.999; 226.999 226.999; 228.999 228.999; 236.999 236.999; 242.463 242.463; 242.463 242.463; 242.487 242.487; 247.999 247.999; 251.999 251.999; 252.483 252.483; 253.315 254.163; 265.999 265.999; 267.999 267.999; 272.358 272.358; 272.410 272.410; 298.999 298.999; 300.999 300.999; 314.037 314.037; 314.047 314.047; 321.813 321.813; 321.813 321.813; 321.833 321.833; 329.117 329.117; 329.117 329.117; 329.156 329.156; 346.999 346.999; 347.999 347.999; 349.320 349.320; 349.329 349.329; 352.528 355.869; 364.040 364.040; 396.278 397.420; 404.865 404.865; 407.905 407.905; 407.905 407.905; 407.954 407.954; 418.454 418.454; 418.454 418.454; 418.486 418.486; 441.999 441.999; 444.999 444.999; 447.999 447.999; 453.550 453.650; 454.931 455.055; 457.999 457.999; 479.964 480.079; 481.999 481.999; 482.949 483.073; 488.948 489.064; 511.999 511.999; 520.999 520.999; 523.999 523.999; 526.954 527.069; 533.999 533.999; 537.999 537.999], ...
635 [82.889 83.198; 206.143 207.424; 403.873 404.096; 406.790 407.485; 413.936 414.075; 415.912 416.112; 536.621 537.532], ...
636 [182.999 182.999; 182.999 182.999; 183.999 183.999; 195.999 195.999; 208.999 208.999; 209.999 209.999; 278.601 278.955; 413.913 414.067; 415.250 417.792; 419.940 420.087; 420.999 420.999; 512.999 512.999; 514.363 516.254; 521.917 522.134; 522.999 522.999; 523.915 524.101; 533.999 533.999; 538.999 538.999; 539.892 540.062; 543.999 543.999; 548.874 549.089; 549.900 550.062; 552.755 554.075; 585.999 585.999; 587.957 588.073; 588.999 588.999; 594.999 594.999; 603.971 604.056; 604.928 605.097; 615.986 618.495; 623.999 623.999], ...
637 [53.215 54.203; 164.227 168.965; 201.433 202.775; 292.973 295.889; 303.961 304.817; 309.354 311.236; 313.123 313.639; 322.547 323.017; 331.141 334.852; 356.637 356.985; 367.855 368.079; 369.995 373.459; 377.849 378.142; 406.072 407.306; 436.164 436.665; 458.033 458.905; 516.889 518.656; 517.812 518.684], ...
638 [113.945 115.225; 198.570 198.863; 264.162 264.795; 383.705 385.641; 396.477 397.064; 399.706 406.457; 452.261 453.179; 486.338 487.102; 498.869 499.579; 507.968 508.925; 546.266 547.285; 558.825 560.128]};
639 BadSegments{10} = {...
640 [235.104 236.184; 324.272 325.028; 330.799 331.401; 541.062 541.826; 564.747 565.072], ...
641 [100.581 101.229; 285.644 285.644; 285.644 285.644; 297.979 298.033; 298.999 298.999; 300.949 301.026; 301.999 301.999; 303.999 303.999; 304.999 304.999; 306.999 306.999; 307.999 307.999; 310.999 310.999; 311.999 311.999; 313.971 314.025; 314.999 314.999; 316.995 317.035; 317.999 317.999; 319.999 319.999; 320.923 321.046; 326.971 327.048; 327.999 327.999; 329.999 329.999; 330.999 330.999; 332.974 333.028; 333.999 333.999; 335.971 336.025; 336.999 336.999; 338.973 339.035; 339.999 339.999; 341.951 342.044; 343.999 343.999; 344.964 345.033; 346.999 346.999; 347.999 347.999; 452.015 453.411; 458.735 459.776; 467.974 468.089], ...
642 [41.266 41.675; 53.905 55.240; 159.159 159.345; 220.255 220.394], ...
643 [66.751 67.190; 201.674 201.812; 294.119 295.168; 303.188 303.528; 316.992 317.037; 317.999 317.999; 317.999 317.999; 319.999 319.999; 320.980 321.042; 322.999 322.999; 325.999 325.999; 326.955 327.048; 328.999 328.999; 333.999 333.999; 334.999 334.999; 342.999 342.999; 343.999 343.999; 351.915 352.038; 352.999 352.999; 407.979 408.056; 409.999 409.999; 411.999 411.999; 415.405 416.154; 435.999 435.999; 436.962 437.046; 468.885 469.100; 469.999 469.999; 470.999 470.999; 516.210 516.357], ...
644 [105.369 106.172; 141.468 142.178; 199.008 199.818; 201.269 201.716; 352.851 353.083; 449.865 452.041; 459.508 459.802], ...
645 [229.948 230.033; 230.999 230.999; 230.999 230.999; 259.588 259.990; 343.970 345.798; 361.999 361.999; 362.970 363.046; 364.999 364.999; 366.999 366.999; 367.962 368.031; 372.521 374.179; 405.999 405.999; 409.999 409.999; 411.999 411.999; 412.999 412.999; 434.999 434.999; 435.970 436.031; 516.999 516.999; 519.999 519.999]};
646 BadSegments{11} = {...
647 [179.045 180.225; 323.999 323.999; 323.999 323.999; 324.999 324.999; 327.978 328.025; 328.966 329.028; 330.999 330.999; 365.999 365.999; 367.999 367.999; 370.999 370.999; 371.999 371.999; 373.999 373.999; 375.965 376.042; 377.999 377.999], ...
648 [52.107 53.156; 521.077 521.155], ...
649 [65.331 66.156], ...
650 [81.579 82.143; 108.565 108.565; 108.566 108.566; 112.176 112.176; 122.377 122.377; 122.565 122.565; 213.882 213.882; 215.305 215.305; 224.851 224.851; 224.919 224.919; 255.815 255.815; 257.893 257.893; 359.952 359.952; 361.630 361.754; 370.285 370.285; 376.853 376.853; 511.600 511.600; 513.631 513.631; 513.988 513.988; 518.215 518.215], ...
651 [63.205 64.154; 170.951 171.036; 176.841 176.841; 176.842 176.842; 177.282 177.282; 223.886 223.971; 259.963 259.963; 261.547 261.547; 341.185 341.302; 368.999 368.999; 370.999 370.999; 374.999 374.999; 382.971 383.033; 383.999 383.999; 386.999 386.999; 388.958 389.035; 390.999 390.999; 391.955 392.031; 394.955 395.040; 396.999 396.999; 398.999 398.999; 400.999 400.999; 402.999 402.999; 404.963 405.033; 406.999 406.999; 429.829 429.829; 430.346 430.346; 465.184 465.184; 470.290 470.290], ...
652 [50.195 51.145; 158.850 159.012]};
653 BadSegments{12} = {...
654 [193.435 194.175; 500.000 501.000; 08.988 509.868; 528.999 528.999; 528.999 528.999; 529.999 529.999; 531.999 531.999; 547.999 547.999], ...
655 [133.484 134.185; 134.911 135.065; 553.999 553.999; 553.999 553.999; 554.999 554.999; 557.999 557.999; 558.999 558.999; 564.999 564.999; 565.977 566.046; 568.999 568.999; 569.951 570.028; 571.999 571.999; 579.959 580.028; 580.977 581.055; 583.971 584.033; 583.999 583.999; 585.999 585.999; 586.999 586.999; 588.980 589.026; 590.999 590.999; 591.958 592.044; 595.999 595.999], ...
656 [46.028 47.216; 129.557 130.236; 200.999 200.999; 200.999 200.999; 201.975 202.075; 203.999 203.999; 204.967 205.053; 213.627 214.515; 218.958 219.044; 222.959 228.252; 309.999 309.999; 311.897 317.537; 311.999 311.999; 351.968 353.234; 399.999 399.999; 446.508 451.454; 487.999 487.999; 512.278 512.926], ...
657 [82.535 84.062; 102.247 102.355; 134.999 134.999; 134.999 134.999; 136.999 136.999; 138.973 139.042; 147.999 147.999; 148.999 148.999; 149.999 149.999; 193.999 193.999; 194.999 194.999; 224.984 226.095; 241.954 242.077; 243.948 244.087; 280.999 280.999; 281.999 281.999; 305.000 305.451; 310.945 311.068; 328.999 328.999; 352.999 352.999; 353.999 353.999; 397.025 397.411; 415.944 418.096; 415.999 415.999; 470.999 470.999; 477.971 478.041; 483.997 484.737; 492.958 493.089; 493.999 493.999; 494.999 494.999; 522.999 522.999; 523.999 523.999; 524.955 525.039], ...
658 [175.183 176.155; 246.112 246.991; 412.005 413.340; 483.915 484.061; 485.959 486.028; 497.360 498.880; 594.457 594.944; 596.641 598.006; 615.277 618.710], ...
659 [66.900 72.232; 75.510 78.688; 543.495 545.999]};
660 BadSegments{13} = {...
661 [307.246 308.179; 627.881 628.089; 629.941 630.049], ...
662 [171.506 172.301; 172.999 172.999; 430.999 430.999; 432.999 432.999; 434.999 434.999; 448.999 448.999; 479.999 479.999; 489.999 489.999; 490.999 490.999; 491.999 491.999], ...
663 [92.059 93.209], ...
664 [52.791 53.231; 54.999 54.999; 65.985 68.138; 91.240 91.386; 92.137 92.207; 105.346 105.493; 121.120 121.398; 146.546 146.955; 194.025 197.652; 242.985 243.571; 247.101 247.556; 269.889 270.083; 270.946 271.070; 274.408 275.866; 294.565 295.267; 319.370 319.879; 365.999 365.999; 375.913 376.044; 376.869 377.055; 377.999 377.999; 390.130 393.225; 403.965 404.035; 404.953 405.075; 419.265 419.565; 427.015 427.154; 427.879 428.118; 427.999 427.999; 428.913 429.067; 480.945 482.095; 484.285 484.571; 516.929 517.099; 517.963 518.079], ...
665 [115.921 116.060; 117.999 117.999; 120.999 120.999; 126.896 127.058; 130.955 131.039; 131.957 132.043; 132.961 133.015; 134.999 134.999; 135.950 136.043; 139.999 139.999; 143.999 143.999; 144.850 145.066; 145.999 145.999; 146.999 146.999; 160.929 161.045; 161.999 161.999; 189.086 189.310; 193.939 194.063; 195.975 196.045; 199.961 200.045; 201.952 202.059; 202.939 205.084; 210.999 210.999; 211.738 213.629; 217.954 218.061; 218.987 219.049; 219.975 220.028; 258.953 259.037; 259.952 260.045; 260.963 261.032; 261.958 262.043; 262.954 263.054; 262.999 262.999; 288.999 288.999; 295.477 295.655; 298.976 299.054; 299.979 300.018; 300.975 301.005; 336.952 337.044; 337.977 338.032; 338.988 339.027; 338.999 338.999; 339.991 340.045; 340.975 341.029; 341.978 342.016; 343.976 344.054; 361.961 362.045; 365.942 366.050; 370.999 370.999; 384.961 385.045; 409.977 410.047; 415.984 416.030; 426.944 427.059; 430.999 430.999; 436.999 436.999; 448.999 448.999; 449.500 452.131; 481.962 482.055; 482.973 483.027; 484.975 485.060; 503.999 503.999; 504.999 504.999; 507.967 508.028; 508.955 509.032; 510.999 510.999; 512.999 512.999; 518.865 520.292; 522.958 523.043; 523.999 523.999; 524.948 525.034; 525.959 526.044; 526.946 527.055; 527.942 528.027; 528.961 529.053; 529.945 530.022; 535.999 535.999; 536.957 537.058; 538.967 539.036; 539.963 540.040; 541.976 542.045; 543.999 543.999; 551.956 552.034; 552.990 553.028; 553.977 554.024; 555.968 556.022; 556.204 557.454; 557.963 558.032; 559.969 560.023; 563.973 564.035; 563.999 563.999], ...
666 [87.308 87.639; 107.039 108.189; 427.943 428.035; 437.965 438.027; 439.973 440.026; 491.275 492.571]};
667 BadSegments{14} = {...
668 [365.982 367.194; 368.999 368.999; 368.999 368.999; 369.999 369.999; 371.999 371.999; 373.999 373.999; 375.942 376.050; 376.999 376.999; 378.999 378.999; 380.999 380.999; 382.999 382.999; 383.999 383.999; 386.999 386.999; 387.999 387.999; 418.999 418.999; 444.999 444.999; 690.178 693.218], ...
669 [101.796 102.183; 218.999 218.999; 219.999 219.999; 221.999 221.999; 222.999 222.999; 227.999 227.999; 229.999 229.999; 230.999 230.999; 232.999 232.999; 233.999 233.999; 235.999 235.999; 237.999 237.999; 238.985 239.031; 240.999 240.999; 242.999 242.999; 243.999 243.999; 246.999 246.999; 247.999 247.999; 253.999 253.999; 258.999 258.999; 259.999 259.999; 260.227 261.123; 264.999 264.999; 265.999 265.999; 268.999 268.999; 269.999 269.999; 274.999 274.999; 276.999 276.999; 277.999 277.999; 280.999 280.999; 282.966 283.021; 284.999 284.999; 285.999 285.999; 360.999 360.999; 362.973 363.026; 364.999 364.999; 365.999 365.999; 367.999 367.999; 368.985 369.024; 370.999 370.999; 371.999 371.999; 378.999 378.999; 379.999 379.999; 381.999 381.999; 383.999 383.999; 385.999 385.999; 387.999 387.999; 388.999 388.999; 398.999 398.999; 415.999 415.999; 417.999 417.999; 418.999 418.999; 421.999 421.999; 422.999 422.999], ...
670 [89.285 90.187], ...
671 [85.991 87.179; 88.999 88.999; 89.999 89.999; 92.966 93.027; 93.999 93.999; 107.999 107.999; 110.999 110.999; 111.971 112.025; 118.999 118.999; 121.999 121.999; 122.999 122.999; 166.999 166.999; 168.957 169.035; 173.999 173.999; 174.999 174.999; 175.957 176.035; 198.999 198.999; 199.999 199.999; 202.999 202.999; 203.999 203.999; 205.999 205.999; 207.999 207.999; 208.999 208.999; 210.974 211.044; 211.999 211.999; 212.999 212.999; 232.984 233.023; 234.999 234.999; 235.999 235.999; 236.999 236.999; 239.958 240.044; 240.999 240.999], ...
672 [71.132 72.227], ...
673 [98.134 98.827; 110.497 110.814; 114.000 117.272; 279.999 279.999; 279.999 279.999; 280.999 280.999; 282.999 282.999; 283.999 283.999; 284.999 284.999; 286.999 286.999; 287.999 287.999; 288.999 288.999; 289.999 289.999; 291.966 292.044; 292.999 292.999; 346.999 346.999; 437.608 437.646; 518.988 519.104; 522.969 523.015; 534.999 534.999; 536.967 537.029; 562.725 563.426]};
674 BadSegments{15} = {...
675 [66.755 67.172; 257.988 258.042; 258.714 259.385; 260.958 261.044; 265.931 266.062; 267.968 268.022; 268.978 269.033; 270.992 271.031; 272.964 273.017; 273.968 274.030; 275.975 276.059; 276.970 277.046; 288.976 289.038; 289.984 290.015; 295.994 296.025; 296.982 297.013; 300.966 301.152; 364.730 364.800], ...
676 [46.531 46.555; 65.156 66.148; 87.535 87.642; 183.896 183.934; 294.342 294.375; 330.004 330.528; 333.785 333.815; 345.177 345.203; 399.544 399.573; 480.964 480.980], ...
677 [53.557 54.175; 148.851 149.136], ...
678 [43.681 44.214; 361.255 361.587; 409.462 411.021], ...
679 [90.868 91.130; 410.935 411.044], ...
680 [64.786 65.141; 132.385 132.786; 177.735 178.252; 278.747 278.902; 313.959 314.028; 314.970 315.031]};
681 BadSegments{16} = {...
682 [38.958 44.335; 305.280 312.826; 324.778 326.067; 393.183 398.367; 403.410 422.854], ...
683 [49.726 50.251; 58.675 65.465; 263.918 264.088; 321.795 322.195; 522.881 523.081], ...
684 [61.245 62.325; 65.628 65.905; 315.686 317.275; 335.184 336.449; 387.870 388.117; 389.892 390.062], ...
685 [56.020 57.255; 60.711 61.096; 64.353 64.924; 73.202 73.757; 76.134 76.844; 95.625 96.859; 171.436 172.625; 178.359 178.715; 212.828 213.415; 352.425 352.617; 367.755 369.097; 488.400 489.426], ...
686 [71.665 72.513; 76.347 77.181; 127.134 128.862], ...
687 [46.190 47.209; 53.461 55.752; 194.317 194.510; 213.101 213.240; 216.356 216.442; 225.962 226.286; 245.163 245.464; 252.041 253.422; 254.648 254.856; 292.966 294.340; 340.928 341.067; 346.298 346.452]};
688 BadSeg = BadSegments{iSubj}{iRun}';
689 end
690
You should note that this is not the result of a fully automated procedure. The bad channels were identified manually and are defined for each run in the script. The bad segments were detected automatically, confirmed manually for each run and saved in external files, then exported as text ans copied at the end of this script.
All the process calls (bst_process) were generated automatically using with the script generator (menu Generate .m script in the pipeline editor). Everything else was added manually (loops, bad channels, file copies).