08-source_estimation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Source estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib qt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pathlib\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import mne_bids\n",
    "import mne\n",
    "\n",
    "matplotlib.use('Qt5Agg')\n",
    "mne.set_log_level('warning')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Start with some fresh epochs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bids_root = pathlib.Path('out_data/sample_BIDS')\n",
    "\n",
    "bids_path = mne_bids.BIDSPath(subject='01',\n",
    "                              session='01',\n",
    "                              task='audiovisual',\n",
    "                              run='01',\n",
    "                              datatype='meg',\n",
    "                              root=bids_root)\n",
    "\n",
    "raw = mne_bids.read_raw_bids(bids_path)\n",
    "raw.load_data()\n",
    "raw.filter(l_freq=0.1, h_freq=40)\n",
    "events, event_id = mne.events_from_annotations(raw)\n",
    "\n",
    "tmin = -0.200\n",
    "tmax = 0.500\n",
    "baseline = (None, 0)\n",
    "\n",
    "epochs = mne.Epochs(raw,\n",
    "                    events=events,\n",
    "                    event_id=event_id,\n",
    "                    tmin=tmin,\n",
    "                    tmax=tmax,\n",
    "                    baseline=baseline,\n",
    "                    preload=True,\n",
    "                    proj=False)\n",
    "epochs.save(pathlib.Path('out_data') / 'epochs_for_source_epo.fif')\n",
    "epochs.info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## View the BEM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subjects_dir = pathlib.Path(mne.datasets.sample.data_path()) / 'subjects'\n",
    "\n",
    "mne.viz.plot_bem(subject='sample', subjects_dir=subjects_dir,\n",
    "                 orientation='coronal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coregistration!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs_fname = pathlib.Path('out_data') / 'epochs_for_source_epo.fif'\n",
    "\n",
    "mne.gui.coregistration(subject='sample', subjects_dir=subjects_dir,\n",
    "                       inst=epochs_fname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trans_fname = pathlib.Path('out_data') / 'sample.trans'\n",
    "info = mne.io.read_info(epochs_fname)\n",
    "fig = mne.viz.plot_alignment(info=info, trans=trans_fname, subject='sample', dig=True,\n",
    "                             subjects_dir=subjects_dir, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute the source space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subject = 'sample'\n",
    "src = mne.setup_source_space(subject=subject,\n",
    "                             spacing='oct4',  # Use oct6 during an actual analysis!\n",
    "                             subjects_dir=subjects_dir,\n",
    "                             add_dist=False)  # Remove this one during an actual analysis!\n",
    "src"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mne.viz.plot_alignment(info=info, trans=trans_fname, subject=subject,\n",
    "                       src=src, subjects_dir=subjects_dir, dig=True,\n",
    "                       surfaces=['head-dense', 'white'], coord_frame='meg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute the forward solution!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conductivity = (0.3,)  # for single layer – used in MEG\n",
    "# conductivity = (0.3, 0.006, 0.3)  # for three layers – used in EEG\n",
    "model = mne.make_bem_model(subject=subject, ico=4,\n",
    "                           conductivity=conductivity,\n",
    "                           subjects_dir=subjects_dir)\n",
    "model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bem_sol = mne.make_bem_solution(model)\n",
    "bem_sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bem_fname = pathlib.Path('out_data') / 'sample_bem.fif'\n",
    "mne.bem.write_bem_solution(bem_fname, bem_sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fwd = mne.make_forward_solution(raw_fname,\n",
    "                                trans=trans_fname,\n",
    "                                src=src,\n",
    "                                bem=bem,\n",
    "                                meg=True, # include MEG channels\n",
    "                                eeg=False, # exclude EEG channels\n",
    "                                mindist=5.0, # ignore sources <= 5mm from inner skull\n",
    "                                n_jobs=1) # number of jobs to run in parallel\n",
    "fwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fwd_fname = pathlib.Path('out_data') / 'sample_fwd.fif'\n",
    "mne.write_forward_solution(fwd_fname, fwd, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute noise covariance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "noise_cov = mne.compute_covariance(epochs, tmax=0.,\n",
    "                                   method=['shrunk', 'empirical'],\n",
    "                                   rank='info')\n",
    "mne.viz.plot_cov(noise_cov, info=info)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "contrast.plot_white(noise_cov)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create the inverse operator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mne.forward import read_forward_solution\n",
    "from mne.minimum_norm import (make_inverse_operator, apply_inverse,\n",
    "                              write_inverse_operator)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fwd_fname = os.path.join(data_path,\n",
    "    'derivatives/meg_derivatives/sub-01/ses-meg/meg/sub-01-meg-fwd.fif')\n",
    "fwd = mne.read_forward_solution(fwd_fname)\n",
    "fwd = mne.convert_forward_solution(fwd, surf_ori=True)\n",
    "\n",
    "# Restrict forward solution as necessary for MEG\n",
    "fwd = mne.pick_types_forward(fwd, meg=True, eeg=False)\n",
    "\n",
    "# make an M/EEG, MEG-only, and EEG-only inverse operator\n",
    "info = contrast.info\n",
    "inverse_operator = make_inverse_operator(info, fwd, noise_cov,\n",
    "                                         loose=0.2, depth=0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Apply the inverse operator -> Calculate the source estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "method = \"dSPM\"\n",
    "snr = 3.\n",
    "lambda2 = 1. / snr ** 2\n",
    "stc = apply_inverse(epochs['Auditory/Left', inverse_operator, lambda2,\n",
    "                    method=method, pick_ori=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Show the result!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "brain = stc.plot(surface='inflated',\n",
    "                 hemi='both',\n",
    "                 subjects_dir=subjects_dir,\n",
    "                 time_viewer=True)"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}