Freesurfer BOLD files
Freesurfer, like SPM, can perform functional analyses on fMRI data. These analyses generally employ the general linear model (GLM) to assess the degree to which the blood-oxygen dependent (BOLD) signal correlates with different predictors (e.g., group, experimental condition, task or even nuisance regressors such as motion).
For Freesurfer and SPM, these files will start off as NIfTI (.nii) files. Each file contains a header, which contains important information directing the software how to interpret the data, and one or more 3D matrices of raw data, representing the signal at each voxel at a particular point in time, as measured by the scanner during the experiment. The data for a single run (e.g., a 5-minute scan of an experimental task) can be stored as a series of individual 3D volume files, with each file representing a single scan of the brain like a single frame of a 5-minute movie. Alternately, multiple volumes can be concatenated together and stored in a single 4D file. Whereas SPM operates on either 3D or 4D data, Freesurfer seems to require 4D data.
The header information is important because it provides the context required to interpret the data matrix. For example, an experiment looking at hippocampal processing might use a 32 x 32 x 24 array of 1mm voxels and look at a relatively small subcortical volume. Another experiment looking at global cortical processing might use a 32 x 32 x 24 array of 3mm voxels to look at the entire brain volume. Without knowing the voxel size, it would be impossible for any piece of software to interpret the data matrix. The NIfTI header includes information about how big the voxels are, the orientation, and many other important bits of information, though not all software makes use of all the header information.
File Fixing
Freesurfer expects the BOLD data to be organized into individual folders for each experimental run, and that each file is called "f.nii". Moreover, Freesurfer (but not SPM) appears to make use of the TR value stored in the NIfTI header. In order to ensure that the filename requirement is met and that the TR value is correctly set, a script called set_TR.m can be found in the ubfs Scripts/Matlab folder. Another issue, specific to our scanner, is the fact that the Toshiba scanner at the CTRC does not (currently) do anything to account for the magnetic field inhomogeneity that invariably accompanies the first several seconds of a functional run. Other scanners (e.g., Siemens) automatically discard the first several seconds of data. Because the BOLD data from the CTRC will have this magnetic field inhomogeneity artifact, we must manually drop the affected volumes and also remember to adjust the onsets of the experimental events accordingly.
Golden Rule
Never directly modify the BOLD data in the raw files. Instead, copy the raw data to your working directory and work on copies of the data. This avoids the potentially catastrophic possibility of procedures being applied multiple times to the same data, resulting in data loss or mangling.
Note that certain operations are fine to apply to BOLD data because they don't change the data itself, but instead just make changes to the NII file headers that are easily detected and can easily be reversed. The only allowable change that you might consider doing to the raw data, however, is to set the TR.
Determining which volumes are important
Before you drop any volumes (the next step), you should first check to see which volumes were associated with any events. A MATLAB function called critVols can be found in the ubfs/Scripts/Matlab folder. Copy it to your MATLAB path if you haven't already done so. This function requires the participant's PsychToolBox MATLAB runtime files -- the same files that are used to generate .par files. The simplest way to run the script is to call it with no parameters:
v=critVols(); %assume the default TR of 2.047
The function assumes a TR of 2.047, which is correct for the LDT and Multisensory Imagery study. If you are using it for any other study, you may need to specify the TR:
v=critVols('TR', 2.00); %overrides the default TR
As usual, there is a built-in help
help critVols
After you run this function on a set of .mat files, you will get a list of start and end volumes for each of the associated runs:
v{:}
Below is the information for the 6th run:
ans =
run: 6
subject: '0202'
startvol: 5
endvol: 242
This would tell us that nothing happened during volumes 1:4 and any volume from 242:end. You will almost certainly want to drop the first 4 volumes. Depending on your requirements, you may wish to drop the volumes at the end as well.
Fixing Files in the Terminal
It turns out there are several useful utilities associated with AFNI and FSL that can let us do everything we need to do in the terminal, without running MATLAB. This is helpful if you want to automate the process with a script.
Setting the TR in the Terminal Using 3drefit
The AFNI utility, 3drefit can do all sorts of things, including the trivial task of modifying the TR field in your BOLD data when you use the -TR switch:
FILENAME=f.nii #the name of the file to modify
TR_DURATION=2.047 #the duration of the TR
3drefit -TR ${TR_DURATION} ${FILENAME}
Running the above command will produce output like the following:
++ 3drefit: AFNI version=AFNI_2011_12_21_1014 (Jan 8 2015) [64-bit] ++ Authored by: RW Cox ++ Processing AFNI dataset f.nii + loading and re-writing dataset f.nii (./f.nii in NIFTI storage) ++ 3drefit processed 1 datasets
Dropping Volumes in the Terminal Using fslsplit and fslmerge
The FSL utility fslsplit will break up a 4D BOLD volume into a series of 3D BOLD volumes. The default behavior is to produce files named vol0000.nii.gz to voln-1.nii.gz from a 4D volume with n time points. For example, if your BOLD volume was 220 time points, this command would produce vol0000.nii.gz, vol0001.nii.gz, ... , vol0219.nii.gz
Conversely, the FSL utility fslmerge does the reverse.
We can chain these two commands together, deleting any unwanted volumes in between. Suppose we wish to drop the first 4 volumes. These correspond to vol0000.nii.gz, vol0001.nii.gz, vol0002.nii.gz and vol0003.nii.gz. Fortunately, the BASH shell uses regular expressions which have a convenient shorthand for specifying wildcard characters in a range of values. In this case, we will want to delete all the files that start with "vol000" followed by the numbers 0 to 3, followed by ".nii.gz". This works out to the following pattern:
vol000[0-3].nii.gz
Now we can put this all together:
BASEFILENAME=f.nii #the file we will drop volumes from
DELETE_TO=3 #we wish to drop the first 4 volumes, but we start counting at 0, so we stop at 3 (0,1,2,3) rather than 4
mkdir old #a directory to store the backup out of the way
cp ${BASEFILENAME} old/ #back up the original file just in case
fslsplit ${BASEFILENAME} #break up the original file
rm ${BASEFILENAME} #delete the original file
rm vol000[0-${DELETE_TO}].nii.gz #delete the files corresponding to first 4 volumes
#merge all the remaining vol*.nii.gz files into a single 4D volume
#it will gzip the file, producing f.nii.gz. That's ok if it's gzipped.
fslmerge -t f.nii vol*
rm vol* #delete all the individual 3D BOLD volumes
#inspect the new f.nii.gz file using fslhd to confirm the header information is what you expect. If so, delete the backup
rm -rf old
Fixing Files in MATLAB
Originally, we were using a series of MATLAB scripts to modify the TR and drop volumes from our BOLD datasets. The sections below describe this process for posterity. These methods will still work; they just don't lend themselves as well to command-line scripting.
Dropping initial volumes using nii_4D_drop_vols.m
A MATLAB function, nii_4D_drop_vols.m can be found along with the rest of our library of MATLAB code. This code relies on Jimmy Shen's NiFTI Tools toolbox. Ensure that this function and the toolbox are in your MATLAB path. If they are, you can type help nii_4D_drop_vols.m, which will tell you that you can call the function with or without a filename (you will be prompted to find a file if none is provided). The second parameter is the list of volumes to drop. Positive values are the volumes, relative to the start of the series, to drop. Values 0 or less are the volumes from the end of the series. Your cwd should be an individual bold folder (e.g. ~/LDT/FS_1234/bold/001). This means you have to edit each *.nii file individually. Shell script should be coming soon to alleviate this annoyance.
nii_4D_drop_vols('./<name of par files>', [<number of initial drop vols> <number of end vols to drop>]);
- note: numbers should be separated by spaces; for the end vols, 0 is the last, -1 is the second to last, etc
For example:
nii_4D_drop_vols('./4D.nii', [1 2 3 4 0 -1 -2]);
%will drop volumes 1 to 4, as well as the last (0), second-last (-1) and third-last (-2)
nii_4D_drop_vols([], 1:4);
%will drop volumes 1 to 4, but you will be prompted to find a file because an empty array was given as a filename
The original file will be saved with a .bak extension.
If this function doesn't work, the first thing you should check is that the NiFTI tools toolbox folder is in your path. Also note that the function will be looking for a fully qualified filename (including the path). If you get an error message indicating that it "Cannot CD to (No directory specified)", this is because you didn't indicate the full path to the target file. The example above uses ./4D.nii, which specifies that the path to the file is the current working directory. Try that out if you're having a problem (if you do, you must run the function from the directory containing the target file).
Dropping the first 8 seconds should suffice. With a TR of just over 2 seconds, this would be the first 4 volumes (for the sake of consistency across- and within-projects, let's assume that we will always drop 4 volumes from our CTRC data unless you are explicitly aware of a plan to do otherwise). The functionality of dropping the end volumes was added because it was trivial to do so.
Remember to account for the dropped volumes when creating your .par files
Combining critVols with nii_4D_drop_vols
Specifying the initial volumes to drop is easy, as it is almost always just the first 4 volumes, 1:4 (1 up to, but not including, the first critical volume). This is usually the case, but make sure you check the startvol from critVols and drop one less than that (for example, if it starts at volume 5, drop 4 volumes). The end volumes might be a bit trickier. Finding the total volumes of each f.nii file is explained above.
Let's make a variable called tail that lists the indices of the tail end of the run:
nvols=256; %determined using fslhd endvol=242+4; %retains 4 volumes after the last event occurs tail=(endvol+1:nvols)-endvol;
initial_drop=5-1 %subtract 1 from the startvol
nii_4D_drop_vols('./f.nii', [1:initial_drop, tail]);
It is not necessary, but I (rebecca) like to make a little matrix in matlab specifying all these data points. I'll make a column for the bold file numbers, one for start vol, endvol, nvol, and use this. Keeps things more organized and prevents you from having to reference crtiVols and then a reference sheet or fslhd for the nvol every time.
Fixing the TR header field using set_TR.m
The first requirement is that the source files for a participant be organized into numbered folders as described in the section below concerning file organization. If this requirement is met, you can use the MATLAB set_TR.m function to simultaneously rename the source files and set the TR to the correct value (you will, of course, have to know what the correct TR value is!). Data coming from the CTRC currently have the TR header field intialized to a value of 1.0, corresponding to a 1-second TR. SPM doesn't make use of the TR value from the header field, so this isn't a problem. However Freesurfer does use this value, so you may want to ensure that the header contains the correct TR value.
If you wish to simply inspect the values in a .nii file header, you can use a command-line FSL utility called fslhd:
fsledithd 4D.nii
This will echo header information to the screen
filename f.nii.gz sizeof_hdr 348 data_type FLOAT32 dim0 4 <--- 4 Dimensional BOLD dim1 64 <---64 voxels in the X (L/R) direction dim2 64 <---64 voxels in the Y (Front/Back) direction dim3 29 <---29 slices (the Z direction) dim4 216 <---216 time points dim5 1 dim6 1 dim7 1 vox_units mm <--voxel dimensions expressed in millimeters time_units s <---time units expressed in seconds datatype 16 nbyper 4 bitpix 32 pixdim0 0.000000 pixdim1 3.312500 pixdim2 3.312500 pixdim3 4.100000 <-- voxels are 3.3 * 3.3 * 4.1 pixdim4 2.047000 <-- time points are 2.047 seconds in duration pixdim5 0.000000 pixdim6 0.000000 pixdim7 0.000000 vox_offset 352 cal_max 0.0000 cal_min 0.0000 scl_slope 1.000000 scl_inter 0.000000 phase_dim 0 freq_dim 0 slice_dim 0 slice_name Unknown slice_code 0 <-- valid values are 0,1,2,3,4 0="unknown". For slice timing slice_start 0 slice_end 0 slice_duration 0.000000 time_offset 0.000000 intent Unknown intent_code 0 intent_name intent_p1 0.000000 intent_p2 0.000000 intent_p3 0.000000 qform_name Scanner Anat qform_code 1 qto_xyz:1 3.312500 0.000000 0.000000 -106.000000 qto_xyz:2 0.000000 3.312500 0.000000 -142.328003 qto_xyz:3 0.000000 0.000000 4.100000 -27.664000 qto_xyz:4 0.000000 0.000000 0.000000 1.000000 qform_xorient Left-to-Right qform_yorient Posterior-to-Anterior qform_zorient Inferior-to-Superior sform_name Scanner Anat sform_code 1 sto_xyz:1 3.312500 0.000000 0.000000 -106.000000 sto_xyz:2 0.000000 3.312500 0.000000 -142.328003 sto_xyz:3 0.000000 0.000000 4.100000 -27.664000 sto_xyz:4 0.000000 0.000000 0.000000 1.000000 sform_xorient Left-to-Right sform_yorient Posterior-to-Anterior sform_zorient Inferior-to-Superior file_type NIFTI-1+ file_code 1 descrip FSL5.0 aux_file
The above output gives information about the voxel dimensions in the x,y and z axes (dx, dy, dz are 3.3125, 3.3125 and 4.1 units, respectively). Those units are interpreted as millimeters, according to the xyz_units field value (2=mm). Similary, for the time dimension field, the value is set to 1 in this header, interpreted as seconds according to the time_units field (8=seconds; 16=milliseconds). Thus, according to this file header, the functional volumes were collected 1-second apart, which is not correct (the TR for this experiment was 2.047).
These values can be edited using fsledithd by simply changing the dt field value from 1 to 2.047 (if you do so, you would then type ^O (ctrl-O) to save the changes, then ^X (ctrl-X) to exit the editor). However, a MATLAB script exists that will update the fields for multiple files belonging to the same individual.
This script, set_TR.m can be found in the ubfs Scripts/Matlab folder. Ensure that this file can be found somewhere in your MATLAB path (e.g., copy it to your Documents/MATLAB folder)
First, open a terminal and start matlab:
matlab &
Your startup.m script should run automatically, setting your path to include the most commonly used toolboxes in the process. Nonetheless, you should ensure that set_TR and other needed functions are indeed in your path:
help load_untouch_nii help set_TR
If either of these commands return a command not found type of message, ask for help in getting them permanently in your path. It could be simple matter of running a startup script (e.g., startup_spm8).
Assuming that these files are in your path, the next step is to use MATLAB to navigate to the directory containing your participant data. For example, working with a data set on ubfs:
cd ~/ubfs/cpmcnorg/openfmri/test/1001
Now you are in the directory for participant 1001. There should be a bold subdirectory with one or more numbered directories below it, as described above. Move into the bold directory:
cd bold
Now you can run the set_TR command:
set_TR(2.047, [1 2 5 6], '4D.nii')
The above command would look for a file called '4D.nii' in directories 001, 002, 005 and 006, modify the NIfTI header file information to have a TR of 2.047 seconds, and then save the changed files to 'f.nii', leaving the originals unchanged (if your files are already called f.nii, they will be overwritten).
File Organization
All FreeSurfer BOLD files for an individual need to have the same name, and be organized in numerically labelled directories. As the example below illustrates, it is not imperative to our usual workflow that these folders be sequentially numbered (e.g., if there is any missing data). By convention, we store the BOLD data in the bold subfolder in the participant's FreeSurfer folder:
- $SUBJECTS_DIR
- $SUBJECT_ID
- bem
- bold
- 001
- f.nii
- 002
- f.nii
- 006
- f.nii
- 001
- label
- mri
- scripts
- src
- stats
- surf
- tmp
- touch
- trash
- $SUBJECT_ID
Handling Multiple Sessions
When data are acquired over multiple sessions, it is likely that the files will be in very different space. If you plan to analyze all the data simultaneously, you will need to ensure that all your BOLD files are coregistered with each other.
Is There a Coregistration Issue?
View a template BOLD file from each of your sessions, with the surfaces overlaid. If the surfaces do not line up with the BOLD data, you will have to force the registration to a BOLD file that does line up with the surfaces:
#inspect the BOLD/SURF registration for one of the preprocessed files (alternately, could check it for the f.nii file)
SUBJID=FS_0217
freeview -v fmcpr.nii.gz -f ${SUBJECTS_DIR}/${SUBJID}/surf/lh.white:edgecolor=green ${SUBJECTS_DIR}/${SUBJID}/surf/rh.white:edgecolor=green -viewport coronal
Assume that your 3D surface was generated from the T1 MPRAGE from one of those sessions. Use the template BOLD file that was acquired closest in time to this MPRAGE image. For example, if a participant was in SESS1 and SESS2, and the SESS2 T1 file was used to generate the surfaces, coregister the SESS1 BOLD files to the first run of SESS2 (in the example below, this is run 007):
cd 001
#generate a template BOLD image from f.nii
mri_convert --mid-frame f.nii.gz template.nii.gz
#What is the reference (target) image to align to?
TARGET=${SUBJECTS_DIR}/${SUBJID}/bold/007/template.nii.gz
mri_robust_register --mov template.nii.gz --dst ${TARGET} --lta coregister.lta --satit
Then apply the lta coregistration file to the BOLD data:
mri_convert -at coregister.lta f.nii.gz f.nii.gz
Repeat for the remaining SESS1 runs