Source reconstruction with pymeg

Pymeg provides a host of helper functions to facilitate analysis of MEG data. It makes heavy use of python-mne and mainly automates specific tasks and provides the glue between separate processing steps. It does not actually provide a lot of new functionality.

This guide focuses exclusively on carrying out source reconstruction using an LCMV beamformer. It describes most of the necessary preprocessing steps.

The core idea of pymeg is that you want to carry out single-trial source reconstruction in a set of regions of interest. Anything else will require various amounts of extra work not covered here.

Preparing a subject for source reconstruction

Getting a subject ready for source reconstruction requires a few steps that can be carried out in parallel:

• Preprocessing/epoching
• Surface generation / recon-all

After this you can continue by:

• Aligning the HCPMMP atlas to your subject
• Aligning the Wang et al. 2015 atlas

(or in general any kind of labels that you need) - Creating a high density head shape model - Creating BEM surfaces with fieldtrip - Creating a MRI<>MEG transformation matrix

Preprocess MEG data

Use MNE/pymeg or fieldtrip to preprocess and epoch your data. Make sure that you do not remove head localization channels from your epochs.

Please save your epochs as MNE fif files. Pymeg contains a fieldtrip->MNE conversion script to load saved fieldtrip structures and save them as fif files. Please note that this requires the raw file that generated each epoch file - this allows us to extract MNE specific info structures.

It is also advisable to not change the channel names if you use fieldtrip, if you do matching of fieldtrip channels to MNE channel info might become difficult.

The conversion script can be found in pymeg/read_fieldtrip.py

Surface generation / recon-all

MNE performs source reconstruction on cortical surfaces. It is easiest to create these surfaces by using freesurfers recon-all pipeline.

If you are @UKE you can use the cluster wide install for this. The first step is to set your freesurfer subject folder:

> export SUBJECTS_DIR=/path/to/subjects_dir

This directory will host the output of freesurfer (e.g. all cortical meshes). Freesurfer will automatically create a subdirectory for each subject here.

At this point you should also copy freesurfers fsaverage and fsaverage_sym from /opt/progs/freesurfer/subjects to your subject dir. We will later add new information for these subjects and linking to /opt/progs makes these subjects write protected.

You will likely have gotten a single T1 MPRAGE image in DICOM format as each subjects MRI. This needs to be converted to NIFTI format. You can use dcm2niix for this (installed in /home/nwilming/bin/). The correct DICOM series will have ‘mprage’ in it’s name once converted to NIFTI.

Now call recon-all:

> recon-all -subjid SUBID -all -i path/to/nifti_T1

And wait. This will take a few hours. It makes sense to parallelize this call over the cluster.

Aligning the HCPMMP atlas to your subject

The cluster wide freesurfer install is quite old. There is a newer version available in /home/nwilming/freesurfer. Setup your bashrc to reflect this install.

To then get the HCPMMP atlas fire up a python >= 3.6 interpreter and import pymeg:

>>> from pymeg import atlas_glasser as ag
>>> ag.get_hcp('subjects_dir') # Needs to be done only once
>>> ag.get_hcp_annotation(subjects_dir, subject)

Now check if this worked by loading the subjects inflated surface and the HCOMMP annotation file in freeview (only works on node031).

Aligning the Wang et al. 2015 atlas

This is a bit more involved. I suggest you use docker. Follow the instructions here:

…. link missing …

Docker does not work on the cluser! You need to do this on your own laptop.

Creating a high density head shape model

Create a scalp surface from the MRI:

>>> mne make_scalp_surfaces -s SUBID -d SUBJECT_DIR

Creating BEM surfaces with fieldtrip

MNE has a tool to create bem surfaces from an MRI (mne watershed_bem). I found that this did not work particularly well with the mprage based images that we get. Using SPM/fieldtrip produced less artifacts.

To create head meshes with fieldtrip use the scripts in mfiles/

Creating a MRI<>MEG transformation matrix

Use the MNE GUI to create transformation matrices. Pymeg provides a little help here… XXX

pymeg.source_reconstruction.add_volume_info(subject, surface, subjects_dir, volume='T1')

Add volume info from MGZ volume

pymeg.source_reconstruction.check_bems(subjects)

Create a plot of all BEM segmentations.

pymeg.source_reconstruction.circumcenter(coil1, coil2, coil3)

Determines position and orientation of the circumcenter of fiducials. Adapted from: http://www.fieldtriptoolbox.org/example/how_to_incorporate_head_movements_in_meg_analysis CIRCUMCENTER determines the position and orientation of the circumcenter of the three fiducial markers (MEG headposition coils).

Parameters:coil1-3 – 3xN array X,y,z-coordinates of the 3 coils [3 X N],[3 X N],[3 X N] where N is timesamples/trials. Returns:X,y,z-coordinates of the circumcenter [1-3 X N], and the orientations to the x,y,z-axes [4-6 X N].

1. Stolk, 2012

pymeg.source_reconstruction.get_ctf_trans(filename)

Get transformation matrix between sensors and head space.

pymeg.source_reconstruction.head_movement(epochs)

Compute head movement from epochs.

Returns the circumcenter of the three fiducials for each time point.

pymeg.source_reconstruction.make_bem_model(subject, ico=4, conductivity=(0.3, 0.006, 0.3), subjects_dir=None, verbose=None, bem_sub_path='bem')

Create a BEM model for a subject.

Copied from MNE python, adapted to read surface from fieldtrip / spm segmentation.

pymeg.source_reconstruction.make_trans(subject, raw_filename, epoch_filename, trans_name)

Create coregistration between MRI and MEG space.

Call MNE gui to create a MEG<>MRI transformation matrix

pymeg.source_reconstruction.replace_fiducials(info, fiducials)

Replace initial fiducial measuremnt with new estimates

CTF systems measure fiducial location at the beginning of the measurement. When used with online head loc over multiple sessions these measurements are not accurate. This is because subjects are guided to the head position of previous sessions.

Parameters:

• info – MNE info structure
• fiducials – dict Dictionary that contains fiducial positions, e.g. see output of get_ref_head_pos.

Returns:

Info structure with updated head position.

pymeg.source_reconstruction.set_fs_subjects_dir(directory)

Set freesurfer subjectdir environment variable

pymeg.lcmv.accumulate(data, time, est_key, est_val, roi, trial)

Transform source reconstruction results to a DataFrane.

Parameters:

• data – ndarray If ntrials x vertices x time in which case the function will average across vertices. If ntrials x time will be directly converted to df.
• time – ndarray time points that match last dimension of data
• est_key – value estimation key for this value
• est_val – value estimation value for this set of data
• roi – str Name of the roi that this comes from
• trial – ndarray Needs to match first dim of data

Returns:

A pandas DataFrame that contains source reconstructed data with hierarchical index to describe each data point.

pymeg.lcmv.apply_lcmv(tfrdata, est_key, est_vals, events, times, info, filters, post_func=None, accumulate_func=None, max_ori_out='signed')

Apply Linearly Constrained Minimum Variance (LCMV) beamformer weights.

Parameters:

• tfrdata – ndarray Data to be reconstructed. Should be either n_trials x n_sensors x Y x n_time or trials x sensors x time. Reconstruction treats epochs and dim Y as independent dimensions.
• est_key – value A key to identify this reconstruction (e.g. F for power)
• est_vals – sequence Values that identify different reconstructions along dimension Y for a single epoch, e.g. the frequency for power reconstructions. Needs to be length Y.
• events – array Identifiers for different epochs. Needs to be of length n_trials.
• times – array Time of entries in last dimension of input data.
• info – mne info structure Info structure of the epochs which are to be reconstructed
• filters – dict Contains ROI names as keys and MNE filter dicts as values.
• post_func – function This function is applied to the reconstructed epochs, useful to convert complex TFR estimates into power values.
• accumulate_func – function Function that is applied after post_func has been applied. Can for example be used to transform the output into a dataframe.
• max_ori_out – str, default ‘signed’ This is passed to the MNE LCMV function which at the moment requires this to be ‘signed’

Returns:

List of source reconstructed epochs transformed by post_func.

pymeg.lcmv.broadband_est(x, time, est_val=[-1], est_key='BB', **kwargs)

Estimate broadband from source reconstructed time course

pymeg.lcmv.complex_tfr(x, time, est_val=None, est_key=None, sf=600.0, foi=None, cycles=None, time_bandwidth=None, n_jobs=1, decim=10)

Estimate power of epochs in array x.

pymeg.lcmv.flip_and_avg_vertices(data)

Correct random sign flips in reconstructed vertices.

Average over vertices but correct for random flips first Correction is done by ensuring positive correlations between vertices

Parameters:data – ndarray A 2D array with vertices in the first dimension and time in the second. Returns:A single time series constructed by averaging over vertices.

pymeg.lcmv.get_filter(info, forward, data_cov, noise_cov, label=None, reg=0.05, pick_ori='max-power')

Comput LCMV filter for one region of interest.

pymeg.lcmv.get_filters(estimator, epochs, forward, source, noise_cov, data_cov, labels)

Compute LCMV filters for a list of regions of interest.

pymeg.lcmv.par_reconstruct(pre_estimator, pre_est_args, epochs, events, times, info, filters, post_func=<function tfr2power_estimator>, accumulate_func=<function accumulate>, njobs=4)

Source reconstruct epochs with flexible transform before and after.

This function performs source reconstruction and can transform the input data before reconstruction (e.g. for TFR) and after reconstruction (e.g. to compute power from complex FFT output). Output data can be passed through yet another function to shape into the desired output.

The flow of data through this function therefore is:

epochs -> pre_estimator -> to source space ->

post_func -> accumulate_func

Parameters:

• pre_estimator – function A function that is applied to the sensor space data before source reconstruction. Use ‘complex_tfr’ to project the sensor space TFR representation into source space.
• pre_est_args – dict A dict that is **pased to pre_estimator, e.g. additional arguments to customize behavior of this function.
• epochs – ndarray Data array of MNE epochs object
• events – ndarray Vector that assigns unique identifiers to different epochs in data array.
• times – ndarray Vector that assigns time points to time dimension in data array
• info – MNE info object
• filters – Dictionary returned by setup_filters
• post_func – function A function that is applied to the source reconstructed data. To get TFR in source space you can pass ‘complex_tfr’ as ‘estimator’ and then use ‘tfr2power_estimator’ here to compute power from complex FFT output.
• accumulate_func – function A function that takes the output of post func, the estimation keys, estimation values, time points, the region of interest and trial identifiers as inputs and returns a pandas DataFrame.
• njobs – int Number of cores to parallelize over.

Returns:

List that contains output for each ROI.

pymeg.lcmv.reconstruct_broadband(filters, info, epochs, events, times, estimator=<function broadband_est>, est_args={}, post_func=None, accumulate_func=<function accumulate>, njobs=4)

Reconstruct broadband activity from a set of regions of interest.

Parallelization is applied across filters, i.e. regions of interest. This function calls par_reconstruct with appropriate default settings for broadband reconstruction.

See reconstruct_tfr for description of arguments

reconstruct_tfr(filters, info, epochs, events, times, estimator=<function complex_tfr>, est_args={'cycles': array([ 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. , 4.5, 5. ,

5.5, 6. , 6.5, 7. , 7.5, 8. , 8.5, 9. , 9.5,

10. , 10.5, 11. , 11.5, 12. , 12.5, 13. , 13.5, 14. , 14.5]), 'est_key': 'F', 'est_val': array([ 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70,

75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130, 135,

140, 145]), 'foi': array([ 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70,

75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130, 135,

140, 145]), 'n_jobs': 1, 'time_bandwidth': 2}, post_func=<function tfr2power_estimator>, accumulate_func=<function accumulate>, njobs=4)

Reconstruct time frequency representation of epochs.

Parallelization is applied across filters, i.e. regions of interest. This function calls par_reconstruct with appropriate default settings for TFR reconstruction. Change the est_args argument to specify parameters for the time frequency conversion.

Parameters:

• filters – Dictionary returned by setup_filters
• info – MNE info object
• epochs – ndarray Data array of MNE epochs object
• events – ndarray Vector that assigns unique identifiers to different epochs in data array.
• times – ndarray Vector that assigns time points to time dimension in data array
• estimator – function A function that is applied to the sensor space data before source reconstruction. Use ‘complex_tfr’ to project the sensor space TFR representation into source space.
• est_args – dict A dict that is **pased to estimator, e.g. parameters for the TFR transformation.
• post_func – function A function that is applied to the source reconstructed data. To get TFR in source space you can pass ‘complex_tfr’ as ‘estimator’ and then use ‘tfr2power_estimator’ here to compute power from complex FFT output.
• accumulate_func – function A function that takes the output of post func, the estimation keys, estimation values, time points, the region of interest and trial identifiers as inputs and returns a pandas DataFrame.
• njobs – int Number of cores to parallelize over.

Returns:

Concatenated outputs across regions of interest.

pymeg.lcmv.tfr2power_estimator(x)

Compute power on source reconstructed FFT results

Indices and tables

• Index
• Module Index
• Search Page