.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/gallery/examples/plot_meg_somato.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_generated_gallery_examples_plot_meg_somato.py: Source localization of somatosensory MEG data ============================================= This example adapts the experiment presented in :footcite:t:`chevalier2020statistical`. We show how to identify which cortical sources are activated during a somatosensory task. To do so, we leverage spatially constrained clustering to effectively reduce the dimensionality of the problem while accounting for the spatial structure of the data. We then perform inference with the desparsified multitask Lasso to perform support recovery from spatio-temporal data. .. admonition:: MultiTaskLasso: The `MultiTaskLasso `_ is a linear model that estimates sparse coefficients for multiple correlated tasks simultaneously. Here, the tasks correspond to the different time points. Instead of fitting a separate Lasso model for each time point, it fits a single model that enforces joint feature selection across all time points. To do so, it uses an L1/L2 mixed-norm regularization, which for a coefficient matrix :math:`W` of shape (n_features, n_tasks) is defined as: .. math:: ||W||_{21} = \sum_{j=1}^p \sqrt{\sum_{t=1}^T W_{j,t}^2} .. GENERATED FROM PYTHON SOURCE LINES 29-33 Load somatosensory MEG data from MNE-Python ------------------------------------------- We use the somatosensory MEG dataset available in MNE-Python. We visualize the evoked response across MEG sensors. .. GENERATED FROM PYTHON SOURCE LINES 33-74 .. code-block:: Python import mne import numpy as np from mne.datasets import somato cond = "somato" data_path = somato.data_path(verbose=True) subject = "01" subjects_dir = data_path / "derivatives/freesurfer/subjects" raw_fname = ( data_path / f"sub-{subject}" / "meg" / f"sub-{subject}_task-{cond}_meg.fif" ) fwd_fname = ( data_path / "derivatives" / f"sub-{subject}" / f"sub-{subject}_task-{cond}-fwd.fif" ) # Read evoked raw = mne.io.read_raw_fif(raw_fname) events = mne.find_events(raw, stim_channel="STI 014") reject = {"grad": 4000e-13, "eog": 350e-6} picks = mne.pick_types(raw.info, meg=True, eeg=True, eog=True) event_id, tmin, tmax = 1, -0.2, 0.25 epochs = mne.Epochs( raw, events, event_id, tmin, tmax, picks=picks, reject=reject, preload=True, ) evoked = epochs.average() evoked = evoked.pick_types("grad") evoked.plot() .. image-sg:: /generated/gallery/examples/images/sphx_glr_plot_meg_somato_001.png :alt: Gradiometers (204 channels) :srcset: /generated/gallery/examples/images/sphx_glr_plot_meg_somato_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Using default location ~/mne_data for somato... Attempting to create new mne-python configuration file: /home/circleci/.mne/mne-python.json Could not read the /home/circleci/.mne/mne-python.json json file during the writing. Assuming it is empty. Got: Expecting value: line 1 column 1 (char 0) Opening raw data file /home/circleci/mne_data/MNE-somato-data/sub-01/meg/sub-01_task-somato_meg.fif... Range : 237600 ... 506999 = 791.189 ... 1688.266 secs Ready. Finding events on: STI 014 111 events found on stim channel STI 014 Event IDs: [1] Not setting metadata 111 matching events found Setting baseline interval to [-0.19979521315838786, 0.0] s Applying baseline correction (mode: mean) 0 projection items activated Loading data for 111 events and 136 original time points ... 0 bad epochs dropped NOTE: pick_types() is a legacy function. New code should use inst.pick(...).
.. GENERATED FROM PYTHON SOURCE LINES 75-83 Preprocessing MEG data for source localization ---------------------------------------------- Before performing source localization, we need to preprocess the MEG data. To do so we rely on the MNE-Python library. We compute the forward model, which describes how the activity of cortical sources is projected onto MEG sensors. We also select the time window of interest, which is the early response. Finally we compute the noise covariance matrix, in order to whiten the data. .. GENERATED FROM PYTHON SOURCE LINES 83-108 .. code-block:: Python from mne.inverse_sparse.mxne_inverse import _prepare_gain # Read forward matrix forward = mne.read_forward_solution(fwd_fname) # Compute noise covariance matrix noise_cov = mne.compute_covariance(epochs, rank="info", tmax=0.0) # We must reduce the whitener since data were preprocessed for removal # of environmental noise with maxwell filter leading to an effective # number of 64 samples. pca = True # Preprocessing MEG data forward, gain, gain_info, whitener, _, _ = _prepare_gain( forward, evoked.info, noise_cov, pca=pca, depth=0.0, loose=0.0, weights=None, weights_min=None, rank=None, ) .. rst-class:: sphx-glr-script-out .. code-block:: none Reading forward solution from /home/circleci/mne_data/MNE-somato-data/derivatives/sub-01/sub-01_task-somato-fwd.fif... Reading a source space... [done] Reading a source space... [done] 2 source spaces read Desired named matrix (kind = 3523 (FIFF_MNE_FORWARD_SOLUTION_GRAD)) not available Read MEG forward solution (8155 sources, 306 channels, free orientations) Source spaces transformed to the forward solution coordinate frame Setting small MEG eigenvalues to zero (without PCA) Reducing data rank from 306 -> 64 Estimating covariance using EMPIRICAL Done. Number of samples used : 6771 [done] Converting forward solution to fixed orientation No patch info available. The standard source space normals will be employed in the rotation to the local surface coordinates.... Changing to fixed-orientation forward solution with surface-based source orientations... [done] Computing inverse operator with 204 channels. 204 out of 306 channels remain after picking Selected 204 channels Whitening the forward solution. Computing rank from covariance with rank=None Using tolerance 2.1e-12 (2.2e-16 eps * 204 dim * 46 max singular value) Estimated rank (grad): 64 GRAD: rank 64 computed from 204 data channels with 0 projectors Setting small GRAD eigenvalues to zero (without PCA) Creating the source covariance matrix Adjusting source covariance matrix. .. GENERATED FROM PYTHON SOURCE LINES 109-110 Selecting relevant time window: focusing on early signal .. GENERATED FROM PYTHON SOURCE LINES 110-122 .. code-block:: Python t_min, t_max = 0.01, 0.05 t_step = 1.0 / 300 # Croping evoked according to relevant time window evoked.crop(tmin=t_min, tmax=t_max) # Choosing frequency and number of clusters used for compression. # Reducing the frequency to 100Hz to make inference faster step = int(t_step * evoked.info["sfreq"]) evoked.decimate(step) y = np.dot(whitener, evoked.data) .. GENERATED FROM PYTHON SOURCE LINES 123-131 Spatially constrained clustering -------------------------------- We then extract the spatial adjacency matrix that is then used in the clustering step to incorporate the spatial structure of the data. For MEG data ``n_clusters = 750`` is generally a good default choice. Taking ``n_clusters > 2000`` might lead to an unpowerful inference. Taking ``n_clusters < 500`` might compress too much the data leading to a compressed problem not close enough to the original problem. .. GENERATED FROM PYTHON SOURCE LINES 131-141 .. code-block:: Python from sklearn.cluster import FeatureAgglomeration # Collecting features' connectivity connectivity = mne.source_estimate.spatial_src_adjacency(forward["src"]) n_clusters = 750 ward = FeatureAgglomeration(n_clusters=n_clusters, connectivity=connectivity) .. rst-class:: sphx-glr-script-out .. code-block:: none -- number of adjacent vertices : 8196 /home/circleci/project/examples/plot_meg_somato.py:135: RuntimeWarning: 0.5% of original source space vertices have been omitted, tri-based adjacency will have holes. Consider using distance-based adjacency or morphing data to all source space vertices. connectivity = mne.source_estimate.spatial_src_adjacency(forward["src"]) .. GENERATED FROM PYTHON SOURCE LINES 142-147 Running clustered inference --------------------------- We can now run the clustered inference using the desparsified multitask Lasso We then select the clusters that are significant at the target FWER level, which is set to 0.1 in this example. .. GENERATED FROM PYTHON SOURCE LINES 147-167 .. code-block:: Python from sklearn.linear_model import MultiTaskLassoCV from hidimstat import CluDL, DesparsifiedLasso # Setting theoretical FWER target fwer_target = 0.1 cludl = CluDL( clustering=ward, desparsified_lasso=DesparsifiedLasso(estimator=MultiTaskLassoCV()), random_state=0, ) cludl.fit_importance(gain, y) selected = cludl.fwer_selection(fwer_target, two_tailed_test=True) # multiplying the -log10(p-value) map by the sign of the estimated coefficients log_pvalues = -np.log10(cludl.pvalues_) * selected .. rst-class:: sphx-glr-script-out .. code-block:: none /home/circleci/project/examples/plot_meg_somato.py:155: FutureWarning: CluDL is deprecated and will be removed in version 0.5.0.Please use class ClusterImportance instead. cludl = CluDL( /home/circleci/project/.venv/lib/python3.13/site-packages/sklearn/cluster/_agglomerative.py:324: UserWarning: the number of connected components of the connectivity matrix is 2 > 1. Completing it to avoid stopping the tree early. connectivity, n_connected_components = _fix_connectivity( Group reid: AR1 cov estimation Using number of clusters for multiple testing correction. .. GENERATED FROM PYTHON SOURCE LINES 168-173 We here used CluDL, which relies on a single clustering. Alternatively, :class:`~hidimstat.EnCluDL` can be used, which relies on aggregating the results over multiple clusterings obtained by running the clustering algorithm on bootstrapped samples of the data. This can lead to more stable results at the cost of a higher computational time. .. GENERATED FROM PYTHON SOURCE LINES 176-181 Visualize the p-value map on the brain surface ---------------------------------------------- We then visualize the identified cortical sources on the brain surface. We plot the -log10(p-value) map to which we assign the sign of the estimated coefficient in the multitask Lasso model. .. GENERATED FROM PYTHON SOURCE LINES 181-214 .. code-block:: Python from pathlib import Path import matplotlib.pyplot as plt from nilearn import datasets, plotting stc_pvals = mne.SourceEstimate( data=log_pvalues[:, np.newaxis], vertices=[forward["src"][0]["vertno"], forward["src"][1]["vertno"]], tmin=0.0, tstep=1.0, subject=f"{subject}", ) # TODO: tmp fix: unlinking fsaverage directory to avoid error fs_dir = Path(subjects_dir) / "fsaverage" if fs_dir.is_symlink(): fs_dir.unlink() mne.datasets.fetch_fsaverage(subjects_dir=subjects_dir) morph = mne.compute_source_morph( src=forward["src"], subject_from=f"{subject}", subject_to="fsaverage", subjects_dir=subjects_dir, spacing=5, ) stc_fsaverage = morph.apply(stc_pvals) start_id = stc_fsaverage.data.shape[0] // 2 rh_stat_map = stc_fsaverage.data[start_id:, 0] .. rst-class:: sphx-glr-script-out .. code-block:: none 0 files missing from root.txt in /home/circleci/mne_data/MNE-somato-data/derivatives/freesurfer/subjects 0 files missing from bem.txt in /home/circleci/mne_data/MNE-somato-data/derivatives/freesurfer/subjects/fsaverage surface source space present ... Computing morph matrix... Left-hemisphere map read. Right-hemisphere map read. 13 smooth iterations done. 12 smooth iterations done. [done] [done] .. GENERATED FROM PYTHON SOURCE LINES 215-245 .. code-block:: Python fsaverage = datasets.fetch_surf_fsaverage("fsaverage5") fig, axes = plt.subplots( 1, 2, subplot_kw={"projection": "3d"}, figsize=(8, 5), gridspec_kw={"hspace": -0.1}, ) # sphinx_gallery_thumbnail_number = 2 for ax, view in zip(axes, ["lateral", "medial"], strict=True): plotting.plot_surf_stat_map( surf_mesh=fsaverage.infl_right, stat_map=rh_stat_map, bg_map=fsaverage.sulc_right, hemi="right", view=view, colorbar=True, cmap="RdBu_r", threshold=-np.log10(fwer_target), vmax=10, axes=ax, title=f"rh-{view} view", ) ax._colorbars[0].set_ylabel("-log10(p-value)") plotting.show() .. image-sg:: /generated/gallery/examples/images/sphx_glr_plot_meg_somato_002.png :alt: rh-lateral view, rh-medial view :srcset: /generated/gallery/examples/images/sphx_glr_plot_meg_somato_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/circleci/project/examples/plot_meg_somato.py:243: UserWarning: You are using the 'agg' matplotlib backend that is non-interactive. No figure will be plotted when calling matplotlib.pyplot.show() or nilearn.plotting.show(). You can fix this by installing a different backend: for example via pip install PyQt6 plotting.show() .. GENERATED FROM PYTHON SOURCE LINES 246-248 Most of the identified sources are located in the somatosensory cortex, which is expected for this somatosensory task. .. GENERATED FROM PYTHON SOURCE LINES 250-253 References ---------- .. footbibliography:: .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 23.123 seconds) **Estimated memory usage:** 3410 MB .. _sphx_glr_download_generated_gallery_examples_plot_meg_somato.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_meg_somato.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_meg_somato.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_meg_somato.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_