MEG: Preprocessing, Source Reconstruction, Parcellation

This tutorial walks through the full MEG processing pipeline step by step:

1. Preprocessing — Filter, downsample, detect bad segments/channels, ICA artefact rejection.

2. Surface Extraction — Extract skull/scalp surfaces from a structural MRI.

3. Coregistration — Align MEG sensor space to MRI space.

4. Forward Model — Compute the lead field matrix.

5. Source Reconstruction — LCMV beamformer to project sensor data to source space.

6. Parcellation — Reduce voxel data to parcel time courses.

The processing code in osl-dynamics is built on top of MNE-Python:

Gramfort, A., Luessi, M., Larson, E., Engemann, D. A., Strohmeier, D., Brodbeck, C., Goj, R., Jas, M., Brooks, T., Parkkonen, L., & Hämäläinen, M. S. (2013). MEG and EEG data analysis with MNE-Python. Frontiers in Neuroscience, 7(267), 1–13. https://doi.org/10.3389/fnins.2013.00267.

Prerequisites

• FSL (needed for surface extraction).

• osl-dynamics (this installs MNE-Python as a dependency). Note, TensorFlow is not required for processing M/EEG. (osl-dynamics can be installed without TensorFlow using the envs/osld.yml environment.)

Input Data

In this tutorial, we will use data from the Wakeman-Henson dataset. This is a visual perception task recorded on an Elekta MEG system with a co-registered structural MRI. The data in BIDS format:

BIDS/
├── sub-01/
│   ├── meg/
│   │   └── sub-01_run-01_task-visual_raw_sss.fif
│   └── anat/
│       └── sub-01_T1w.nii.gz
├── ...

Output is written to derivatives/.

Download the dataset

We will download example data hosted on OSF.

import os

def get_data(name):
    os.system(f"osf -p by2tc fetch data/{name}.zip")
    os.system(f"unzip -o {name}.zip")
    os.remove(f"{name}.zip")
    return f"Data downloaded to: {name}"

# Download the dataset
get_data("wakeman_henson_raw_1_subject")

Setup and Configuration

from pathlib import Path

import mne
import numpy as np
%matplotlib inline

from osl_dynamics import files
from osl_dynamics.meeg import preproc, rhino, source_recon, parcellation
from osl_dynamics.utils.filenames import OSLFilenames

Edit the cell below to match your data.

# Session info
subject = "01"
run = "01"
task = "visual"
id = f"sub-{subject}_run-{run}_task-{task}"

# Paths
bids_dir = Path("BIDS")
raw_file = bids_dir / f"sub-{subject}/meg/sub-{subject}_run-{run}_task-{task}_raw_sss.fif"
mri_file = bids_dir / f"sub-{subject}/anat/sub-{subject}_T1w.nii.gz"
output_dir = Path("derivatives")
plots_dir = Path("plots")

# Preprocessing parameters
resample_freq = 250  # Hz
bandpass = (1, 45)  # Hz
notch_freqs = [50, 100]  # Hz (mains frequency and harmonics)

# Source reconstruction parameters
gridstep = 8  # dipole grid resolution in mm
chantypes = ["mag", "grad"]  # Elekta has both magnetometers and gradiometers
rank = {"meg": 60}  # effective rank after MaxFilter

# Parcellation
parcellation_file = "atlas-Glasser_nparc-52_space-MNI_res-8x8x8.nii.gz"

print(f"Session ID: {id}")
print(f"Raw file: {raw_file}")
print(f"MRI file: {mri_file}")

Step 1: Preprocessing

We clean the sensor-level MEG data by:

• Notch filtering to remove mains frequency contamination.

• Bandpass filtering to remove slow drifts and high-frequency noise.

• Resampling to reduce data size.

• Detecting and annotating bad segments (outlier temporal windows).

• Detecting bad channels (outlier sensors).

Load raw data

raw = mne.io.read_raw_fif(raw_file, preload=True)

Filter and resample

raw = raw.notch_filter(notch_freqs)
raw = raw.filter(
    l_freq=bandpass[0],
    h_freq=bandpass[1],
    method="iir",
    iir_params={"order": 5, "ftype": "butter"},
)
raw = raw.resample(sfreq=resample_freq)

Bad segment detection

We use the Generalized Extreme Studentized Deviate (G-ESD) algorithm to identify windows with outlier standard deviation. We run this separately for magnetometers and gradiometers, in both standard and diff modes (the latter catches segments where the signal changes abruptly within a window).

raw = preproc.detect_bad_segments(raw, picks="mag")
raw = preproc.detect_bad_segments(raw, picks="mag", mode="diff")
raw = preproc.detect_bad_segments(raw, picks="grad")
raw = preproc.detect_bad_segments(raw, picks="grad", mode="diff")

Bad channel detection

Channels with outlier standard deviations are marked as bad.

raw = preproc.detect_bad_channels(raw, picks="mag")
raw = preproc.detect_bad_channels(raw, picks="grad")

ICA artefact rejection

Automatic ICA artefact rejection to remove ECG and EOG artefacts. Two approaches are available:

Option A: ECG/EOG correlation — Uses MNE-Python’s built-in correlation-based detection.

raw, ica, ic_labels = preproc.ica_ecg_eog_correlation(raw, picks="meg")

Option B: MEGNet automatic labelling — Uses mne-icalabel (included in the osl-dynamics conda environments) to classify components with a pre-trained deep learning model. Note, MEGNet was trained on 'mag' sensor topographies but we’ll pass 'meg' to include both magnetometers and gradiometers.

raw, ica, ic_labels = preproc.ica_label(raw, picks="meg", method="megnet")

Save preprocessing QC plots (PSD, sum-of-squares time series, channel standard deviations, ICA component topographies). Check that bad segments, channels, and ICA components are being correctly identified.

preproc.save_qc_plots(raw, plots_dir / id, show=True, ica=ica, ic_labels=ic_labels)

Save preprocessed data

preproc_file = output_dir / "preprocessed" / f"{id}_preproc-raw.fif"
preproc_file.parent.mkdir(parents=True, exist_ok=True)
raw.save(preproc_file, overwrite=True)
print(f"Saved: {preproc_file}")

Step 2: Surface Extraction

We extract three surfaces from the structural MRI using FSL BET:

• Inner skull (inskull) — used for the forward model.

• Outer skull (outskull) — boundary between skull and scalp.

• Outer skin (outskin) — scalp surface (used for coregistration).

• Outer skin (outskin_plus_nose) — scalp surface including the nose, only generated if include_nose=True.

These surfaces define the geometry needed for coregistration and source reconstruction.

Note

No structural MRI? If you don’t have a subject-specific MRI, you can skip this step and use the standard MNI152 brain bundled with osl-dynamics. Set surfaces_dir = files.mni152_surfaces.directory in Step 3 and pass allow_mri_scaling=True during coregistration.

The output plots overlay each extracted surface (yellow line) on the structural MRI. Check that each surface matches the corresponding anatomical boundary. If they don’t, consider using the standard MNI152 brain as a fallback.

surfaces_dir = str(output_dir / "anat_surfaces" / f"sub-{subject}")

rhino.extract_surfaces(
    mri_file=str(mri_file),
    outdir=surfaces_dir,
    include_nose=False,
    do_mri2mniaxes_xform=False,
    show=True,
)

Step 3: Coregistration

Coregistration aligns the MEG sensor coordinate system (“head” space) to the MRI coordinate system. We use the Polhemus fiducials (nasion, LPA, RPA) and headshape points recorded during the MEG session, and fit them to the MRI surfaces using the Iterative Closest Point (ICP) algorithm.

First, let’s create an OSLFilenames container to keep track of all the pipeline output files.

fns = OSLFilenames(
    outdir=str(output_dir / "osl"),
    id=id,
    preproc_file=str(preproc_file),
    surfaces_dir=surfaces_dir,
    # If using standard brain: surfaces_dir=files.mni152_surfaces.directory
)
print(fns)

Extract fiducials and headshape

The Polhemus digitisation points are stored inside the FIF file. We extract them to text files that RHINO uses for coregistration.

rhino.extract_fiducials_and_headshape_from_fif(fns)

Fix stray Polhemus headshape points

Sometimes we want to remove some Polhemus headshape points, we have to do this manually. Note, this step has to be adapted to your specific dataset of interest.

hs = np.loadtxt(fns.coreg.head_headshape_file)
nas = np.loadtxt(fns.coreg.head_nasion_file)
lpa = np.loadtxt(fns.coreg.head_lpa_file)
rpa = np.loadtxt(fns.coreg.head_rpa_file)

remove = np.logical_and(hs[1] > max(lpa[1], rpa[1]), hs[2] < nas[2])
hs = hs[:, ~remove]

print(f"Overwriting: {fns.coreg.head_headshape_file}")
np.savetxt(fns.coreg.head_headshape_file, hs)

Run coregistration

Key parameters:

• use_nose=False — Should we use the nose in the coregistration? (the nose is often removed when an MRI is defaced, in this scenario you shouldn’t use the nose).

• allow_mri_scaling=False — Don’t scale the MRI to fit the headshape. Set True if using the MNI152 standard brain, which needs to be scaled to match the subject’s head size.

rhino.coregister_head_and_mri(
    fns,
    use_nose=False,
    allow_mri_scaling=False,  # set True if using MNI152 standard brain
    show=True,
)

The coregistration plot is saved automatically to fns.coreg_dir/coreg.png. The 3D plot shows the MEG sensors (blue), headshape points (red dots), fiducials, and MRI surfaces. Check that the headshape points sit on the scalp surface and the sensors surround the head correctly. The QC report copies this plot automatically when generated.

If the coregistration looks off, you can try:

• Removing stray headshape points with rhino.remove_stray_headshape_points.

• Manually editing the headshape/fiducial text files in fns.coreg_dir.

• Setting use_nose=False.

Step 4: Forward Model

The forward model (lead field matrix) describes how a dipole at each source location projects onto the MEG sensors. We use a Single Layer (Single Shell) head model based on the inner skull surface and a volumetric dipole grid.

• model="Single Layer" — Single shell head model (standard for MEG).

• gridstep=8 — 8 mm dipole grid spacing. Smaller values give finer resolution but are slower.

rhino.forward_model(fns, model="Single Layer", gridstep=gridstep)

Step 5: Source Reconstruction

We use an LCMV (Linearly Constrained Minimum Variance) beamformer to project the sensor data into source space. The beamformer computes a spatial filter for each voxel that passes activity from that location while suppressing interference from elsewhere.

Key parameters:

• chantypes=["mag", "grad"] — Use both magnetometers and gradiometers (Elekta systems). For CTF, use ["mag"].

• rank={"meg": 60} — Effective rank of the data after MaxFilter. Elekta data is typically rank 60-64.

Compute beamformer weights

source_recon.lcmv_beamformer(fns, raw, chantypes=chantypes, rank=rank)

Apply beamformer

This applies the spatial filters to the sensor data to produce voxel time courses in MNI space. Bad segments are automatically excluded.

voxel_data, voxel_coords = source_recon.apply_lcmv_beamformer(fns, raw)
print(f"Voxel data shape: {voxel_data.shape} (voxels x time)")
print(f"Voxel coords shape: {voxel_coords.shape} (voxels x 3, in MNI mm)")

Step 6: Parcellation

We reduce the high-dimensional voxel data to a smaller number of parcel time courses using a brain atlas. This makes the data more manageable for downstream analysis. See the parcellations page for the full list of available parcellations.

• method="spatial_basis" — Weight voxels by their loading on each parcel (from the atlas) because calculate PCA.

• orthogonalisation="symmetric" — Apply symmetric orthogonalisation to reduce spatial leakage between parcels.

parcel_data = parcellation.parcellate(
    fns,
    voxel_data,
    voxel_coords,
    method="spatial_basis",
    orthogonalisation="symmetric",
    parcellation_file=parcellation_file,
)
print(f"Parcel data shape: {parcel_data.shape} (parcels x time)")

Save parcellated data

We save the parcel time courses as a FIF file. The extra_chans="stim" option preserves any stimulus channels from the original recording.

parc_fif = str(output_dir / "osl" / id / "lcmv-parc-raw.fif")
parcellation.save_as_fif(
    parcel_data,
    raw,
    extra_chans="stim",
    filename=parc_fif,
)
print(f"Saved: {parc_fif}")

A good sanity check is to plot the PSD of each parcel. We expect to see an alpha peak (~10 Hz) that is strongest in posterior parcels. If this is absent or the spectra look unusual, something may have gone wrong in the pipeline. You can also pass power_maps=True to generate band-limited power maps on the cortical surface.

parcellation.save_qc_plots(parc_fif, parcellation_file, show=True)

Summary and Next Steps

We have completed the full MEG preprocessing pipeline:

Step	Name	Output
1	Preprocessing	derivatives/preprocessed/{id}_preproc-raw.fif
2	Surface extraction	derivatives/anat_surfaces/sub-{subject}/
3	Coregistration	derivatives/osl/{id}/coreg/
4	Forward model	derivatives/osl/{id}/coreg/model-fwd.fif
5	Source reconstruction	derivatives/osl/{id}/src/filters-lcmv.h5
6	Parcellation	derivatives/osl/{id}/lcmv-parc-raw.fif

Loading the parcellated data

The final output can be loaded for downstream analysis:

from osl_dynamics.data import Data

data = Data(parc_fif, picks="misc", reject_by_annotation="omit")

Where to go from here

• Canonical HMM analysis — See here for notebooks applying a canonical HMM to the parcellated data.

• Batch processing — See the batch scripts for processing many sessions in parallel.

• Analysis — See the documentation for further tutorials for static and dynamic network analysis.

Notes on CTF and OPM data

This tutorial uses Elekta MEG data. If you are working with different MEG systems (or EEG) see Canonical-HMM-Networks. Note:

• CTF: Use mne.io.read_raw_ctf() to load data. Set chantypes=["mag"] and adjust rank (CTF data is not MaxFiltered, so the rank is typically higher). Fiducials/headshape points may need to be extracted from a .pos file rather than the data file — pass pos_file to OSLFilenames.

• OPM: Loading depends on the OPM system. Adjust chantypes and rank accordingly.