This tutorial steps you through the analysis of an fMRI data set with the FreeSurfer Functional Analysis Stream (FSFAST), from organizing the data to group analysis.
- FreeSurfer Slides
- Tutorial Data Description
- Getting and Organizing the Tutorial
- Quick Visualization Tutorial (tkmedit/tksurfer)
- Assembling the Data into the FSFAST Hierarchy
- Preprocessing of fMRI Data (preproc-sess)
- Function-Structure Registration
- First-Level Analysis
- Higher-Level (Group) Analysis
- Play with Your Data
Tutorial Data Description
The data being analyzed are part of the fBIRN Phase I data set (www.nbirn.net).
- 5 subjects
- Each scanned twice
- Sensory-Motor Block Task; 15 sec OFF, 15 sec ON; 8 cycles
- 4 runs/visit
- TR = 3 sec, 85 time points/run (run length = 255 sec)
- 10 sites in full data set (tutorial data only from 1st visit at MGH site)
Getting and Organizing the Tutorial
If you do not have the tutorial data set up, then consult the FsFastTutorialData page. You will need to set the FSFTUTDIR environment variable. NOTE: if you are taking the Multi-Modal Imaging Class at MGH, the data have already been set up.
Cd into the tutorial data directory and set SUBJECTS_DIR to use this data:
cd $FSFTUTDIR setenv SUBJECTS_DIR $PWD
It will be assumed that you are in this directory (or subdirectory thereof) throughout the tutorial. This directory will have five folders:
- fb1-raw - raw fMRI data
- fb1-raw-study - raw data organized as an FSFAST study but unanalyzed
- fb1-preproc-study - The same data preprocessed.
- fb1-analysis-study - The same data analyzed
subjects - FreeSurfer reconstruction of anatomical data (plus the fsaverage subject).
If you look in fb1-raw you will see 20 NIFTI data sets:
ls fb1-raw f.mgh-101.1-r003.nii f.mgh-103.1-r003.nii f.mgh-104.1-r003.nii f.mgh-105.1-r003.nii f.mgh-106.1-r003.nii f.mgh-101.1-r005.nii f.mgh-103.1-r005.nii f.mgh-104.1-r005.nii f.mgh-105.1-r005.nii f.mgh-106.1-r005.nii f.mgh-101.1-r007.nii f.mgh-103.1-r007.nii f.mgh-104.1-r007.nii f.mgh-105.1-r007.nii f.mgh-106.1-r007.nii f.mgh-101.1-r010.nii f.mgh-103.1-r010.nii f.mgh-104.1-r010.nii f.mgh-105.1-r010.nii f.mgh-106.1-r010.nii
These are the 20 fMRI runs mentioned above. In the name, f.mgh-10X.1-rYYY.nii, X indicates the subject number (1, 3, 4, 5, 6), and YYY indicates the run number. The sensory-motor task happened to be runs 3, 5, 7, and 10.
Quick Visualization Tutorial (tkmedit/tksurfer)
The purpose of this tutorial is to familiarize you with how to use FreeSurfer volume viewer (tkmedit) and surface viewer (tksurfer) in the context of viewing functional data. You should otherwise already know how to use tkmedit and tksurfer. See the pages below for a more detailed handling of tkmedit and tksurfer.
More info on tkmedit
More info on tksurfer
Viewing a single functional overlay in the volume
cd into the study with already analyzed data:
Run tkmedit-sess (this is an FSFAST wrapper for tkmedit):
tkmedit-sess -s mgh-101.1 -a sm-gamma-fwhm5 -c odd-v-0 -aparc+aseg
Don't worry about what all the arguments mean, as this part is only about visualization. This command will bring up two windows: one with a brain image, the other a control panel.
The brain image is the FreeSurfer anatomical for this subject.
The slightly pale colors on the anatomical indicate the FreeSurfer automatic segmentation.
- The bright red/yellow/blue are super-threshold voxels in the functional overlay.
- As you click on different points, you will see the "Functional value" field in the Tools window change as well as the "Sgmtn label".
- The interpretation of the value depends on what is being viewed. This is a significance map, so the value is -log10(pvalue)*sign (i.e., for a pvalue = .01, the functional value will be +2).
- The sign is a functional direction. The red/yellow are positive, and the blue are negative. As a functional value gets more positive, its color will change from red to yellow. As it gets more negative, it will change from blue to cyan.
- Note that areas that are not above threshold will still have non-zero functional values.
Toggle the functional overlay on and off by hitting the button with the red/yellow blob in the Tools window (it's in the top row - if you mouse over it, you'll see a "Show Functional Overlay" tooltip).
Configure the functional overlay by clicking on "View..." in the Tools window, then "Configure->Functional Overlay...". You should see the following interface:
- The thresholds are currently set at 2 (Min) and 4 (max).
- The Min threshold is the minimum absolute value needed for a voxel to show an overlay color.
- The maximum is the value beyond which the voxel will stop changing color.
Try changing these values, then hit the Apply button to see their effect.
So, why does it look like there's activation OUTSIDE of the brain? Try hitting View->Mask Functional Overlay to Aux button to see a marked improvement. And why is it so blocky? There are big voxels in functional data, but if you'd like to pretend otherwise, try hitting the "Trilinear" button on the Configure Functional Overlay window.
Viewing multiple functional overlays and time courses in the volume
Exit from your current tkmedit window. Run tkmedit-sess (this is an FSFAST wrapper for tkmedit):
tkmedit-sess -s mgh-101.1 -a sm-fir-fwhm5 -c odd-v-0 -aparc+aseg
This will bring up the two windows as you saw before (the brain image window will have a different overlay). It will also bring up a window called "Time Course". Click anywhere in the volume, and you will see a time course associated. If you happened to have clicked in some activation, then the time course will look similar to the one below.
The meaning of the time course depends upon the nature of the time course loaded. In this case it is the hemodynamic response to the task block averaged over all blocks over all runs (this is an "FIR" model). The horizontal axis is time (0 means the onset of the block). You can visualize the values of any multi-frame data set in this way (it does not have to be a "time" course).
Notice that there is not much activation in the functional overlay. This is good! You are looking at an overlay that corresponds 3 seconds prior to stimulus onset, so there should be no activation.
To see the other time points, bring up the functional overlay configuration window (View->Configure->Functional Overlay). Notice that the "Time Point" field now has a range of 0-8 indicating the 9 time points in the Time Course window. If you hit the + button next to "Time Point" then hit Apply, you will see the overlay change.
You will also see a vertical dashed line in the Time Course window. You are now looking at the map associated with the time between 0 and 3 sec after stimulus onset.
Keep hitting the +/Apply buttons to see different time points. You can hit the "|>" button to start a movie of the activation (then "" to stop it).
Viewing a single functional overlay on the surface
To view the same data on the surface, start a new sheel and run:
cd $FSFTUTDIR/fb1-analysis-study tksurfer-sess -s mgh-101.1 -a sm-gamma-fwhm5 -c odd-v-0 -hemi lh -aparc
Again you will see two windows, "Tksurfer Tools" and a surface image window like the one below. As with tkmedit, you can configure the overlay with View->Configure...->Overlay. The second command will load the surface parcellation which will hide the overlay. Change it to outline mode with: View->Label Style->Outline.
You can view the time course with:
tksurfer-sess -s mgh-101.1 -a sm-fir-fwhm5 -c odd-v-0 -hemi lh -aparc
Assembling the Data into the FSFAST Hierarchy
Before analyzing data in FSFAST, it needs to be in a "Hierarchy", i.e., a certain directory structure and naming convention. By adhering to a hierarchy, you allow FSFAST to analyze large amounts of data automatically because it knows where to find the input and where to put the output.
- Each step in the analysis modifies the hierarchy in some way.
- All the data collected in one functional session (or visit) are placed under the folder (the "session directory").
- All the sessions are arranged under one folder (called the "Study Directory").
The full hierarchy is shown in the image below:
The red boxes indicate a directory with raw or preprocssed functional data. E.g., f.nii is a raw fMRI 4D file, and fmc.nii is motion corrected.
Creating a Hierarchy from the Tutorial Data
There are five subjects of data, each with four runs. We will create a an FSFAST session from one of the subjects. A properly created hierarchy can be found in $FSFTUTDIR/fb1-raw-study.
Create a Study Directory and cd into it:
cd $FSFTUTDIR mkdir -p my-study cd my-study
Create a session directory for the first subject:
mkdir -p mgh-101.1
Create a bold and four run directories:
mkdir -p mgh-101.1/bold mkdir -p mgh-101.1/bold/003 mkdir -p mgh-101.1/bold/005 mkdir -p mgh-101.1/bold/007 mkdir -p mgh-101.1/bold/010
Copy the raw data:
cp ../fb1-raw/f.mgh-101.1-r003.nii mgh-101.1/bold/003/f.nii cp ../fb1-raw/f.mgh-101.1-r005.nii mgh-101.1/bold/005/f.nii cp ../fb1-raw/f.mgh-101.1-r007.nii mgh-101.1/bold/007/f.nii cp ../fb1-raw/f.mgh-101.1-r010.nii mgh-101.1/bold/010/f.nii
Note that they are all named "f.nii" but are distinguished by the directory they are in.
Link to FreeSurfer Anatomical Analysis
FSFAST makes extensive use of the FreeSurfer anatomical analysis. To link a subject's functional and structural data, create a "subjectname" file in the session directory. This file is literally called "subjectname" and its contents have the FreeSurfer subject name assigned during reconstruction and found in $SUBJECTS_DIR. These subjects have already been reconstructed; this subject's name is "fbph1-101". You can create this file with any text editor (emacs, vi, pico), or just:
cd $FSFTUTDIR/my-study echo fbph1-101 > mgh-101.1/subjectname
Create a Stimulus Schedule (Paradigm File)
A "paradigm" file is a record of which stimulus was presented when and for how long. Each paradigm file has four columns (extras are ignored). The columns are:
- 1. Stimulus onset time (sec)
- 2. Condition ID code (0, 1, 2, ...)
- 3. Stimulus Duration (sec)
- 4. Stimulus Weight (usually 1)
The Stimulus onset is with respect to the time of the first image. The Condition is a unique number that identifies the event or block type that the stimulus belongs to. They must be contiguous, and Condition 0 is reserved for the Null Event (usually Fixation). The weight can be used to model various effects of learning or adaptation (just set them to 1 for now).
For this exercise, you will create two paradigm files based on the fBIRN sensory-motor experiment. The timing for the experiment was described above. The task timing is shown graphically below.
Using a text editor (emacs, vi, pico, etc), create a paradigm file called "sensory-motor0.par". Eg: cd $FSFTUTDIR/fb1-raw-study emacs sensory-motor0.par&
The correct paradigm file is here SensMotPar (This example file has a 5th column with the condition name. This is ignored in the analysis).
To make things a little more interesting, we are going to create another paradigm file pretending that the odd and even blocks are different. Assign Odd blocks Condition 1 and Even blocks Condition 2. cd $FSFTUTDIR/fb1-raw-study emacs sensory-motor.par&
Call this paradigm file "sensory-motor.par". The correct paradigm file here is SensMotOddEvenPar.
Now copy this file into each run directory:
cd $FSFTUTDIR/my-study cp ../fb1-raw-study/sensory-motor.par mgh-101.1/bold/003/sensory-motor.par cp ../fb1-raw-study/sensory-motor.par mgh-101.1/bold/005/sensory-motor.par cp ../fb1-raw-study/sensory-motor.par mgh-101.1/bold/007/sensory-motor.par cp ../fb1-raw-study/sensory-motor.par mgh-101.1/bold/010/sensory-motor.par
While each of these have the same contents, they can just as easily be different.
Create a Session ID File
You would normally do the above steps for each of the subjects, but you've probably got the point by now. If you were to create a session for each, you would then have five sessions: mgh-101.1, mgh-103.1, mgh-104.1, mgh-105.1, mgh-106.1. To help keep track of your sessions, you can create a session ID file with your favorite text editor. This is just a simple file with a list of your sessions. The contents should look like (the order is unimportant):
mgh-101.1 mgh-103.1 mgh-104.1 mgh-105.1 mgh-106.1
This file can be passed to many of the FSFAST commands with the -sf argument. Individual sessions can be passed with -s, as you will see below.
Preprocessing of fMRI Data (preproc-sess)
What preprocessing stages do you want to run?
There are at least seven types of preprocessing that are run on fMRI data:
- 1. motion correction (MC)
- 2. spatial smoothing
- 3. brain masking
- 4. slice-timing correction (STC)
- 5. intensity normalization
- 6. resampling to common space
- 7. B0 distortion correction.
Not all packages run all of these seven, and they are not always run in the same order, and some stages are sometimes run as post-processing. In FSFAST, we can run 1. MC, 2. smoothing 3. STC (optional), and 4. intensity normalization (optional). We create a brain mask, but we do not mask the functional data.
For these exercises, cd into the fb1-raw-study directory:
Note that the fully preprocessed data can be found in fb1-preproc-raw.
To run MC and spatial smoothing by 5 mm FWHM along with brain mask creation on one session run:
cd $FSFTUTDIR/my-study preproc-sess -s mgh-101.1 -fwhm 5
Note that this has already been done with all subjects in the fb1-preproc-study directory. To preprocess all of the session, you would run "preproc-sess -sf sessid -fwhm 5". Also note that if preprocessing has already been performed on a session, it will automatically skip it and move on to the next (unless you add -force to the command-line).
Examine additions to FSFAST hierarchy
As mentioned above, each processing step adds something to the FSFAST hierarchy. To see what has been created, run:
cd $FSFTUTDIR/my-study ls mgh-101.1/bold/003
You will see:
- f.nii (the raw data)
- fmc.nii (motion corrected)
- fmcsm5.nii (motion corrected and smoothed).
- fmc.mcdat - text file with the motion correction parameters (AFNI)
You will also see mcextreg.bhdr. This is a binary file with the orthogonalized motion correction parameters which can later be used as nuisance regressors when analyzing the data. These files will exist in each of the runs (i.e., 005, 007, 010).
cd $FSFTUTDIR/fb1-raw-study ls mgh-101.1/bold/masks
You will see a brain.nii volume in the masks folder. This is a binary mask of the brain as found by FSL's BET program. The functional data themselves are not masked.
To view the translation components, run
cd $FSFTUTDIR/my-study plot-twf-sess -s mgh-101.1 -mc
This will bring up a plot with the translations for each of the runs. Notice that run 007 has a head jerk at about time point 11.
In order to render the functional results on the anatomical background as well as to map the functional results into a common space for group analysis, it is necessary to register/align the functional volume with a structural volume. In FSFAST, we:
First register to the same-subject FreeSurfer anatomical with a 6 DOF registration.
- Map the functional to Talairach/MNI305/fsaverage space by concatenating the within-subject function/structure registration with the Talairach registration (talairach.xfm) created when the subject was reconstructed.
- For mapping to surface-based space, we concatenate the within-subject function/structure registration with the surface-based registration.
Since we are only dealing with the functional analysis here, we will just consider the within-subject function/structure registration.
For these series of exercises, cd into fb1-preproc-study:
View unregistered (tkregister-sess)
Run the following command to see how the functional and structural are aligned prior to performing any automatic registration (delete the old one, if there).
cd $FSFTUTDIR/fb1-preproc-study rm -f mgh-101.1/bold/register.dat tkregister-sess -s mgh-101.1 -regheader -inorm
You should see the anatomical image (left, below). The green line is the cortical surface for that subject. Hit the compare button to see the functional middle image. The cortical surface shows that the volumes are not in line.
Run automatic registration (spmregister-sess)
cd $FSFTUTDIR/fb1-preproc-study fslregister-sess -s mgh-101.1
When this command is complete, you will see a register.dat file in mgh-101.1/bold. This is the only change. The functional data are not resampled! Note that this command requires FSL. You can also run:
cd $FSFTUTDIR/fb1-preproc-study spmregister-sess -s mgh-101.1
This does the same thing but requires matlab and spm in your matlab path (in startup.m), otherwise use fslregister-sess.
Check automatic registration (tkregister-sess)
cd $FSFTUTDIR/fb1-preproc-study tkregister-sess -s mgh-101.1 -inorm
When you hit the Compare button, you should see the image on the right, above.
Check talairach registration
When performing a volume-based group analysis, it is important that the talairach transform be accurate (or as accurate as talairach can be). You should have done this during the FreeSurfer reconstruction for this subject. The command below is just a reminder of how to do it.
tkregister2 --s fbph1-101 --fstal --surf
The First-Level Analysis (FLA) consists of:
- Setting up models of the task-related components
- Setting up models of the nuisance components
- Defining contrasts
- Fitting the model
- Making inferences.
In FSFAST, the FLA is done in two stages. In the first stage, the FLA is configured (with mkanalysis-sess). This is done ONCE regardless of how many data sets you have (you do not even need to have any data to run the configuration).
In the second stage, you actually perform the analysis/inference with selxavg3-sess by passing it the configuration and the session that you want to analyze. selxavg3-sess customizes the analysis for that session based on what it finds in the hierarchy, builds the design matrix and performs the analysis. Breaking the FLA up into these two stages assures that all sessions are analyzed in the same way.
In this exercise, we will configure two analyses, one assuming an HRF and the other using a Finite Impulse Response (FIR) model. You will then launch the analysis and view the results. The fully analyzed data can be found in fb1-analysis-study. These exercises will be done in fb1-preproc-study:
Configure Analysis and Contrasts I: Gamma HRF Model
Configuring the FLA is performed with mkanalysis-sess. When you run:
cd $FSFTUTDIR/fb1-preproc-study mkanalysis-sess -gui
You will see the following window:
You will use this window to specify the input of the analysis, the hemodynamic response model, contrasts and nuisance regressors. The red fields are fields that you must enter before you can save the analysis. There is a lot going on with this GUI, so we'll break it down. Note that many of the components have "tooltips" that will show when you pause the mouse pointer over them.
In the upper left corner is a panel called "FS-FAST Hierarchy". The "Func Stem" is the input to the analysis. You should specify the output from the preprocessing. For this exercise, we are going to use the motion corrected and 5mm smoothed data. This functional volume is called fmcsm5.nii in the hierarchy which makes its stem "fmcsm5" (i.e., just strip off the nii). Enter "fmcsm5" into the field. When you hit return, it changes from red to white. Next, enter the TR (sec). For this experiment it was 3 sec. This will be checked against the TR found in the input nifti file. Leave INorm checked.
Func Stem = fmcsm5 TR = 3
Turn your attention to the "Noise and Nuisance Variables" panel. Low frequency noise so prevalent in fMRI is compensated for in a combination of three ways. Drift components are modeled with polynomial regressors. The order can be adjusted, but leave it at 2 for now. The motion correction parameters can be used as regressors by checking the "MC Regressors" box (do so now). Finally, the remaining noise is modeled as time-invariant linear AR1 process when the "Temporal Whitening" box is checked (leave it so). There is one additionally way to compensate for noise through the use of a "Time Point Exclude File", but we will not consider that here.
Check "MC Regressors" Box
You will specify the model of the task-related signal in the "Event Related/Block Design" panel (leave that box checked). Choose the number of conditions by clicking on the "NConditions" slider. This is the number of TASK conditions (do not include the Null/Fixation condition). In this example, we have two conditions (Odd and Even), so set this to 2. To the right of this is the "Paradigm File". Enter "sensory-motor.par". Note that the number of task conditions in the paradigm file must match that specified with "NConditions". Below, you will specify the Hemodynamic Response Model. There are three choices: Gamma, SPM HRF, and FIR. Choose Gamma for now. If you hit the "Plot" button it will show the Gamma and SPM HRF. As you change the Gamma parameters (Delay, Dispersion (Tau), and Exponent (Alpha)), the Gamma plot will change. Make sure that they are at Delay=2.25, Tau=1.25, and Alpha=2.
Set NConditions to 2 Paradigm File = sensory-motor.par Press "Plot" button
At this point, you have specified the model of the BOLD signal including HRF, nuisance, and noise. The GUI should look like the image below:
Now you are ready to specify contrasts. A contrast is an instantiation of a hypothesis and is represented by a contrast matrix (i.e., a linear summation of the regression coefficients). Contrasts are managed through a separate GUI accessed through the "Contrast" list box. When you click on "Add Contrast", you will see the following screen:
Click on "Add Contrast"
There are several things going on here, but the most important is the list of conditions in the middle of the GUI (i.e., "Condition 1", "Condition 2") shown as green, red, and black radio buttons.
- Green indicates an "active" condition.
- Red means a "control" condition
- Black means to ignore the condition in the contrast.
Active conditions are given a weight of +1; controls are given -1; ignores get 0. The weight is given to the right of the buttons. All contrasts are implicitly computed against the Null or Fixation condition. If you want to test the null hypothesis that Condition 1 is no different than the Null condition, then you would make Condition 1 active and ignore the rest. To test the null hypothesis that Condition 1 is no different than Condition 2, then you would make Condition 1 active and Condition 2 control.
For this exercise, we are going to test four NULL hypotheses:
- Odd == Fixation (odd-v-0)
- Even == Fixation (even-v-0)
- Odd == Even (odd-v-even)
- Odd+Even == Fixation (odd+even)
The last one tests whether the average of the responses to odd and even are different than fixation. Remember that, according to the
Edit conflict - other version:
Paradigm File, Condition 1 is Odd, and Condition 2 is Even.
Click on the green button next to "Condition 1". (see its weight change from 0 to 1).
Condition 2 is already black (ignored). This our first contrast, so there is nothing more we need to do except give the contrast a name. You should give your contrasts meaningful but terse names.
Specify "odd-v-0" for the name of this contrast. Hit the "Done/Save" button.
You will now see "odd-v-0" appear in the Contrast list box in the mkanalyiss GUI.
Click on "Add Contrast" again to bring up the contrast GUI again. Click on the green button next to Condition 2 (see its weight change from 0 to 1).
This activates condition 2. Condition 1 is still black (ignored).
Change the name to "even-v-0" Click Done/Save
"even-v-0" will appear in the list box.
Click on "Add Contrast" again Click the red button next to Condition 2 (see its weight change from 0 to -1). Change the name to "odd-v-even" Click Done/Save.
Click on "Add Contrast" one more time Click the green button next to Condition 1. lick the green button next to Condition 2. Change the name to "odd+even" Click Done/Save.
Edit conflict - your version:
Paradigm File, Condition 1 is Odd, and Condition 2 is Even. When "Add Contrast" is clicked, "Condition 1" will be active and Condition 2 will be ignored. This corresponds to our first contrast, so there is nothing we need to do except give the contrast a name. You should give your contrasts meaningful but terse names. Specify "odd-v-0" for this contrast. Hit the "Done/Save" button. You will now see "odd-v-0" appear in the Contrast list box in the mkanalyiss GUI.
Click on "Add Contrast" again to bring up the contrast GUI again. This time, click on the green button next to Condition 2 (see its weight change from 0 to 1). Then click on the black button next to Condition 1 (see weight change from 1 to 0). Change the name to "even-v-0", then click Done/Save. "even-v-0" will appear in the list box.
Click on "Add Contrast" again, and click the red button next to Condition 2 (see its weight change from 0 to -1). Change the name to "odd-v-even", then click Done/Save.
Click on "Add Contrast" one more time, and click the green button next to Condition 1 and Condition 2 (see its weight change from 0 to +1). Change the name to "odd+even", then click Done/Save.
End of edit conflict
You can go back and view and/or edit an contrast by clicking on it in the list box.
The last thing you have to do is to give your analysis a name. Like the contrasts, it should be terse but descriptive. It cannot have any
Edit conflict - other version:
spaces or blanks.
Specify "sm-gamma-fwhm5" as the contrast name.
sm = sensory-motor, gamma = Gamma HRF, and fwhm5 for the input. The interface should now look
Edit conflict - your version:
spaces or blanks. Specify "sm-gamma-fwhm5" (sm = sensory-motor, gamma = Gamma HRF, and fwhm5 for the input). The interface should now look like:
End of edit conflict
Hit the "Save" button, then "Quit".
After you hit Quit, control will be returned to the shell that you ran mkanalysis-sess from. If you type "ls", you will see a new folder called "sm-gamma-fwhm5". If you "ls sm-gamma-fwhm5", you will see analysis.info, analysis.cfg, odd-v-0.mat, even-v-0.mat, odd-v-even.mat, odd+even.mat. Your configuration is stored in these files. You can browse/edit your configuration by running:
mkanalysis-sess -gui -analysis sm-gamma-fwhm5
Configure Analysis and Contrasts II: FIR HRF Model
Now we are going to use a Finit Impulse Response (FIR) to model the hemodynamic response. The FIR does not make any assumptions about the shape of the HRF but is also less interpretable. Again, run
Set the Func Stem, TR, NConditions, and Paradigm File as above, but now click on the "FIR" checkbox. This will enable the "Total Time Window", "PreStim", and "TER" entry boxes. The Time Window is the window within which we will estimate the HRF. Given that the task is 15 sec long and the rest is 15 sec, let's choose 27 sec. The PreStim is the amount of time before stimulus onset to start estimating the HRF. A non-zero PreStim gives us an idea of what the baseline is at stimulus onset. Set it to 6.
Set up the same contrasts as you did above, then name the analysis "sm-fir-fwhm5", hit Save, then Quit.
Analyze First Level (selxavg3-sess)
You are now ready to analyze some data! Note that the fully analyzed data (along with correctly configured analyses) can be found in fb1-analysis-study. To analyze the data for session mgh-101.1 with the sm-gamma-fwhm5 analysis, run:
selxavg3-sess -s mgh-101.1 -analysis sm-gamma-fwhm5
Note that if you want to analyze all the sessions, you can run "selxavg3-sess -sf sessid -analysis sm-gamma-fwhm5", but this might take a while.
To analyze the data for session mgh-101.1 with the sm-fir-fwhm5 analysis, run:
selxavg3-sess -s mgh-101.1 -analysis sm-fir-fwhm5
Examine additions to the FSFAST hierarchy
ls mgh-101.1/bold ls mgh-101.1/bold/sm-gamma-fwhm5
The bold folder now has sm-gamma-fwhm5 (as well as the fir too). This folder has all of the results for this analysis, including:
- beta.nii - regression coefficients
- rvar.nii - residual error variance
- mask.nii - mask (copy of bold/masks/brain.nii)
- meanfunc.nii - mean functional image
- fsnr.nii - functional SNR map
- X.mat - design matrix (in matlab format)
- dof - text file with degrees of freedom
- fwhm.dat - text file with smoothness estimate (Full-Width/Half-Max)
- contrast folders: odd-v-0, even-v-0, odd+even, odd-v-even
Look in one of the contrasts:
- ces.nii - contrast effect size (contrast matrix * regression coef)
- cesvar.nii - variance of CES
- sig.nii - significance map (-log10(p))
Visualization of the first-level analysis will be done in the folder with the fully analyzed data set:
Volume-based visualization (tkmedit-sess)
View the result of the Gamma HRF analysis on the FreeSurfer anatomical volume for mgh-101.1 with:
tkmedit-sess -s mgh-101.1 -aparc+aseg -analysis sm-gamma-fwhm5 \ -c odd-v-0 -c even-v-0 -c odd+even -c odd-v-even
This will automatically load sig.nii. When you configure the functional overlay with View->Configure->Functional Overlay, you will see that there are 4 "Time Points". Each point is different contrast (i.e., 0 is odd-v-0, 1 is even-v-0, etc). Scroll through each one.
View the result of the HRF HRF analysis for mgh-101.1 with:
tkmedit-sess -s mgh-101.1 -aparc+aseg -analysis sm-fir-fwhm5 -c odd-v-0
Here we view only one contrast at a time because each contrast has multiple time points for each point in the FIR time window. When you click on a point, you will see the HRF for both Condition 1 (Odd) and Condition 2 (Even) blocks.
Finally, view the overlay maps from the Gamma with the HRF from the FIR:
tkmedit-sess -s mgh-101.1 -aparc+aseg -analysis sm-fir-fwhm5 \ -mapanalysis sm-gamma-fwhm5 -c odd-v-0 -c even-v-0 -c odd+even -c odd-v-even
When you configure the overlay, you will see that there are 4 "Time Points" -- these correspond to the 4 contrasts from the Gamma analysis.
Surface-based visualization (tksurfer-sess)
View the result of the Gamma HRF analysis on the FreeSurfer anatomical surface for mgh-101.1 with:
tksurfer-sess -hemi lh -aparc -s mgh-101.1 -analysis sm-gamma-fwhm5 \ -c odd-v-0 -c even-v-0 -c odd+even -c odd-v-even
Note that the command-line is nearly identical to that of tkmedit above. The difference is that the hemisphere is specified ("-hemi lh"), and "-aparc+aseg" is replaced with "-aparc" to load the surface-based segmentation.
Higher-Level (Group) Analysis
Higher-Level is where you make inferences about the population that your subjects are drawn from. It is a bit confusing at times because both use GLMs, so at both levels you are constructing design matrices, contrasts, etc.
Traditionally, fMRI group analysis has been done in a standard volume space (i.e., Talairach/MNI152/MNI305). With FreeSurfer, we also have the option to analyze group data in the surface space. Volume-based analyses are done in MNI305 space (which is the same as the fsaverage subject).
For the next exercises, we will work in the fb1-analysis-study directory
All the first level analyses have been done here. The group-level analyses have been done in the group-analysis-tut directory.
Volume-based (MNI305/fsaverage) Group Analysis
Assemble the Data (isxconcat-sess) (Volume)
The first step in the group analysis is to "assemble" the data. This means creating a single 4D file with where the 4th "time" dimension is actual all the subjects concatenated together in a common space. There is a different command, depending upon whether the common space is volume- or surface-based.
To run the volume-based concatenation, run the command below. Note that the data from this command already exist in the group-analysis-tut directory.
isxconcat-sess -sf sessid -analysis sm-gamma-fwhm5 -c odd-v-0 -o group-analysis
This command will:
- Go through each session in the sessid file
- Find the odd-v-0 contrast in the sm-gamma-fwhm5 analysis
- Use the register.dat for that session to resample to MNI305/fsaverage space.
- Concatenate all subjects' data into one file
The concatenated data are saved in group-analysis/sm-gamma-fwhm5/odd-v-0/tal.ces.nii file. You can run other contrasts by adding -c arguments. In addition to this output, several other files are created. To see them,
cd $FSFTUTDIR/fb1-analysis-study ls group-analysis-tut/sm-gamma-fwhm5 cat group-analysis-tut/sm-gamma-fwhm5/sessid cat group-analysis-tut/sm-gamma-fwhm5/ffxdof.dat
In the output directory, you will see a series of files that start with "tal":
- tal.meanfunc.nii is a stack where each "time point" is the mean functional image of each subject sampled in the MNI305 space.
- tal.masks.nii are the binary masks for all the subjects
- tal.fsnr.nii are the functional SNR maps from each subject.
- tal.mask.nii is a single binary mask made from the intersection of the individuals.
- ffxdof.dat is the fixed-effects DOF across all subjects.
- sessid.txt is the list of sessions, the corresponding freesurfer subject name, and the DOF contributed by each subject.
You will also see some files that being with "lh". These are the same thing in the surface-based space (see below).
Now look in the directory for odd-v-0 the contrast
cd $FSFTUTDIR/fb1-analysis-study ls group-analysis-tut/sm-gamma-fwhm5/odd-v-0
You will see:
- tal.ces.nii - the contrast maps for each of the subjects
- tal.cesvar.nii are the variance of the contrast for each subject (i.e., the square of the standard error). This variance is needed for fixed-effects and weighted random-effects analysis.
Quality Assurance (Volume)
There are three important quality assurance steps that can be performed here. First, view the mean functionals to make sure that all are registered together properly. To do this run,
tkregister2 --s fsaverage --surf --regheader --check-reg \ --mov group-analysis-tut/sm-gamma-fwhm5/tal.meanfunc.nii
The image window will show the MNI305 brain. Hit the "Compare" button to show the average functional of the first session. Click in the image window, then hit the 'a' key. Each time you hit the 'a' key, it will advance to the next subject.
The next QA step is to check the individual masks. This can be done with:
tkmedit fsaverage orig.mgz \ -overlay group-analysis-tut/sm-gamma-fwhm5/tal.masks.nii -fthresh 0.5
The threshold of 0.5 is appropriate because these masks are binary (ie 0-1). When you View->Configure->Functional Overlay, you will see that there are 5 "Time Points" (0-4) corresponding to the 5 subjects. Advance through each one to assure that the masks are in the proper place.
The final step is to look at the functional SNR maps with
tkmedit fsaverage orig.mgz \ -overlay group-analysis-tut/sm-gamma-fwhm5/tal.fsnr.nii \ -timecourse group-analysis-tut/sm-gamma-fwhm5/tal.fsnr.nii \ -fthresh 50 -fmax 250
When you View->Configure->Functional Overlay, you will see that there are 5 "Time Points" (0-4) corresponding to the 5 subjects. Advance through each one to assure that the SNR maps are "consistent".
When you click on a voxel, you will see the SNR for each subject plotted in the "Time Course" window. The actual value of the FSNR will vary depending upon how much smoothing you've done and the details of the acquisition. You are looking for outliers here.
Group GLM Analysis (Volume)
When you perform a group analysis, you are looking for effects of the task across the population. In this example, we are only going to test whether the mean effect of Odd vs Fixation (odd-v-0) is different than 0. This is known as a one-sample-group-mean (OSGM) and corresponds to a group design matrix that is simply a column of 1s. One can perform considerably more elaborate tests (e.g., differences in groups). Now, cd to where the data are:
Random Effects (RFx, OLS) (Volume)
Random effects is a test of whether the mean of the population that the subjects were drawn from is 0. It is done with mri_glmfit:
cd $FSFTUTDIR/fb1-analysis-study/group-analysis-tut/sm-gamma-fwhm5/odd-v-0 mri_glmfit --y tal.ces.nii --osgm --mask ../tal.mask.nii --glmdir tal.rfx.osgm --nii
- tal.ces.nii are the inputs for this contrast for all subjects.
- --osgm tells it to create the simple OSGM design matrix and to create a group contrast to test the group mean against 0.
- The mask is used to exclude extra-brain voxels.
The results will be stored in tal.rfx.osgm:
cd $FSFTUTDIR/fb1-analysis-study/group-analysis-tut/sm-gamma-fwhm5/odd-v-0 ls tal.rfx.osgm ls tal.rfx.osgm/osgm
You will see several files in the output. beta.nii is the map of regression coefficients. There is also an "osgm" folder with sig.nii. This is a map of the significances of the test. We will view it below.
Weighted Random Effects (WRFx, WLS) (Volume)
In a real experiment, some subjects are noisier than others, and it is a good idea to take this into account since we have information about how noisy a subject is through the lower-level analysis. In weighted least squares (WLS), this is handled by weighting each subject by the inverse of their noise (i.e., noisier subjects get lower weight). To do this with mri_glmfit, run:
cd $FSFTUTDIR/fb1-analysis-study/group-analysis-tut/sm-gamma-fwhm5/odd-v-0 mri_glmfit --y tal.ces.nii --osgm --glmdir tal.wrfx.osgm --nii --mask ../tal.mask.nii \ --wls tal.cesvar.nii
The command-line is almost the same as the RFx except that the first-level noise variances (tal.cesvar.nii) are passed with the --wls option. The results are saved in tal.wrfx.osgm and have the same naming convention as the RFx.
Fixed Effects (FFx) (Volume)
A fixed effects analysis tests whether an effect exists within the given subject pool (as opposed to the population that the subjects were pulled from). It is like doing one big fist-level analysis. The big difference with RFx or WRFx is that the variance is computed from the lower level analysis variance. This is done in mri_glmfit with:
cd $FSFTUTDIR/fb1-analysis-study/group-analysis-tut/sm-gamma-fwhm5/odd-v-0 mri_glmfit --y tal.ces.nii --osgm --glmdir tal.ffx.osgm --nii --mask ../tal.mask.nii \ --yffxvar tal.cesvar.nii --ffxdofdat ../ffxdof.dat
Again, the command-line is similar to above, except that the lower level variances are passed with --yffxvar, and the total lower level DOFs are passed with --ffxdofdat.
Output and visualization (Volume)
To visualize these three analyses, we will first concatenate them into one file to make it easier to view in tkmedit:
cd $FSFTUTDIR/fb1-analysis-study/group-analysis-tut/sm-gamma-fwhm5/odd-v-0 mri_concat tal.rfx.osgm/osgm/sig.nii \ tal.wrfx.osgm/osgm/sig.nii \ tal.ffx.osgm/osgm/sig.nii \ --o all.sig.nii) tkmedit fsaverage orig.mgz -aux brain.mgz -bc-main-fsavg \ -overlay all.sig.nii -fthresh 2 -fmax 4
Configure the overlay with View->Configure->Functional Overlay. Scroll through the "time" points to see the differences in the analyses.
Correction for Multiple Comparisons/Cluster Analysis (Volume)
With so many voxels in fMRI maps, it is very likely that many voxels will appear to be active purely by random chance (ie, a false positive). The is known as the "Problem of Multiple Comparisons". One way around this is to do a cluster analysis in which active voxels are eliminated unless they appear in a cluster, the idea being that false positives will not appear next to each other. To do this in FSFAST, run the following command on the RFx analysis:
mri_volcluster --in tal.rfx.osgm/osgm/sig.nii --thmin 2 \ --fwhmdat tal.rfx.osgm/fwhm.dat --cwpvalthresh 0.1 \ --mask tal.rfx.osgm/mask.nii --fsaverage \ --sum tal.rfx.osgm/osgm/cluster.sum --sign pos\ --cwsig tal.rfx.osgm/osgm/cwsig.cluster.nii \ --out tal.rfx.osgm/osgm/sig.cluster.nii \ --ocn tal.rfx.osgm/osgm/ocn.cluster.nii
In this command we use the significance map as input thresholded at 2, where 2 is -log10(p), so p < .01. The --fwhmdat option tells mri_volcluster to use the spatial smoothness estimate from GLM analysis (saved in tal.rfx.osgm/fwhm.dat). The "--cwpvalthresh 0.1" tells it to ignore clusters unless their p-value is less than 0.1 (you can threshold more stringently later on).
One of the outputs is a cluster table (cluster.sum). This is a text file with the list of clusters. View it with:
You can also see it here: FsFastTutorial/ClusterSummary. The table shows the size of each cluster in voxels and mm^3, the Talairach coordinate, the maximum significance in the cluster, and the clusterwise p-value (CWP). Some of them are 0, meaning that the p-value was so low that it could not be computed (remember, this is a good thing!).
The other outputs are volumes. cwsig.cluster.nii is a map of the clusters with the voxel value equal to the -log10(pvalue) of the cluster that the voxel is associated with. sig.cluster.nii is the original sig map with non-cluster voxels removed. ocn.cluster.nii is a map where the value at each voxel is replaced by the number of the cluster that the voxel is associated with. View the clusterwise significance:
tkmedit fsaverage orig.mgz -aux brain.mgz -bc-main-fsavg \ -overlay tal.rfx.osgm/osgm/cwsig.cluster.nii -fthresh 2 -fmax 4
Surface-based Group Analysis
Assemble the Data (isxconcat-sess) (Surface)
The surface-based analysis proceeds in much the same way as the volume-based. The assembly command is the same except you add "-hemi lh" (or "-hemi rh") to the command:
cd $FSFTUTDIR/fb1-analysis-study isxconcat-sess -sf sessid -analysis sm-gamma-fwhm5 -c odd-v-0 -o group-analysis -hemi lh
This creates the same files as the volume-based command, except that they have "lh" or "rh" preprepended instead of "tal".
ls group-analysis/sm-gamma-fwhm5 ls group-analysis/sm-gamma-fwhm5/odd-v-0
You should see the following files: lh.meanfunc.nii, lh.mask.nii, lh.fsnr.nii, lh.ces.nii, lh.cesvar.nii.
Group GLM Analysis (Surface)
The commands for the GLM analysis are the same as for the volume-based except that "--surf fsaverage lh" is added to the command line:
Cd into the directory with the assembled data:
Random Effects (RFx, OLS) (Surface)
mri_glmfit --y lh.ces.nii --osgm --mask ../lh.mask.nii --glmdir lh.rfx.osgm --nii \ --surf fsaverage lh
Weighted Random Effects (WRFx, WLS) (Surface)
mri_glmfit --y lh.ces.nii --osgm --glmdir lh.wrfx.osgm --nii --mask ../lh.mask.nii \ --wls lh.cesvar.nii --surf fsaverage lh
Fixed Effects (FFx) (Surface)
mri_glmfit --y lh.ces.nii --osgm --glmdir lh.ffx.osgm --nii --mask ../lh.mask.nii \ --yffxvar lh.cesvar.nii --ffxdofdat ../ffxdof.dat --surf fsaverage lh
Output and visualization (Surface)
For convenience, we concatenate the different analyses into a single file, then visualize with tksurfer:
mri_concat lh.rfx.osgm/osgm/sig.nii lh.wrfx.osgm/osgm/sig.nii lh.ffx.osgm/osgm/sig.nii \ --o lh.all.sig.nii tksurfer fsaverage lh inflated -annot aparc -overlay lh.all.sig.nii -fthresh 2
Correction for Multiple Comparisons/Cluster Analysis (Surface)
In principle, correction for multiple comparisons on the surface is done in the same way as volume-based, but in implementation is considerably more difficult. We have several methods for correcting on the surface, including False Discovery Rate (FDR), Monte-Carlo Simulation, and Permutation Simulation. The Guassian Random Field (GRF) is still being developed (as are the tutorials for it).
Play with Your Data
In this section, you are encouraged to experiment with the data. These analyses have not already been pre-computed, but they should not take long to run.
Change the smoothing level
Use this series of commands to analyze the data smoothed by 10 mm FWHM. The preproc-sess command will create fmcsm10.nii. In mkanalysis-sess, use the -clone option to copy the parameters from an existing analysis, then just change the Func Stem to fmcsm10.
cd $FSFUTDIR/fb1-analysis-study preproc-sess -s mgh-101.1 -fwhm 10 mkanalysis-sess -gui -clone sm-gamma-fwhm5 -analysis sm-gamma-fwhm10 selxavg3-sess -s mgh-101.1 -a sm-gamma-fwhm10 tkmedit-sess -s mgh-101.1 -analysis sm-gamma-fwhm10 \ -c odd-v-0 -c even-v-0 -c odd+even -c odd-v-even -aparc+aseg
Change noise/nuisance model
Use this series of commands to analyze the data with little compensation for noise and nuisance variables. Change the analysis by turning off whitening, set the Polynomial Order to 0, and uncheck the MC Regressors.
cd $FSFUTDIR/fb1-analysis-study mkanalysis-sess -gui -clone sm-gamma-fwhm5 -analysis sm-gamma-fwhm5B selxavg3-sess -s mgh-101.1 -a sm-gamma-fwhm5B tkmedit-sess -s mgh-101.1 -analysis sm-gamma-fwhm5B \ -c odd-v-0 -c even-v-0 -c odd+even -c odd-v-even -aparc+aseg