Adept level

Welcome to the adept level! You have heard our previous speakers talk about continuous MEG data, epochs, evokeds, inverse solutions, etc. In this hands-on session, we're going to use MNE-Python to analyze some data to see how this all works in practice.

I'm going to assume you have some programming experience (for example, by completing the beginners level), but are not necessarily familiar with the Python programming language. As you progress through the session, you will learn both the basics of Python and the basic of data analysis side by side.

The experiment

I have prepared a bit of MEG data for you. It was recorded while a brave volunteer was in the scanner, listening to auditory beeps and looking at visual checkerboards. The volunteer was staring at a "fixation cross" in the center of a screen in front of him/her. From time to time, beeps sounded from either the left or the right side of his/her head. At other times, a checkerboard would be displayed either to the left or right side of to the cross. The volunteer was instructed to keep his/her eyes fixed on the cross and not directly look at the checkerboards. This way, the checkerboard was presented either to the left or right visual field of the volunteer.

By analyzing the MEG signal, we should be able to see the activity in the auditory and visual cortices, perhaps even see a difference whether the stimulus was presented on the left or on the right side.

Loading the data, or, how to import and call a function in Python

The MEG data is stored on disk as it was produced by the Vectorview MEG scanner, namely in its "FIF" file format. We refer to this as "raw" data, as we haven't done anything to it yet.

Loading raw MEG data into memory is easy if you know how: just call the read_raw_fif function and store the result in some variable. There's one catch though. Since there are tons of functions available to you, to keep track of them all, they are organized in different "modules". You can think of a module as a toolbox that has a nice label on it so we know what is inside.

The read_raw_fif function we want to use is kept inside the io module that itself is kept inside the mne module, giving it the name: mne.io (dots indicate "inside of"). In order for us to use it, we must first "import" the function from the mne.io module. In our metaphor, importing a function is like taking a tool from a toolbox and placing it on our workbench. Here is the line of Python code that accomplishes this:

from mne.io import read_raw_fif

From that point on, the read_raw_fif function is available for you to use.

Below is the code to read in the MEG data. The box it is in is called a "cell". To run it, click in the cell and press Ctrl+Enter (both keys at the same time) or click the "Run" button in the toolbar at the top of the this webpage.


In [ ]:
from mne.io import read_raw_fif
raw = read_raw_fif('data/sample-raw.fif', preload=True)
print(raw)

Dissecting the program above

The first line of the program above imports the read_raw_fif function from the mne.io module.

The second line of the program is the most complicated. A lot of stuff is going on there:

The read_raw_fif function is called with two parameters. The first parameter is a piece of text (a "string") containing the name of the FIF file to load. Literal text (strings) must always be enclosed in ' quotes. The second parameter is a "named" parameter, which we will use a lot during this session (see below). This parameter is set to the special value True. Python has three special values: True, False, and None, which are often used to indicate "yes", "no", and "I don't know/care" respectively. Finally, the result is stored in a variable called raw.

The last line of the program calls the print function, which is used to display things. Here, it is called with the raw variable as parameter, so it displays the data contained in this variable, namely the data we loaded with read_raw_fif.

Named parameters

Many functions of MNE-Python take dozens of parameters that fine-tune exactly how to perform some operation. If you had to specify them all every time you want to call a function, you'd spend ages worrying about little details and get nothing done. Luckily, Python allows us to specify default values for parameters, which means these parameter may be omitted when calling a function to use the default. In MNE-Python, most parameters have a default value, so while a function may have 20 parameters, you only have to specify one or two. The rest of the parameters are like little knobs and buttons you can use to fine tune things, or just leave alone. This allows MNE-Python to keep simple things simple, while making complicated things possible.

Parameters with default values are called "named" parameters, and you specify them with name=value. The preload parameter that you saw in the program above is such a named parameter. It controls whether to load all of the data in memory, or only read the "metadata" of the file, i.e., when it was recorded, how long it is, how many sensors the MEG machine had, etc. By default, preload=False, meaning only the metadata is read. In the example above, we set it to True, indicating we wish to really load all the data in memory.

Visualizing the data

Raw data can be visualized (or "plotted") by the plot_raw function that is kept inside the mne.viz module. It needs one parameter: the variable containing the data you wish to plot. (It also has a lot of named parameters, but you can leave them alone for now.)

I'm going to let you write the visualization code, but first, there is a little housekeeping that we need to do. We need to tell the visualization engine to send the results to your browser and not attempt to open a window on the server where this code is running. Please run the cell below:


In [ ]:
%matplotlib notebook
print('From now on, all graphics will send to your browser.')

Now, it's your turn! Write the Python code that will visualize the raw MEG data we just loaded. Keep the following things in mind:

  1. The function is called plot_raw and is kept inside the mne.viz module. Remember to import to function first!
  2. Call the function with one parameter, namely the raw variable we created above that contains the MEG data.
  3. Assign the result of the plot_raw function to a variable (pick any name you want), otherwise the figure will show twice.

Use the cell below to write your code:


In [ ]:

If you wrote the code correctly, you should be looking at a little interface that shows the data collected on all the MEG sensors. Click inside the scrollbars or use the arrow keys to explore the data.

Browsing through the channels, you will notice there are several channel types:

  • 204 MEG gradiometers (102 pairs of two)
  • 102 MEG magnetometers
  • 60 EEG sensors
  • 9 STIM channels
  • 1 EOG sensor

All these channels record different information about the volunteer and the environment. We will first take a look at the STIM channels.

Events, or, how to read the documentation

In the figure you just made, scroll down and take a look at channel STI 014. On this channel, the computer that is presenting the stimuli was sending timing information to the MEG equipment. Whenever a stimulus (checkerboard or beep) was presented, the signal at this channel jumps briefly from 0 to either 1, 2, 3 or 4, indicating the type of stimulus.

We can use this channel to create an "events" matrix: a table listing all the times a stimulus was presented, along with the time of the event and the type of stimulus. The function to do this is called find_events, and is kept inside the mne module.

In this document, all the function names are links to their documentation. Click on find_events to pull up its documentation. It will open a new browser tab. It should look like this:

Looking at the function "signature" reveals that many of the parameters have default values association with them. This means these are named parameters and we can ignore them if we want. There is only a single required parameter, named raw. Looking at the parameter list, it seems we need to set it to the raw data we just loaded with the read_raw_fif function. If we called the function correctly, it should provide us with an "array" (don't worry about what an array is for now) with all the events.

Now, call the function and find some events! Keep the following things in mind:

  1. The function is called find_events and is kept inside the mne module. Remember to import to function first!
  2. Call the function. Use the documentation to find out what parameters it needs.
  3. Assign the result to a variable called events.

In [ ]:

If you called the function correctly, running the cell below should display the found events on top of the raw data. It should show as cyan lines, with a number on top indicating the type of event. These numbers are referred to as "event codes".


In [ ]:
fig = plot_raw(raw, events=events)

Epochs, or, how to create a dictionary

Now that we have the information on what stimulus was presented at one time, we can extract "epochs". Epochs are little snippets of signal surrounding an event. These epochs can then be averaged to produce the "evoked" signal.

In order to create epochs, we need a way to translate the event codes (1, 2, 3 and 4) to something more descriptive. This can be done using a new type of variable called a "dictionary". A Python dictionary, allows us (or rather the computer) to "look things up". The following piece of code creates a dictionary called event_id. Take a look and run it:


In [ ]:
event_id = {'audio/left': 1, 'audio/right': 2, 'visual/left': 3, 'visual/right': 4}

A dictionary is created by using curly braces { } and colons :. I've spread out the code over multiple lines to make things a little clearer. The way you create a dictionary is to say {this: means that} and you use commas if you want to put more than one thing in the dictionary.

Armed with our event_id dictionary, we can move on to creating epochs.

For each event, let's cut a snippet of signal from 0.2 seconds before the moment the stimulus was presented, up until 0.5 seconds after it was presented. If we take the moment the stimulus was presented as time 0, we will cut epochs from -0.2 until 0.5 seconds.

The function to do this is called Epochs (with a capital E). Click on the function name to open its documentation and look at the parameters it needs.

That's a lot of parameters!
Ok, how many are optional?
Almost all of them, phew!

Here are the points of the documentation that are most relevant to us right now:

  • There are two required arguments, the raw data and the events array we created earlier.
  • The next optional parameter is the event_id dictionary we just created.
  • The next two optional parameters, tmin and tmax, specify the time range to cut epochs for. They are set to -0.2 to 0.5 seconds by default. Neat! That's exactly what we wanted. We can leave this two alone then.
  • We can also leave the rest of the parameters alone.

Go ahead and import the Epochs function from the mne module. Then call it with the correct parameters and store the result in a variable called epochs (small e):


In [ ]:

MNE-Python has plotting functions for almost everything. The cell below will create an interactive visualization of your epochs. Click inside the scrollbars or use the arrow keys to scroll through the data.


In [ ]:
from mne.viz import plot_epochs
fig = plot_epochs(epochs)

Visualizing the evoked field, or, working with objects

Throughout this exercise, we have created many variables: raw, events, event_id, and epochs. Up to now, we've treated these as simple boxes that hold some data, or rather "objects" as they are called in Python. However, the box/object metaphor is not really a good one. Variables are more like little machines of their own. They can do things!

Take for example our event_id dictionary. We can look things up in it. You do this using square brackets [ ]. (Python uses many different kinds of brackets!) For example, to look up the event code for the event 'visual/left', you would write:

event_id['visual/left']

Try it! Get the event code for the event 'audio/right' (should be 2):


In [ ]:

Both the raw and the epochs variables are very powerful objects. If you want, you can look at the documentation for Raw and Epochs to see the long list of things they can do.

A very useful thing that epochs does is behave like a dictionary: you can look up stuff in it. More precisely, we can select all the epochs corresponding to a certain type of event by using the square backets [ ]. For example, to select all the epochs in which a checkerboard was shown in the left visual field of the volunteer, we would write:

epochs['visual/left']

Another useful thing is that both raw and epochs objects know how to visualize (i.e. "plot") themselves. You already know that modules hold functions, but objects can hold functions too. Functions that are kept inside objects are called "methods", to distinguish them from "functions" that are kept inside modules.

Instead of using the plot_raw and plot_epochs functions, we can use the plotting methods of the objects, like this:

fig = raw.plot()
fig = epochs.plot()

Notice how the .plot() method calls don't need any parameters: they already know which object they need to visualize, namely the object they belong to. In MNE-Python, many objects have such a plot method.

We can string together the "lookup" functionality and the visualization functionality, like this:

fig = epochs['visual/left'].plot()

Go ahead and try it: visualize the epochs belonging to the different stimuli. Remember, we called them: 'audio/left', 'audio/right', 'visual/left' and 'visual/right'.


In [ ]:


In [ ]:


In [ ]:


In [ ]:

Let's move on from epochs to evoked fields.

As you have heard in the talks earlier today, evoked fields can be analyzed by taking all the epochs that correspond to a certain experimental condition and computing their average. The epochs object provides a method named average to compute, you guessed it, the average of all the epochs. This method produces an Evoked object that represents the evoked field. Like most objects, it has a plot method of its own to visualize it. For example, here's how to compute the evoked field in response to the checkerboards presented to the left visual field:

evoked_visual_left = epochs['visual/left'].average()

Now, everything you've learned so far needs to come together. You are given four code cells below. In each code cell, plot the evoked field corresponding to one of the four event types, keeping the following in mind:

  • Select the epochs that correspond to the event type using the square brackets [ ].
  • Either:
    • store the result in some variable and call the average method of the new variable, or
    • string the square brackets together with the average method like in the example.
  • Assign the result of the average method, i.e. the evoked field, to a new variable (pick a good name, for example evoked_visual_left)
  • Use the plot method of the new variable containing the evoked field to visualize it. Store the result in a variable (e.g. fig) or the figure will show up twice.

In [ ]:


In [ ]:


In [ ]:


In [ ]:

In the visualizations you just made, all the channels are drawn on top of each other. (This is called a "butterfly" plot.) These visualizations are great to determine when things are happening. Around 0.1 seconds after an auditory beep was played, we get a burst of activity. The same is true for the checkerboards, with an additional peak at around 0.18 seconds.

Making topographic maps

Now let's look at where these bursts of activity originate from. However, we're not going to go all the way and perform source estimation to trace the activity back to the cortex of the brain. You can learn how to do that in the expert level. For now, we can get some rough idea by looking at which sensors are picking up most of the activity: are they at the front or the back of the head? This question can be nicely answered by topographic maps: maps of the sensors around the head.

Making topographic maps is really easy. Objects of the type Evoked, like the ones you have been making to create the butterfly plots, have a method called plot_topomap. For example, try running the cell below:


In [ ]:
evoked_visual_left = epochs['visual/left'].average()
fig = evoked_visual_left.plot_topomap()

By default, plot_topomap draws topographic maps for 4 different points in time. It shows the head of the subject, along with the positions of all the sensors of the MEG scanner. Colors indicate the strength of the magnetic field and whether the field is exiting the head (red) or entering the head (blue). Remember from the talks you've heared today that the magnetic field is generated by a current inside the brain:

If you take a look at the documention of the plot_topomap method, you will see that the method takes some parameters. All of them are optional, but one that is of interest to us is the first one, called times. This parameters lets us specify at which times we want a topographic map to be drawn. Intriguingly, we can set it to the string 'interactive' to create an interactive plot where we can use a slider to change the time of the map.

Try plotting a topographic map, with the times parameter set to 'interactive':


In [ ]:

Now you can plot topographic maps for all conditions and analyze the sensor data. See if you can spot how the checkerboard activate the visual cortex at the back of the head, while the auditory beeps activate the auditory cortices at the left and right side of the head.

Continue onward to the expert level!

Great job making it this far. You've taken some steps along the path towards becoming a data analyst. Take a look at the clock. If we still have some time, I'd like to invite you to move on to the next level.

Move on to the expert level