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.
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.
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)
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
.
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.
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:
plot_raw
and is kept inside the mne.viz
module. Remember to import
to function first!raw
variable we created above that contains the MEG data.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:
All these channels record different information about the volunteer and the environment. We will first take a look at the STIM channels.
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:
find_events
and is kept inside the mne
module. Remember to import
to function first!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)
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.
Here are the points of the documentation that are most relevant to us right now:
raw
data and the events
array we created earlier.event_id
dictionary we just created.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.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)
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:
[ ]
.average
method, i.e. the evoked field, to a new variable (pick a good name, for example evoked_visual_left
)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.
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.
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.