-
Notifications
You must be signed in to change notification settings - Fork 26
MEG preprocessing
Only resting state for now... Check the example here.
For MEG resting state analysis, first, we need to create a time course per cortical label. This can be done using one call to the MEG preprocessing script:
python -m src.preproc.meg -s subject-name -a atlas-name -t rest -f rest_functions --l_freq low-freq --h_freq high-freq --windows_length window-length --windows_shift windows-shift
Where low-freq and high-freq are the band limit for filtering. Instead of reconstructing the source for all the raw data at once, the script divides the raw data into epochs (windows), and reconstruct the time course separately for each window. In case you are interested only in calculating a global measurement of the resting state data, you can set the windows-shift to be equal to the window-length (10s for example), and then just concatenate the windows. If you are interested also in calculating the resting state dynamics, set the windows-shift to a 1/n of the window-length, and calculate the connectivity measurement for each window. The consecutive steps are the following:
- Raw data (sensors level) filtering. On default, the raw data is also filtered to remove the power line noise using a notch filter (The default frequency is 60Hz. You can change it using the power_line_freq flag). If you don't want to remove the power line noise, set the remove_power_line_noise flag to 0. The raw data can be filtered also by setting the flags low-freq and high-freq:
- low_freq < high_freq: band-pass filter
- low_freq > high_freq: band-stop filter
- only low_freq: high-pass filter
- only h_freq: low-pass filter
- The data is segmented into windows of window-length length (ms) and windows-shift (ms). On default, certain windows can be rejected according to the flags reject_grad, reject_mag, and reject_eog. If you prefer to take the risk and not to reject any of the windows, set the reject to 0.
- The forward model and inverse operator are calculated (if they don't exist)
- The source is reconstructed for each window separately using the selected inverse method. The default is dSPM.
- For each label, according to the atlas-name (the default is aparc.DKTatlas40), the time course is extracted using the selected extract-mode (the default is mean_flip)
For example, you wish to calculate the resting state dynamics of 0.5s windows length in the alpha band, using laus125 as atlas, MNE as inverse method and calculate the pca-flip of each label:
python -m src.preproc.meg -s subject-name -a laus125 -t rest -f rest_functions -i MNE --l_freq 8 --h_freq 10 --windows_length 500 --windows_shift 100 --extract_mode mean_flip
The end result of this call is two files, labels_data_rh.npz, and labels_data_lh.npz. Each one contains two variables, data (numpy matrix, #cortical labels x T), and names, which is the cortical labels names. If you wish to run a resting state analysis, take a look here.
If you wish to call this script from python, take a look at the example in src.preproc.examples.meg, function calc_rest:
from src.preproc import meg
from src.utils import utils
args = meg.read_cmd_args(dict(
subject=args.subject,
mri_subject=args.mri_subject,
atlas='laus125',
function='rest_functions',
task='rest',
l_freq=3, h_freq=80,
windows_length=500,
windows_shift=100,
inverse_method='MNE',
))
call_main(args)