{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n# Preprocessing optically pumped magnetometer (OPM) MEG data\n\nThis tutorial covers preprocessing steps that are specific to :term:`OPM`\nMEG data. OPMs use a different sensing technology than traditional\n:term:`SQUID` MEG systems, which leads to several important differences for\nanalysis:\n\n- They are sensitive to :term:`DC` magnetic fields\n- Sensor layouts can vary by participant and recording session due to flexible\n sensor placement\n- Devices are typically not fixed in place, so the position of the sensors\n relative to the room (and through the DC fields) can change over time\n\nWe will cover some of these considerations here by processing the\n`UCL OPM auditory dataset `\n:footcite:`SeymourEtAl2022`\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Authors: The MNE-Python contributors.\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\nimport nibabel as nib\nimport numpy as np\n\nimport mne\n\nsubject = \"sub-002\"\ndata_path = mne.datasets.ucl_opm_auditory.data_path()\nopm_file = (\n data_path / subject / \"ses-001\" / \"meg\" / \"sub-002_ses-001_task-aef_run-001_meg.bin\"\n)\nsubjects_dir = data_path / \"derivatives\" / \"freesurfer\" / \"subjects\"\n\n# For now we are going to assume the device and head coordinate frames are\n# identical (even though this is incorrect), so we pass verbose='error'\nraw = mne.io.read_raw_fil(opm_file, verbose=\"error\")\nraw.crop(120, 210).load_data() # crop for speed" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Examining raw data\n\nFirst, let's look at the raw data, noting that there are large fluctuations\nin the sub 1 Hz band. In some cases the range of fields a single channel\nreports is as much as 600 pT across this experiment.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "picks = mne.pick_types(raw.info, meg=True)\n\namp_scale = 1e12 # T->pT\nstop = len(raw.times) - 300\nstep = 300\ndata_ds, time_ds = raw[picks[::5], :stop]\ndata_ds, time_ds = data_ds[:, ::step] * amp_scale, time_ds[::step]\n\nfig, ax = plt.subplots(layout=\"constrained\")\nplot_kwargs = dict(lw=1, alpha=0.5)\nax.plot(time_ds, data_ds.T - np.mean(data_ds, axis=1), **plot_kwargs)\nax.grid(True)\nset_kwargs = dict(\n ylim=(-500, 500), xlim=time_ds[[0, -1]], xlabel=\"Time (s)\", ylabel=\"Amplitude (pT)\"\n)\nax.set(title=\"No preprocessing\", **set_kwargs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Denoising: Regressing via reference sensors\n\nThe simplest method for reducing low frequency drift in the data is to\nuse a set of reference sensors away from the scalp, which only sample the\nambient fields in the room. An advantage of this method is that no prior\nknowldge of the locations of the sensors is required. However, it assumes\nthat the reference sensors experience the same interference as scalp\nrecordings.\n\nTo do this in our current dataset, we require a bit of housekeeping.\nThere are a set of channels beginning with the name \"Flux\" which do not\ncontain any evironmental data, these need to be set to as bad channels.\nAnother channel -- G2-17-TAN -- will also be set to bad.\n\nFor now we are only interested in removing artefacts seen below 5 Hz, so we\ninitially low-pass filter the good reference channels in this dataset prior\nto regression\n\nLooking at the processed data, we see there has been a large reduction in the\nlow frequency drift, but there are still periods where the drift has not been\nentirely removed. The likely cause of this is that the spatial profile of the\ninterference is dynamic, so performing a single regression over the entire\nexperiment is not the most effective approach.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# set flux channels to bad\nbad_picks = mne.pick_channels_regexp(raw.ch_names, regexp=\"Flux.\")\nraw.info[\"bads\"].extend([raw.ch_names[ii] for ii in bad_picks])\nraw.info[\"bads\"].extend([\"G2-17-TAN\"])\n\n# compute the PSD for later using 1 Hz resolution\npsd_kwargs = dict(fmax=20, n_fft=int(round(raw.info[\"sfreq\"])))\npsd_pre = raw.compute_psd(**psd_kwargs)\n\n# filter and regress\nraw.filter(None, 5, picks=\"ref_meg\")\nregress = mne.preprocessing.EOGRegression(picks, picks_artifact=\"ref_meg\")\nregress.fit(raw)\nregress.apply(raw, copy=False)\n\n# plot\ndata_ds, _ = raw[picks[::5], :stop]\ndata_ds = data_ds[:, ::step] * amp_scale\n\nfig, ax = plt.subplots(layout=\"constrained\")\nax.plot(time_ds, data_ds.T - np.mean(data_ds, axis=1), **plot_kwargs)\nax.grid(True, ls=\":\")\nax.set(title=\"After reference regression\", **set_kwargs)\n\n# compute the psd of the regressed data\npsd_post_reg = raw.compute_psd(**psd_kwargs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Denoising: Regressing via homogeneous field correction\n\nRegression of a reference channel is a start, but in this instance assumes\nthe relatiship between the references and a given sensor on the head as\nconstant. However this becomes less accurate when the reference is not moving\nbut the subject is. An alternative method, Homogeneous Field Correction (HFC)\nonly requires that the sensors on the helmet stationary relative to each\nother. Which in a well-designed rigid helmet is the case.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# include gradients by setting order to 2, set to 1 for homgenous components\nprojs = mne.preprocessing.compute_proj_hfc(raw.info, order=2)\nraw.add_proj(projs).apply_proj(verbose=\"error\")\n\n# plot\ndata_ds, _ = raw[picks[::5], :stop]\ndata_ds = data_ds[:, ::step] * amp_scale\n\nfig, ax = plt.subplots(layout=\"constrained\")\nax.plot(time_ds, data_ds.T - np.mean(data_ds, axis=1), **plot_kwargs)\nax.grid(True, ls=\":\")\nax.set(title=\"After HFC\", **set_kwargs)\n\n# compute the psd of the regressed data\npsd_post_hfc = raw.compute_psd(**psd_kwargs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Comparing denoising methods\n\nDiffering denoising methods will have differing levels of performance across\ndifferent parts of the spectrum. One way to evaluate the performance of a\ndenoising step is to calculate the power spectrum of the dataset before and\nafter processing. We will use metric called the shielding factor to summarise\nthe values. Positive shielding factors indicate a reduction in power, whilst\nnegative means in increase.\n\nWe see that reference regression does a good job in reducing low frequency\ndrift up to ~2 Hz, with 20 dB of shielding. But rapidly drops off due to\nlow pass filtering the reference signal at 5 Hz. We also can see that this\nmethod is also introducing additional interference at 3 Hz.\n\nHFC improves on the low frequency shielding (up to 32 dB). Also this method\nis not frequency-specific so we observe broadband interference reduction.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "shielding = 10 * np.log10(psd_pre[:] / psd_post_reg[:])\n\nfig, ax = plt.subplots(layout=\"constrained\")\nax.plot(psd_post_reg.freqs, shielding.T, **plot_kwargs)\nax.grid(True, ls=\":\")\nax.set(xticks=psd_post_reg.freqs)\nax.set(\n xlim=(0, 20),\n title=\"Reference regression shielding\",\n xlabel=\"Frequency (Hz)\",\n ylabel=\"Shielding (dB)\",\n)\n\n\nshielding = 10 * np.log10(psd_pre[:] / psd_post_hfc[:])\n\nfig, ax = plt.subplots(layout=\"constrained\")\nax.plot(psd_post_hfc.freqs, shielding.T, **plot_kwargs)\nax.grid(True, ls=\":\")\nax.set(xticks=psd_post_hfc.freqs)\nax.set(\n xlim=(0, 20),\n title=\"Reference regression & HFC shielding\",\n xlabel=\"Frequency (Hz)\",\n ylabel=\"Shielding (dB)\",\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Filtering nuisance signals\n\nHaving regressed much of the high-amplitude, low-frequency interference, we\ncan now look to filtering the remnant nuisance signals. The motivation for\nfiltering after regression (rather than before) is to minimise any filter\nartefacts generated when removing such high-amplitude interfece (compared\nto the neural signals we are interested in).\n\nWe are going to remove the 50 Hz mains signal with a notch filter,\nfollowed by a bandpass filter between 2 and 40 Hz. From here it becomes clear\nthat the variance in our signal has been reduced from 100s of pT to 10s of\npT instead.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# notch\nraw.notch_filter(np.arange(50, 251, 50), notch_widths=4)\n# bandpass\nraw.filter(2, 40, picks=\"meg\")\n# plot\ndata_ds, _ = raw[picks[::5], :stop]\ndata_ds = data_ds[:, ::step] * amp_scale\nfig, ax = plt.subplots(layout=\"constrained\")\nplot_kwargs = dict(lw=1, alpha=0.5)\nax.plot(time_ds, data_ds.T - np.mean(data_ds, axis=1), **plot_kwargs)\nax.grid(True)\nset_kwargs = dict(\n ylim=(-500, 500), xlim=time_ds[[0, -1]], xlabel=\"Time (s)\", ylabel=\"Amplitude (pT)\"\n)\nax.set(title=\"After regression, HFC and filtering\", **set_kwargs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generating an evoked response\n\nWith the data preprocessed, it is now possible to see an auditory evoked\nresponse at the sensor level.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "events = mne.find_events(raw, min_duration=0.1)\nepochs = mne.Epochs(\n raw, events, tmin=-0.1, tmax=0.4, baseline=(-0.1, 0.0), verbose=\"error\"\n)\nevoked = epochs.average()\nt_peak = evoked.times[np.argmax(np.std(evoked.copy().pick(\"meg\").data, axis=0))]\nfig = evoked.plot()\nfig.axes[0].axvline(t_peak, color=\"red\", ls=\"--\", lw=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visualizing coregistration\nBy design, the sensors in this dataset are already in the scanner RAS coordinate\nframe. We can thus visualize them in the FreeSurfer MRI coordinate frame by computing\nthe transformation between the FreeSurfer MRI coordinate frame and scanner RAS:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mri = nib.load(subjects_dir / \"sub-002\" / \"mri\" / \"T1.mgz\")\ntrans = mri.header.get_vox2ras_tkr() @ np.linalg.inv(mri.affine)\ntrans[:3, 3] /= 1000.0 # nibabel uses mm, MNE uses m\ntrans = mne.transforms.Transform(\"head\", \"mri\", trans)\n\nbem = subjects_dir / subject / \"bem\" / f\"{subject}-5120-bem-sol.fif\"\nsrc = subjects_dir / subject / \"bem\" / f\"{subject}-oct-6-src.fif\"\nmne.viz.plot_alignment(\n evoked.info,\n subjects_dir=subjects_dir,\n subject=subject,\n trans=trans,\n surfaces={\"head\": 0.1, \"inner_skull\": 0.2, \"white\": 1.0},\n meg=[\"helmet\", \"sensors\"],\n verbose=\"error\",\n bem=bem,\n src=src,\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting the inverse\nNow we can compute a forward and inverse:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fwd = mne.make_forward_solution(\n evoked.info,\n trans=trans,\n bem=bem,\n src=src,\n verbose=True,\n)\nnoise_cov = mne.compute_covariance(epochs, tmax=0)\ninv = mne.minimum_norm.make_inverse_operator(evoked.info, fwd, noise_cov, verbose=True)\nstc = mne.minimum_norm.apply_inverse(\n evoked, inv, 1.0 / 9.0, method=\"dSPM\", verbose=True\n)\nbrain = stc.plot(\n hemi=\"split\",\n size=(800, 400),\n initial_time=t_peak,\n subjects_dir=subjects_dir,\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n.. footbibliography::\n\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.2" } }, "nbformat": 4, "nbformat_minor": 0 }