First-level analyses (using FSL)#
This notebook provides the first part (first-level analysis) of an introduction to multilevel models and how to implement them in FSL. The next notebook, run_level_analysis.ipynb, will continue with run-level analysis, and next week we will conclude this topic with group-level analyses. If you haven’t done so already, please go throught linux_and_the_command_line notebook first!
You actually know most about “first-level analyses” already, as it describes the process of modelling single-subject (single-run) timeseries data using the GLM, which we discussed at length in week 2. In this notebook, you’ll learn how to use FSL for (preprocessing and) first-level analyses. Note that there are several other excellent neuroimaging software packages, such as SPM (Matlab-based), AFNI, BROCCOLI, Freesurfer, and Nilearn. In fact, in week 7, there is a notebook on how to use Nilearn to fit statistical models on fMRI data.
We like FSL in particular, as it’s a mature and well-maintained package, free and open-source, and provides both a graphical interface and a command line interface for most of its tools.
What you’ll learn: after this lab, you’ll be able to …
visualize and inspect (f)MRI in FSLeyes;
set up a first-level model in FSL FEAT
Estimated time needed to complete: 2-4 hours
# Import functionality
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Using FSLeyes#
In the previous notebook about Linux and the command line, we discussed how to use built-in Linux commands and FSL-specific commands using the command line. Now, let’s focus on using a specific FSL tool called “FSLeyes”, which is a useful program to view (f)MRI files (basically any nifti-file). We can start this graphical tool from the command line with the command: fsleyes &. We append the & to the FSLeyes command because this makes sure FSLeyes will run “in the background” and thus allows us to keep using the command line (if we would not use the &, FSLeyes would run be active “in your terminal”, and we wouldn’t be able to use that anymore).
This should open a new window with FSLeyes. Now, let’s open a file!
Now, FSLeyes should look similar to the figure below.
As you can see, FSLeyes consists of multiple views, panels, and toolbars. You should see an anatomical MRI image shown in three perspectives (sagittal, coronal, and axial) with green crosshairs in the “Ortho view”.
Notice the labels on the sides of the image: P (posterior), A (anterior), I (inferior), S (superior), L (left), and R (right).
In the “Overlay toolbar”, you can adjust both the brightness and contrast of the image. One could argue that the default settings are a little bit too dark (for this scan at least).
Importantly, while moving the location of the crosshairs, you see in the lower right corner (under the “Location panel”) that some numbers change. The left side of the panel shows you both the location of the crosshairs in “world coordinates” (relative to the isocenter of the scanner, in millimeters) and in “voxel coordinates” (starting from 0).
The right side of the panel shows you the voxel coordinates of the crosshairs together with the signal intensity, e.g., \([156\ 123\ 87]: 600774.75\). Try changing the crosshairs from a voxel in the CSF, to a voxel in gray matter, and a voxel in white matter and observe how the signal intensity changes.
# Replace the Nones with the correct numbers
X = None
Y = None
Z = None
# YOUR CODE HERE
raise NotImplementedError()
''' Tests the above ToDo.'''
from niedu.tests.nii.week_5 import test_how_many_slices_fsleyes
test_how_many_slices_fsleyes(X, Y, Z)
Now, let’s check out some functional data (fMRI) in FSLeyes!
FSLeyes neatly plots the functional data on top of the anatomical image. Note that the resolution differs between the anatomical (\(1 \times 1 \times 1\) mm) and functional (\(2.7 \times 2.7 \times 2.97\) mm)! FSLeyes still allows you to view these images on top of each other by interpolating the images.
Now, let’s make our view a little bit nicer.
Setting a minimum value will threshold the image at this value, which will mask most of the background. Now, let’s also change the color scheme for the func file to be able to distinguish the func and anat file a little better:
Now you can clearly see that setting a threshold did not mask all non-brain tissue, as some skull seems to still show “activity”!
Again, just like with the anatomical scan, the “Location panel” shows you both the (voxel and world) coordinates as well as the signal intensity.
Instead of manually adjusting the volume number, you can press the “Video” icon (next to the big “+” icon, which you can use to toggle the crosshairs on or off) to start “Movie mode” in which FSLeyes will loop through all volumes. Due to the lag on the server (or a software bug in FSLeyes, I’m not sure), however, FSLeyes is not fast enough to render the volumes, so you’ll likely see a black screen…
That said, looping through your volumes is usually a great way to visually check for artifacts (e.g. spikes) and (extreme) motion in your data! We highly recommend doing this for all your data!
Another useful way to view your fMRI data is using the “Time series” view.
FSLeyes has many more features than we discussed in this section (e.g., working with standard atlases/parcellations). We’ll demontrate these features throughout the rest of this lab where appropriate!
Skullstripping#
When viewing the anatomical scan in the previous section, you might have noticed that the scan still contained the subject’s skull. Usually, the skull is removed because it improves registration to the subject’s anatomical T1-scan (especially if you use BBR-based registration). This process is often called “skullstripping” and FSL contains the skullstripping-tool bet (which stands for “Brain Extraction Tool”). We’ll demonstrate how to use bet using its command-line interface.
Like all FSL tools with command-line interfaces, you can get more info on how to use it by calling the command with the flag --help.
You see that the “bet” command should be called as follows:
bet <input> <output> [options/flags]
We recommend visually checking the result of skullstripping for imperfections (e.g. part of the brain that are excluded or parts of non-brain that are included). Let’s do this in FSLeyes!
Then, load the non-skullstripped anatomical file (sub-02_desc-preproc_T1w.nii.gz) first. Then, overlay the skullstripped image (the sub-02_desc-preproc+bet_T1w.nii.gz file you just created). The difference between the non-skullstripped and the skullstriped image is, now, however not clear; to make this clearer, we can pick a different colormap for the skullstripped brain! For example, change the “Grayscale” colormap to “Red-Yellow”. Also, change the opacity to about 50% to better see the underlying non-skullstripped anatomical scan.
Scroll through the image to see whether bet did a good job.
Also, it failed to remove some non-brain tissue. For example, what non-brain structure do you think is located at \([X=143, Y=178, Z=72]\)? And at \([X=140, Y=124, Z=44]\)?
You could, of course, just skullstrip all your participants without looking at the resulting skullstripped brain images. However, as we just saw, this may lead to suboptimal results, which may bias subsequent preprocessing steps, such as functional → anatomical registration and spatial normalization to a standard template. As such, we recommend visually inspecting the skullstrip result for each participant in your dataset!
If you find that bet did not do a proper job, you can adjust several parameters to improve the skullstripping operation. The two most useful parameters to adjust are “robust brain centre estimation” (flag -R), “fractional intensity threshold” (flag -f <value>), and the vertical gradient (flag -g <values>). We recommend always enabling “robust brain centre estimation” (which iterates the operation multiple times, which usually gives more robust results).
The proper value for fractional intensity threshold (-f) and the vertical gradient (-g), however, you need to “tune” yourself (although the default, 0.5, is a good starting value). Smaller values for fractional intensity threshold give you larger brain outline estimates (i.e., “threshold” to be used for image is lower). With the vertical gradient you can tell bet to be relatively more “aggressive” (with values < 0) or “lenient” (with values > 0) with removing tissue at the bottom (inferior side) of the brain. So, e.g., using -g -0.1 will remove relatively more tissue at the bottom of the brain.
To include these parameters, you can call bet with, for example, a “fractional intensity threshold” of 0.4, a vertical gradient of -0.1, and robust brain centre estimation enabled as follows:
bet <input> <output> -f 0.4 -g -0.1 -R
For this ToDo, try out some different parameters for fractional intensity threshold and vertical gradient, and enable robust brain centre estimation. Note that the output file will be overwritten if you use the same output-name (which is fine if you don’t want to save the previous result). Also, there may not be a combination of settings that creates a “perfect” skullstrip result (at least for this brain).
Once you’re happy with the segmentation, overlay it on the original anatomical file in FSLeyes (and use a different colormap). Then, make a snapshot using the camera icon and save it in your week_5 directory (name it “my_skullstripped_image.png”). After you created and saved the snapshot, add the following line to the answer-cell below:
After adding this to the text-cell below and running the cell, the snapshot should be visible in the notebook. (This allows us to check whether you actually did the ToDo; there is not necessarily a “correct” answer to this ToDo — we just want you to mess around with bet and the command line).
YOUR ANSWER HERE
Note that adjusting the parameters of bet may in fact not completely solve all the imperfections in the skullstripping procedure. If you really care about accurate segmentation (e.g., when doing anatomical analyses such as voxel-based morphometry or surface-based analyses), you may consider editing your skullstrip result manually. FSLeyes has an “edit mode” (Tools → Edit mode) in which can manually fill or erase voxels based on your own judgment what is actual brain tissue and what is not.
Skullstripping is actually the only preprocessing operation you have to do manually using the command line (which might be a good idea, anyway, as you should inspect the result manually). The rest of the preprocessing steps within FSL are actually neatly integrated in the FSL tool “FEAT” (FMRI Expert Analysis Tool), which is FSL’s main analysis workhorse (implementing first-level and grouplevel analyses). We’ll discuss doing first-level and fixed-effect analyses with FEAT later in this tutorial, but first let’s talk about the dataset that we’ll use for the upcoming lab!
The “NI-edu” dataset#
So far in this course, we have used primarily simulated data. From this week onwards, we’ll also use some real data from a dataset (which we’ll call “NI-edu”) acquired at our 3T scanner at the University of Amsterdam. The NI-edu study has been specifically “designed” for this (and the follow-up “pattern analysis”) course. The experimental paradigms used for this dataset are meant to investigate how the brain processes information from faces and the representation of their subjective impressions (e.g., their perceived attractiveness, dominance, and trustworthiness). In addition to tasks that investigate face processing specifically, the dataset also includes several runs with a “functional localizer” task, which can be used to identify regions that preferentially respond to a particular object category (like faces, bodies, characters, etc.). In this course, we’ll mostly use the localizer data, because the associated task (design) is relativelly simple and generates a robust signal.
Below, we download the data needed for this tutorial, which is the anatomical data and functional data from a single run of the flocBLOCKED task (see below) from a single subject (sub-03). The data will be saved in your home directory in a folder named “NI-edu-data”.
import os
data_dir = os.path.join(os.path.expanduser("~"), 'NI-edu-data')
print("Data (+- 110MB) will be saved in %s ..." % data_dir)
print("Go get some coffee. This may take a while.\n")
# Using the CLI interface of the `awscli` Python package, we'll download the data
# Note that we only download the preprocessed data from a single run from a single subject
# This may take about 2-10 minutes or so (depending on your download speed)
# Download anatomical files
!aws s3 sync --no-sign-request s3://openneuro.org/ds003965 {data_dir} --exclude "*" --include "sub-03/ses-1/anat/*"
# And functional files
!aws s3 sync --no-sign-request s3://openneuro.org/ds003965 {data_dir} --exclude "*" --include "sub-03/ses-1/func/*task-flocBLOCKED*"
# And a brain mask (needed for FEAT, later)
!aws s3 sync --no-sign-request s3://openneuro.org/ds003965 {data_dir} --exclude "*" --include "derivatives/fmriprep/sub-03/ses-1/anat/sub-03_ses-1_acq-AxialNsig2Mm1_desc-brain_mask.nii.gz"
print("\nDone!")
Data (+- 110MB) will be saved in /Users/lukas/NI-edu-data ...
Go get some coffee. This may take a while.
Done!
Anyway, let’s check out the structure of the directory with the data. The niedu package has a utility function (in the submodule utils) that returns the path to the root of the dataset:
if os.path.isdir(data_dir):
print("The dataset is located here: %s" % data_dir)
else:
print("WARNING: the dataset could not be found here: %s !" % data_dir)
The dataset is located here: /Users/lukas/NI-edu-data
Contents of the dataset#
Throughout the course, we’ll explain the contents of the dataset and use them for exercises and visualizations. For now, it suffices to know that the full dataset contains data from 13 subjects, each scanned twice (i.e., two sessions acquired at different days). Each subject has it’s own sub-directory, which contains an “anat” directory (with anatomical MRI files) and/or a “func” directory (with functional MRI files).
We can “peek” at the structure/contents of the dataset using the listdir function from the built-in os package (which stands for “Operating System”):
import os
print(os.listdir(data_dir))
['derivatives', 'sub-03']
The sub-03 directory contain the “raw” (i.e., not-preprocessed) data of subject 3. Let’s check out its contents. We’ll do this using a custom function show_tree from the niedu package:
from niedu.global_utils import show_directory_tree
# To create the full path, we use the function `os.path.join` instead of just writing out
# {data_dir}/sub-03. Effectively, `os.path.join` does the same, but takes into account your
# operating system, which is important because Windows paths are separated by backslashes
# (e.g., `sub-03\ses-1\anat`) but Mac and Linux paths are separated by forward slashes
# (e.g., `sub-03/ses-1/anat`).
sub_03_dir = os.path.join(data_dir, 'sub-03')
show_directory_tree(sub_03_dir, ignore='physio')
sub-03/
ses-1/
anat/
sub-03_ses-1_acq-AxialNsig2Mm1_T1w.json
sub-03_ses-1_acq-AxialNsig2Mm1_T1w.nii.gz
sub-03_ses-1_acq-AxialNsig2Mm1_T1w_brain.nii.gz
func/
sub-03_ses-1_task-face_run-1_events.tsv
sub-03_ses-1_task-flocBLOCKED_bold.json
sub-03_ses-1_task-flocBLOCKED_bold.nii.gz
sub-03_ses-1_task-flocBLOCKED_condition-body_events.txt
sub-03_ses-1_task-flocBLOCKED_condition-character_events.txt
sub-03_ses-1_task-flocBLOCKED_condition-face_events.txt
sub-03_ses-1_task-flocBLOCKED_condition-object_events.txt
sub-03_ses-1_task-flocBLOCKED_condition-place_events.txt
sub-03_ses-1_task-flocBLOCKED_condition-response_events.txt
sub-03_ses-1_task-flocBLOCKED_events.tsv
As you can see, the sub-03 directory contains a single session directory, “ses-1”. “Sessions” ususally refer to different days that the subject was scanned. Here, we only downloaded data from a single session (while, in fact, there is data from two sessions) in order to limit the amount of data that needed to be downloaded.
Within the ses-1 directory, there is a func directory (with functional MRI files) and an anat directory (with anatomical MRI files). These directories contain a bunch of different files, with different file types:
(Compressed) nifti files: these are the files with the extension
nii.gz, which contain the actual MRI data. These are discussed extensively in the next section;Json files: these are the files with the extension
json, which are plain-text files containing important metadata about the MRI-files (such as acquisition parameters). We won’t use these much in this course.Event files: these are the files ending in
_events.tsv, which contain information about the onset/duration of events during the fMRI paradigm.
This particular organization of the data is in accordance with the “Brain Imaging Data Structure” (BIDS; read more about BIDS here). If you are going to analyze a new or existing MRI dataset, we strongly suggest to convert it to BIDS first, not only because makes your life a lot easier, but also because there exist many tools geared towards BIDS-formatted datasets (such as fmriprep and mriqc)!
MRI acquisition parameters#
For each participant, we acquired one T1-weighted, high-resolution (\(1 \times 1 \times 1\) mm.) anatomical scan. Our functional scans were acquired using a (T2*-weighted) gradient-echo echo-planar imaging (GE-EPI) scan sequence with a TR of 700 milleseconds and a TE of 30 milliseconds and a spatial resolution of \(2.7 \times 2.7 \times 2.7\) mm. These functional scans used a “multiband acceleration factor” of 4 (which is a “simultaneous multislice” technique to acquire data from multiple slices within a single RF pulse, greatly accellerating acquisition time).
Experimental paradigms#
In this dataset, we used two types of experimental paradigms: a “functional localizer” (which we abbreviate as “floc”) and a simple face perception task (“face”). (We also acquired a resting-state scan, i.e. a scan without an explicit task, but we won’t use that for this course.) The two different floc paradigms (flocBLOCKED and flocER) were each run twice: once in each session. Each session contained 6 face perception runs (so 12 in total).
Functional localizer task#
The functional localizer paradigm we used was developed by the Stanford Vision & Perception Neuroscience Lab. Their stimulus set contains images from different categories:
characters (random strings of letters and numbers);
bodies (full bodies without head and close-ups of limbs);
faces (adults and children);
places (hallways and houses);
objects (musical instruments and cars).
These categories were chosen to target different known, (arguably) category-selective regions in the brain, including the occipital face area (OFA) and fusiform face area (FFA; face-selective), parahippocampal place area (PPA; place/scene selective), visual word form area (VWFA; character-selective) and extrastriate body area (EBA; body-selective).
In the functional localizer task, participants had to do a “one-back” task, in which they had to press a button with their right hand when they saw exactly the same image twice in a row (which was about 10% of the trials).
Importantly, we acquired this localizer in two variants: one with a blocked design (“flocBLOCKED”) and one with an event-related design (“flocER”; with ISIs from an exponential distribution with a 2 second minimum). The blocked version presented stimuli for 400 ms. with a fixed 100 ms. inter-stimulus interval in blocks of 12 stimuli (i.e., 6 seconds). Each run started and ended with a 6 second “null-blocks” (a “block” without stimuli) and additionally contained six randomly inserted null-blocks. The event-related version showed each stimulus for 400 ms. as well, but contained inter-stimulus intervals ranging from 2 to 8 seconds and no null-blocks apart from a null-block at the start and end of each run (see figure below for a visual representation of the paradigms).
These two different versions (blocked vs. event-related) were acquired to demonstrate the difference in detection and estimation efficiency of the two designs. In this course, we only use the blocked version.
Face perception task#
In the face perception task, participants were presented with images of faces (from the publicly available Face Research Lab London Set). In total, frontal face images from 40 different people (“identities”) were used, which were either without expression (“neutral”) or were smiling. Each face image (from in total 80 faces, i.e., 40 identities \(\times\) 2, neutral/smiling) was shown, per participant, 6 times across the 12 runs (3 times per session). We won’t discuss the data from this task, as it won’t be used in this course.
Using FSL FEAT for preprocessing and first-level analysis#
Most univariate fMRI studies contain data from a sample of multiple individuals. Often, the goal of these studies is to (1) estimate some effect (e.g., \(\beta_{face} >0\)) at the group-level, i.e., whether the effect is present within our sample, and subsequently, taking into account the variance (uncertainty) of this group-level effect, (2) to infer whether this effect is likely to be true in the broader population (but beware of the tricky interpretation of null-hypothesis significance tests).
A primer on the “summary statistics” approach#
When you have fMRI data from multiple participants, you are dealing with inherently hierarchical (sometimes called “multilevel”, “mixed”, or simply “nested”) data: you have data from multiple participants, which may be split across different runs or sessions, which each contain multiple observations (i.e., volumes across time).
Often, people use so called “mixed models” (or “multilevel models”) to deal with this kind of hierarchical data. These type of models use the GLM framework for estimating effects, but have a special way of estimating the variance of effects (i.e., \(\mathrm{var}[c\hat{\beta}]\)), which includes variance estimates from all levels within the model (subject-level, run-level, and group-level). In these mixed models, you can simply put all your data (\(y\) and \(\mathbf{X}\)) from different runs, sessions, and subjects, in a single GLM and estimate the effects of interest (and their variance).
Unfortunately, this “traditional” multilevel approach is not feasible for most univariate (whole-brain) analyses, because these models are quite computationally expensive (i.e., they take a long time to run), so doing this for all > 200,000 voxels is the brain is not very practical. Fortunately, the fMRI community came up with a solution: the “summary statistics” approach, which approximates a proper full multilevel model.
This summary statistics approach entails splitting up the analysis of (hierarchical) fMRI data into different steps:
Analysis at the single-subject (time series) level
Analysis at the run level (optional; only if you have >1 run!)
Analysis at the group level
The GLM analyses you did in week 2-4 were examples of “first-level” analyses, because they analyzed the time series of voxels from a single subject. This is the first step in the summary statistics approach. In this section, we’re going to demonstrate how to implement first-level analyses (and preprocessing) using FSL FEAT.
Later in this lab, we’re going to demonstrate how to implement analyses at the run-level. Next week, we’ll discuss group-level models. For now, let’s focus on how to run first-level analyses in FSL using FEAT.
First-level analyses (and preprocessing) in FEAT#
Alright, let’s get started with FEAT. You can start the FEAT graphical interface by typing Feat & in the command line (the & will open FEAT in the background, so you can keep using your command line).*
* Note that the upcoming ToDos do not have test-cells; at the very end of this section, you’ll save your FEAT setup, which we’ll test (with only hidden tests) afterwards.
The “Data” tab#
After calling Feat & in the command line, the FEAT GUI should pop up. This GUI has several “tabs”, of which the “Data” tab is displayed when opening FEAT. Below, we summarized the most important features of the “Data” tab.
print(sub_03_dir)
/Users/lukas/NI-edu-data/sub-03
We are going to analyze the flocBLOCKED data from the first session (ses-1): sub-03/ses-1/sub-03_ses-1_task-flocBLOCKED_bold.nii.gz.
sub_03_file = os.path.join(sub_03_dir, 'ses-1', 'func', 'sub-03_ses-1_task-flocBLOCKED_bold.nii.gz')
print("We are going to analyse: %s" % sub_03_file)
We are going to analyse: /Users/lukas/NI-edu-data/sub-03/ses-1/func/sub-03_ses-1_task-flocBLOCKED_bold.nii.gz
Let’s fill in the Data tab in FEAT. Click on “Select 4D data”, click on the “folder” icon, and navigate to the file (sub-03_ses-1_task-flocBLOCKED_bold.nii.gz). Once you’ve selected the file, click on “OK”. Note: within FEAT, you can click on .. to go “up” one folder. This way, you can “navigate” to the location of the file that you want to add.
Then, set the output directory to your home folder + sub-03_ses-1_task-flockBLOCKED. Note: this directory doesn’t exist yet! By setting the output directory to this directory, you’re telling FSL to create a directory here in which the results will be stored.
After selecting the 4D data, the “Total volumes” and “TR” fields should update to 325 and 0.700000, respectively. You can leave the “Delete volumes” field (0) and high-pass filter (100) as it is (no need to remove volumes and 100 second high-pass is a reasonable default).
The “Pre-stats” tab#
Now, let’s check out the “pre-stats” tab, which contains most of the preprocessing settings. Most settings are self-explanatory, but if you want to know more about the different options, hover you cursor above the option (the text, not the yellow box) for a couple of seconds and a green window will pop up with more information about that setting!
Note that the checkbox being yellow means that the option has been selected (and a gray box means it is not selected).
Ignore the “Alternative reference image”, “B0 unwarping”, “Intensity normalization”, “Perfusion subtraction”, and “ICA exploration” settings; these are beyond the scope of this course!
Let’s fill in the Pre-stats tab! Select the following:
Select “MCFLIRT” for motion-correction (this is the name of FSL’s motion-correction function);
Set “Slice timing correction” to “None” (we believe that slice-time correction does more harm than good for fMRI data with short TRs, like ours, but “your mileage may vary”);
Turn on “BET brain extraction” (this refers to skullstripping of the functional file! This can be done automatically, because it’s less error-prone than skullstripping of anatomical images);
Set the FWHM of the spatial smoothing to 4 mm;
Turn off “intensity normalization” (you should never use this for task-based fMRI analysis);
Turn on “highpass filtering”; leave “perfusion subtraction” unchecked;
Leave “ICA exploration” unchecked
The “Registration” tab#
Alright, let’s check out the “Registration” tab, in which you can setup the co-registration plan (from the functional data to the T1-weighted anatomical scan, and from the T1-weighted anatomical scan to standard space; we usually don’t have an “expanded functional image”, so you can ignore that).
Now, for the “Main structural image”, you actually need a high-resolution T1-weighted anatomical file that has already been skullstripped. You could do that with, for example, bet for this particular partipant. Below, we do this slightly differently, by applying an existing mask, computed by Fmriprep, to our structural image. We use the Python package Nilearn for this (you will learn more about this package in week 7)! No need to fully understand the code below.
from nilearn import masking, image
# Define anatomical file ...
anat = os.path.join(sub_03_dir, 'ses-1', 'anat', 'sub-03_ses-1_acq-AxialNsig2Mm1_T1w.nii.gz')
# ... and mask
mask = os.path.join(data_dir, 'derivatives', 'fmriprep', 'sub-03', 'ses-1', 'anat', 'sub-03_ses-1_acq-AxialNsig2Mm1_desc-brain_mask.nii.gz')
# Affine of mask is different than original anatomical file (?),
# so let's resample the mask to the image
mask = image.resample_to_img(mask, anat, interpolation='nearest')
# Set everything outside the mask to 0 using the apply_mask/unmask trick
anat_sks = masking.unmask(masking.apply_mask(anat, mask), mask)
# And finally save the skullstripped image in the same directory as the
# original anatomical file
out = anat.replace('.nii.gz', '_brain.nii.gz')
if not os.path.isfile(out):
anat_sks.to_filename(out)
You can leave the file under “Standard space” as it is (this is the MNI152 2mm template, one of the most used templates in the fMRI research community). Turn on the “Nonlinear” option in the “Standard space” box and do not change the default options of non-linear registration.
The “Stats” tab#
Almost there! Up next is the “Stats” tab, in which you specify the design (\(\mathbf{X}\)) you want to use to model the data.
Then, to specify our predictors-of-interest, we need to click the “Full model setup” button. This creates a new window in which we can create our predictors; see the image below for an explanation of the options.
The most often-used option to define the event onsets and durations is the “Custom (3 column format)” option. This allows you to select a standard text-file with three columns (separated by tabs), in which the first column indicates the event onsets (in seconds), the second column indicates the event duration (in seconds), and the third column indicates the event “weight” (which may vary in case of parametric designs; for non-parametric designs, you should set the weight-value to 1). Different rows define different events of a particular condition.
Our data only has a BIDS-compatible “events” file, which contain the onset, duration, and trial type of trials from all conditions, while FSL requires a separate events file per condition (and without a header). Below, we run a function that creates these separate condition files for us:
from niedu.utils.nii import create_fsl_onset_files
create_fsl_onset_files(data_dir)
Let’s set up the design! For our flocBLOCKED run, we should model each of the image conditions (“face”, “body”, “place”, “character”, and “object”) as separate predictors, assuming that they have a different effect on voxels (in different brain regions). Let’s model the responses (button presses) as well.
So, create six predictors with the names “face”, “body”, “place”, “character”, “object”, and “response” (in this exact order). To save you some work, we have already created the onset-files associated with these six different conditions (which end in: *condition-{name of condition}_events.txt).
For each predictor:
set the convolution to “Double-Gamma HRF”
turn off the “Add temporal derivative” and “Orthogonalise” options
leave the “Apply temporal filtering” option checked (this applies the high-pass filter to the design, \(X\), as well).
YOUR ANSWER HERE
If you did it correctly, you can press the “View design” button and a new window with the design should pop up (with the different predictors as squiggly lines from top top bottom). The red line of the left side of the design matrix represents the (length of the) high-pass filter.
YOUR ANSWER HERE
Now let’s create some contrasts! In the “Contrasts and F-tests” tab, you can specify the contrasts that you want to test. Importantly, FSL makes a difference between “Original EVs” and “Real EVs”. Original EVs are the predictors-of-interest that you defined in the “EVs” tab. Real EVs are the original EVs + the extra predictors that are generated due to the specific HRF-model. For example, if you’d choose a gamma basis set instead of a double-gamma HRF, the original canonical HRF of each predictor would be “Original EVs” while all predictors making up the HRF-components (canonical and temporal and dispersion derivatives) would be the “Real EVs”. In our case, the original and real EVs are the same.
To create contrasts (and F-tests, i.e., collections of contrasts), click on the “Contrasts & F-tests” tab. Here, you define the actual contrast vectors. For each contrast, you can give it a name (under “Title”) and set the contrast weight for each predictor (or “EV”). To create F-tests, set the number of F-tests that you want to do (i.e., increase the number next to the “F-tests” field), and check, per F-test, which contrasts to want to include in them. For example, if you want to include contrasts 1 and 2 in F-test, make sure that these checkboxes are selected under “F1” (or “F2”, “F3”, etc.).
Given the six original EVs (“face”, “body”, “place”, “character”, “object”, and “response”), define the contrasts corresponding to the following null (\(H_{0}\)) and alternative (\(H_{a}\)) hypotheses (give them whatever name you want, but adhere to this order):
\(H_{0}: \beta_{face} = 0\)
\(H_{a}: \beta_{face} > 0\)\(H_{0}: \beta_{place} = 0\)
\(H_{a}: \beta_{place} > 0\)\(H_{0}: \beta_{response} = 0\)
\(H_{a}: \beta_{response} > 0\)\(H_{0}: 4\times\beta_{face} - \beta_{body} - \beta_{place} - \beta_{character} - \beta_{object} = 0\)
\(H_{a}: 4\times\beta_{face} - \beta_{body} - \beta_{place} - \beta_{character} - \beta_{object} > 0\)\(H_{0}: 4\times\beta_{place} - \beta_{body} - \beta_{face} - \beta_{character} - \beta_{object} = 0\)
\(H_{a}: 4\times\beta_{place} - \beta_{body} - \beta_{face} - \beta_{character} - \beta_{object} > 0\)\(H_{0}: \beta_{face} = \beta_{body} = \beta_{place} = \beta_{character} = \beta_{object} = 0\)
\(H_{a}: (\beta_{face} > 0)\ \&\ (\beta_{body} > 0)\ \&\ (\beta_{place} > 0)\ \&\ (\beta_{character} > 0)\ \&\ (\beta_{object} > 0)\)\(H_{0}: \beta_{face} = \beta_{body} = \beta_{place} = \beta_{character} = \beta_{object} = 0\)
\(H_{a}: (\beta_{face} < 0)\ \&\ (\beta_{body} < 0)\ \&\ (\beta_{place} < 0)\ \&\ (\beta_{character} < 0)\ \&\ (\beta_{object} < 0)\)
Also, define the following F-test (where the vertical line stands for “or”):
\(H_{0}: \beta_{face} = \beta_{place} = 0\)
\(H_{a}: (\beta_{face} \neq 0)\ |\ (\beta_{place} \neq 0)\)
After you’re done defining the contrasts and F-test, you can press “View design” again, which will now also include the contrasts/F-tests that you defined. Moreover, if you click on “Efficiency”, a new window opens with information about the effiency of the design. On the left, a correlation matrix of the predictors is shown, where darker off-diagonal cells means lower correlation between predictors (which is a good thing!). On the right, the diagonal of the matrix shows you the “eigenvalues” of the singular value decomposition of the design matrix, which should decrease from top-left to bottom-right. In general, the brighter the values, the more efficient the design matrix. In the lower-right corner, it lists the estimated effect size (as % signal change) required for a significant effect. These estimates depend on your design efficiency (\(\mathbf{c}(\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{c}^{T}\)) and the assumed level of noise (\(\hat{\sigma}^2\)): the lower these values, the better!
The “Post-stats” tab#
Alright, the last tab! The “Post-stats” tab contains options to threshold the resulting statistics-maps. The only interesting options for our purposes are those in the “Thresholding” box. Technically, we don’t need to threshold any of our statistics, because we are only interested in inference at the group-level (not the first-level analysis; we’re not aiming to infer anything from the results of a single-subject!). But using a p-value threshold of 0.05 “uncorrected” (which means: not corrected for multiple comparisons — a topic we’ll discuss next week!) is usually used to detect any potential problems during the experimental, preprocessing, and analysis process. For example, if we don’t see significant visual cortex activation in response to a task with visual stimuli at \(p < 0.05\) (uncorrected), then it is very likely that something went wrong somewhere in the preprocessing/analysis pipeline (e.g., stimulus onsets are incorrect).
Alright, we’re done setting up FEAT for our first-level analysis! One last thing: let’s save the setup (which we’ll check and grade).
We actually ran the first-level analysis for this participant already, so you don’t have to do this yourself! We’ll check it out in the next section. (The following cells all test your setup_feat.fsf file; as they all have hidden test, you may ignore them.)
''' Tests your setup_feat.fsf file: initial check + Data tab'''
import os
import os.path as op
par_dir = op.basename(op.dirname(op.abspath('.')))
if par_dir != 'solutions': # must be a student-account
fsf = 'setup_feat.fsf'
if not op.isfile(fsf):
raise ValueError("Couldn't find a 'setup_feat.fsf' file in your week_5 directory!")
print("Only hidden tests!")
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
/var/folders/0c/2qv1yrms7f56lxzzw28nt_440000gn/T/ipykernel_15549/4285188295.py in <cell line: 6>()
7 fsf = 'setup_feat.fsf'
8 if not op.isfile(fsf):
----> 9 raise ValueError("Couldn't find a 'setup_feat.fsf' file in your week_5 directory!")
10
11 print("Only hidden tests!")
ValueError: Couldn't find a 'setup_feat.fsf' file in your week_5 directory!
''' Tests your setup_feat.fsf file: Pre-stats tab'''
print("Only hidden tests!")
Only hidden tests!
''' Tests your setup_feat.fsf file: Registration tab'''
print("Only hidden tests!")
Only hidden tests!
''' Tests your setup_feat.fsf file: Stats tab (filt) '''
print("Only hidden tests!")
Only hidden tests!
''' Tests your setup_feat.fsf file: Stats tab (EVs) '''
print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Stats tab, (con 1)'''
print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Stats tab, (con 2)'''
print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Stats tab, (con 3)'''
print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Stats tab, (con 4)'''
print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Stats tab, (con 5)'''
print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Stats tab, (con 6)'''
print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Stats tab, (con 7)'''
print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Stats tab (F-test)'''
print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Post-stats tab'''
print("Only hidden tests!")
Only hidden tests!
Making sense of FEAT output#
When FEAT is running the analyis, it will create a results-directory (at the location you specified) ending in .feat, which contains all the information and files necessary for the analysis as well as the results. Fortunately, it also generates an informative report.html file with a nice summary of the preprocessing and analysis results. We already ran a first-level analysis like with the setup that you specified in the previous section, which we download below:
print("Starting download of FEAT directory (+- 285MB) ...\n")
!aws s3 sync --no-sign-request s3://openneuro.org/ds003965 {data_dir} --exclude "*" --include "derivatives/fsl/sub-03/flocBLOCKED/ses-1.feat/*"
print("\nDone!")
Starting download of FEAT directory (+- 285MB) ...
Done!
The downloaded data is located here:
feat_dir = os.path.join(data_dir, 'derivatives', 'fsl', 'sub-03', 'flocBLOCKED', 'ses-1.feat')
print(feat_dir)
/Users/lukas/NI-edu-data/derivatives/fsl/sub-03/flocBLOCKED/ses-1.feat
firefox report.html
Checking the report.html file#
The “Registration” tab shows you several figures related to the (two-step) co-registration process.
In each figure, the red outline shows the target (i.e., where the data is registered to). Each figure contains two subplots, which show the co-registration in different “directions”. The reason they plot both “directions” is that a single “direction” may give you a little information about whether the registration completed without problems.
The first image shows you the result from the complete registration process: from our functional file (“FMRI”) to the MNI152 template (“standard space”). Here, the target (in red) is the outline of the MNI152 template. The other figures show the intermediate results in the two-step co-registration process. The second figure (“Registration of example_func to highres”) shows you first registration step: the BBR-based registration of (a reference volume of) our functional data to the high-resolution T1-weighted scan. Note that, in the upper row, the red outline represents the white matter boundaries of our T1-weighted scan, because this is wat BBR uses as the registration target (rather than the entire volume).
The third image show you the second registration step: the affine + non-linear registration of our high-resolution T1-weighted scan to the MNI152 template. The fourth image shows you the result of the entire two-step registration procedure (functional → template; essentially the same as the first image).
This overview is an excellent way to check whether the results are roughly what you’d expect. For example, you’d expect to see clear activation in visual cortex for our the contrasts against baseline (zstat1: face > 0 and zstat2: place > 0 and zstat6: stim).
Apart from the z-statistic maps, FSL shows you for each contrast the time series of the voxel that has the highest z-value. The plots show the data (i.e., \(y\); in red), the model-fit (i.e., \(\hat{y}\); in blue/purple), and the “partial model fit” from that contrast (i.e., \(X_{c}\hat{\beta}_{c}\) for that particular contrast \(c\); in green).
Checking out the results from FEAT (in FSLeyes)#
While the report.html file provides a nice summary of the results, it’s also good to visualize and look at output of FEAT is more detail. First, we’ll walk you through the contents of the FEAT directory and then we’ll take a look at the results in FSLeyes.
You should see a large set of different files: text-files, nifti-files, and subdirectories. Many of these files represent intermediary output from FEAT or data/images necessary for the report*.html files (which we inspected earlier). The most important subdirectories are the reg and stats directories. These contain the registration parameters and GLM results, respectively, that we later need for our fixed-effects and group-level analyses.
The .mat files contain affine matrices for the different registration steps: from functional to the high-resolution anatomical scan (example_func2highres.mat) and from the high-resolution anatomical scan to the standard template (highres2standard.mat) (and the inverse matrices: highres2example_func.mat and standard2example_func.mat). In addition, it contains the “warpfield” (files ending in _warp.nii.gz) computed by the non-linear registration step, i.e., how much each voxel needs to be moved in each direction (\(X, Y, Z\)).
Basically, all the nifti files you see here are 3D files (“brain maps”), representing different parts of the t-value formula, as is visualized in the figure below:
In fact, to show you that the these files are simply parts of the \(t\)-value formula, let’s calculate a t-stat map ourselves using the cope and varcope files.
Below, we load the cope1.nii.gz (containing the \(\mathbf{c}\hat{\beta}\) values per voxel for contrast 1) and varcope1.nii.gz (containing the \(\hat{\sigma}^{2}\mathbf{c}(\mathbf{X}^{T}\mathbf{V}^{-1}\mathbf{X})^{-1}\mathbf{c}^{T}\) values per voxel for contrast 1). Calculate the corresponding 3D tstat map (containing all t-values per voxel for contrast 1) using numpy array operations, so no for loop! Store this in a variable named tstat1.
Hint: all voxels outside the brain (in both the cope1.nii.gz and the varcope.nii.gz files) will have the value 0. When calculating the t-values, this will lead to nan values (“not a number) because you’re dividing 0 by 0 (which will give a RuntimeWarning). Make sure set these nan values back to 0 (the function np.isnan might be useful ….).
# Implement your ToDo here!
import os
import nibabel as nib
import numpy as np
cope1_path = os.path.join(feat_dir, 'stats', 'cope1.nii.gz')
cope1 = nib.load(cope1_path).get_fdata()
varcope1_path = os.path.join(feat_dir, 'stats', 'varcope1.nii.gz')
varcope1 = nib.load(varcope1_path).get_fdata()
# YOUR CODE HERE
raise NotImplementedError()
''' Tests the above ToDo. '''
# Some (visible) tests to check minor mistakes
if np.nan in tstat1:
raise ValueError('There are still nans in your data!')
from niedu.tests.nii.week_5 import test_tvalue_from_fsl
test_tvalue_from_fsl(tstat1, feat_dir)
Alright, now let’s continue with the more interesting part: visualizing the results.
FSLeyes will automatically open the filtered_func_data.nii.gz (the preprocessed functional data) file in the Ortho view.
The time series viewer now shows you both the data (\(y\)) and the “full model fit” (\(\hat{y}\),) for the voxel highlighted by the crosshairs.
Let’s load some statistic files to see the (whole-brain) effects of a particular contrast.
In the time series viewer, we can also add additional time series, such as the fitted “COPE” corresponding to our first contrast. The fitted “COPE” basically represents the predicted time series for a particular contrast, i.e., \(\mathbf{X}\mathbf{c}^{T}\mathbf{c}\hat{\beta}\). (Usually, however, the intercept added to the predicted time series regardless of whether it’s included in the contrast.)
For example, suppose that we want to compute the fitted COPE corresponding to the \(\beta_{face} > 0\) contrast, we would do the following:
# Let's load in the actual design matrix from the FEAT dir
dm_path = os.path.join(feat_dir, 'design.mat')
X = np.loadtxt(dm_path, skiprows=5)
X = X[:, :6] # remove motion preds for clarity of example
# Define face > 0 contrast vector
# Add new axis to make sure matrix product works
# Order: face, body, place, character, object, response
c = np.array([1, 0, 0, 0, 0, 0])[np.newaxis, :]
# Create some hypothetical beta vector
est_betas = np.array([3.3, 2.5, 1.4, -2.1, 0.8, 0.1])
fitted_cope_face = X @ c.T @ c @ est_betas
plt.figure(figsize=(15, 3.5))
plt.plot(fitted_cope_face)
plt.xlim(0, X.shape[0])
plt.ylabel("Activity (A.U.)", fontsize=20)
plt.xlabel("Time (volumes)", fontsize=20)
plt.title(r"Fitted COPE for $\beta_{face} > 0$", fontsize=25)
plt.ylim(-2.5, 3.5)
plt.grid()
In the time series viewer, you should see a new line corresponding to “COPE1 fit: face”.
This is exactly the reason why we also defined the \(4\times\beta_{face} - \beta_{body} - \beta_{place} - \beta_{character} - \beta_{object} > 0\) contrast (contrast nr. 4 in our analysis)! This contrast allows us to identify regions in the brain that are (relatively more) specific for faces than for other types of images. In fact, in the context of our functional localizer task, you could argue that the \(\beta_{face} > 0\) contrast is largely useless.
Many voxels throughout the brain seem to specifically activate in response to images of faces (but beware: these brain images are not yet “properly” thresholded, so many of these voxels might actually be false positives, but that’s the topic of next week!). Two cluster of voxels, however, seem especially strongly activated: one centered at \([24, 19, 21]\) and another one centered at \([23, 30, 13]\). Note that these clusters are present in both hemispheres, but the cluster in the right hemisphere is clearly more strongly activated, which is consistent with the apparent right-hemisphere bias in face processing (see e.g., Kanwisher et al., 1997).
YOUR ANSWER HERE
This contrast shows you which voxels are more active during periods in which stimuli were absent than when stimuli were presented (regardless of condition). The pattern of voxels that are often found to be more active during “rest” than during (external) “stimulation” has been named the default mode network (DMN), which may reflect daydreaming, episodic or autobiographical information processing, or self reflection (amongst others; the exact “function” of the DMN is still unclear).
YOUR ANSWER HERE
This sigmasquareds map (and the underlying “raw” residual time series) can given you a good overview of the noise in your data and potential violations of GLM (e.g., unequal or correlated error variance). This is a great way to check for unexpected issues (e.g., artifacts) with your data, that will, if unmodelled, end up in your noise term.
This particular voxel lies most likely within the Circle of Willis: a “hub” of multiple arteries supplying blood to the brain. As BOLD MRI measures the (relative) amount of oxygenated blood in the brain, it is no surprise that arteries (and veins!) produce a strong heartbeat-related signal, which will, if left unmodelled, make up much of the unexplained variance. Indeed, in this sigmasquareds map, most of voxels with high values (i.e., those that are relatively yellow) are located within or near arteries or veins. For example, you can clearly “see” the sagittal sinus (see image below).
From Wikimedia commons (link)
Apart from the veins/arteries, the cerebrospinal fluid (CSF) also seem to contain a lot of unmodelled noise. In fact, most of this unmodelled variance also likely represents the pulsative cerebral bloodflow, as shown (exaggerated) in the image below:
Now, move your crosshair across different voxels with high values (i.e., relatively yellow ones), for example voxel \([37, 48, 17]\) (located in the CSF). You can see that this voxel’s residual time series is highly (but not perfectly) correlated to the original voxel in the Circle of Willis, the common cause being the pulsative cerebral blood flow. (Note, though, that different voxels may show time series that look slightly shifted in time relative to the original Circle of Willis voxel, especially in more superior regions of the brain, which is caused by different slice onsets!)
Also note that the absolute signal amplitude of voxels within veins/arteries is much higher than voxels in gray matter!
Alright, time for something else: let’s delve into how you can “combine” run or session data in the next notebook (run_level_analyses.ipynb).