fcMRI Processing Tools

The fcMRI processing script is a wrapper around several SPM tools. The processing script is paired with a parameter file, and enables simple batching of data preprocessing.

The steps involved in the preprocessing include:

1. Dropping initial dummy volumes

2. Slice-time correction

3. Realignment

4. Normalization

5. Smoothing

6. Motion regression

7. Temporal filtering

8. Regression of nuisance signal

9. Automated quality checks on the data

In this implementation functional runs can be processed separately or together. For instance we could do

cd Session1
fn = {'Run1.nii' 'Run2.nii'};
fcPreproc(fn,'fcParams')

Alternatively we can analyze each run separately, but this will require that each run is in its own folder

cd Session1/Run1
fn = 'Run1.nii';
fcPreproc(fn,'fcParams')
cd ../Run2
fn = 'Run2.nii';
fcPreproc(fn,'fcParams')

Importantly, you must be in the same folder as the files being analyzed.

All configuration options are set in the Parameter function.

The parameter is specified as a function. This allows us to call P = fcParams, and get the full listing of all the processing parameters

The Parameter function is organized as follows):

function P = fcParams

Here we specify the default order of preprocessing steps.

P.PO = [];
P.PO{1} = 'nfn = drop_vols(fn)';
P.PO{2} = 'nfn = slice_time(nfn);';
P.PO{3} = 'nfn = realign(nfn);';
P.PO{4} = 'nfn = normalize(nfn);';
P.PO{5} = 'nfn = smooth(nfn);';
P.PO{6} = 'nfn = motion_regress(nfn);';
P.PO{7} = 'nfn = filter_data(nfn);';
P.PO{8} = 'nfn = nuisance_regress(nfn);';
P.PO{9} = 'Bad_Vols(nfn)';

First we specify the TR for the functional acquisitions.

P.TR = 3;  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Below are the parameters specific to each of the sub-functions. First up is dropping volumes, which we employ for acquisitions where the dummy volumes were not automatically discarded.

We can specify whether we want this step to be performed, which volumes should be dropped, and the filename prefix for the Nifti file with the dropped volumes.

% Drop Specified Volumes
P.dv.do = 1;
P.dv.dropVols = [1 2 3 4];
P.dv.prefix = 'dv_';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Next comes slice timing, and again we specify whether this step should be performed, provide the parameters for the slice timing procedure, and specify the output file prefix

% Slice Timing Paramters: See SPM8 documentation for more info.
P.st.do = 1;
P.st.sliceorder = [1:2:47 2:2:46]; % This is the order in which the slices were acquired
P.st.refslice = 1;
P.st.timing = [(P.TR/numel(P.st.sliceorder)) (P.TR/numel(P.st.sliceorder))-((P.TR/numel(P.st.sliceorder))/numel(P.st.sliceorder))];
P.st.prefix = 'st_';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Next are the realignment parameters. You can choose between inria_realign and spm_realign. If you pass in multiple runs, the runs will be realigned to one another, and it will produce a single mean realigned image.

Note also that the realigned files are not resliced, rather the changes in orientation are written to the headers of the input files.

%%% Realign: See SPM8 documentation for more info.
% Realign and Reslice Parameters
P.rr.do = 1;  %%% 1 =  do;  0 =  skip;
P.rr.RealignmentType = 2;  % 1 = use inria_realign;         2 = use spm_realign

%%% inria_realign specific parameters
P.rr.RealignParsI.quality = .95;  % Quality
P.rr.RealignParsI.fwhm    = 5;  % pre alignment smoothing kernel in mm.
P.rr.RealignParsI.sep     = 3;  % not 100% sure what this does.
P.rr.RealignParsI.interp  = 5;  % order of interpolation function for resampling.
P.rr.RealignParsI.wrap    = [0 0 0]; %% No wrapping.
P.rr.RealignParsI.rho_func = 'geman';
P.rr.RealignParsI.cutoff = 2.5;
P.rr.RealignParsI.hold = -9;

%%% spm_realign realignment paramters
P.rr.RealignParsS.quality = .50;  % Quality
P.rr.RealignParsS.fwhm    = 5;  % pre alignment smoothing kernel in mm.
P.rr.RealignParsS.sep     = 3;  % not 100% sure what this does.
P.rr.RealignParsS.interp  = 5;  % order of interpolation function for resampling.
P.rr.RealignParsS.wrap    = [0 0 0]; %% No wrapping.
P.rr.RealignParsS.rtm     = 0;  % no Register to mean (second pass alignment).
P.rr.RealignParsS.PW      = ; % Weigthing?

%%% Reslice: See SPM8 documentation for more info.
P.rr.ReslicePars.mask   = 1; %% if one timepoint in a voxel is 0 all  timepoints in that voxel are set to 0;
P.rr.ReslicePars.mean   = 1; %% Write a Mean image.
P.rr.ReslicePars.interp = 5; %% interpolation order for resampling;
P.rr.ReslicePars.which  = 3; %% Reslice all images including the first image
P.rr.ReslicePars.prefix = 'rr_';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Then we provide the normalization parameters. The default implementation is to normalize the mean image to the MNI 152 EPI template.

% Normalization Parameters: See SPM8 documentation for more info.
P.nn.do = 1;
P.nn.source = 'mean*.nii';
tmp = which('spm');
ind = find(tmp==filesep);
whe = tmp(1:ind(end)-1);
P.nn.template = [whe '/templates/epi.nii'];

%%% Normalize data to MNI epi template: See SPM8 documentation for more info.
P.nn.NormPars.smosrc = 8;
P.nn.NormPars.smoref = 0;
P.nn.NormPars.regtype = 'mni';
P.nn.NormPars.cutoff = 25;
P.nn.NormPars.nits = 16;
P.nn.NormPars.reg = 1;
%%% Warping Paramters:  See SPM8 documentation for more info.
P.nn.rflags.preserve = 0;
P.nn.rflags.bb = [-78 -112 -70; 78 76 90]; 
P.nn.rflags.vox = [3 3 3];
P.nn.rflags.interp = 5; 
P.nn.rflags.wrap = [0 0 0];
P.nn.rflags.prefix =  'nn_';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Here we specify the parameters for spatial smoothing.

% Smoothing Parameters % 1 = Perform Smoothing. 0 = Skip it.
P.ss.do = 1;
P.ss.prefix = 'ss_';
P.ss.kernel = [6 6 6]; % 3D smoothing kernel in mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Next are the parameters for the automated QA. Automated QA looks at various data quality metrics to figure out if the data should be used, and if there are specific volumes that should be excluded from the analysis.

The below parameters will junk a run if the Global SNR is less that 115, if the subject changed position by more than 5mm over the course of the run, or rotated by more than % degrees over the course of the run, or if the mean movement across the run is greater than 0.2 mm/TR.

We will also flag volumes as bad if the mean bold signal for the volume is outside 2 standard deviations, if the movement between successive TRs is greater than .75mm, or if there is more than 1.5 degrees of rotation between TRs.

These parameters were arrived at by looking at a study specific sample of elderly subjects and estimating what appeared to be normal and abnormal for the population. You will likely want to adjust these parameters for your own data.

P.bv.do = 1;

P.bv.CheckGlobalSNR = 1;
P.bv.SNR_prefix = 'Global_SNR_Report';
P.bv.SNR_suffix = '.txt';
P.bv.SNRthresh = 115;       %%% The global SNR threshold that must be met.  If a run has less than this value it will be junked.
                            %%% liberal value = 115  moderate value = 175 conservative value =  235
P.bv.RealignPrefix = 'realignment_';
P.bv.GlobalSigThresh = 2.5; %% Units are Standard Deviations to assess outliers of global signal velocity
P.bv.MoveThresh = .75;      %% Units are mm/TR
P.bv.RotThresh = 1.5;       %% Units are degrees/TR
P.bv.BadVolThresh = 20;     %% The maximum number of Bad Volumes allowed before the entire run is junked.
P.bv.MeanMovementThresh = .2; %% average change in position between TRs
P.bv.TotalMovementThresh = 5; %% maximum difference in Position in mm
P.bv.TotalRotationThresh = 5; %% maximum difference in Rotation in degrees
P.bv.LogFileName = 'AnalysisLog.txt';
P.bv.BadVolRegsName = 'ExtraRegressors.mat';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Next we regress out motion. MotionDeriv is a flag the specifies whether or not to also add the first derivative of the motion parameters as well.

% Motion Regression Params
P.mot.do = 1;
P.mot.MotionDeriv = 1;
P.mot.prefix = 'motRes_';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Next we filter the data. There are options to perform linear detrending, high pass filter only, or band-pass filtering.

The parameters below will perform band-pass filtering with a pass-band of 0.01-0.08 Hz. Filtering is performed by a butterworth filter with a specified filter order of 4.

% Filtering Parameters
P.ff.do = 1; % 1 = Perform Filtering. 0 = Skip it.
P.ff.Detrend = 0; % 1 = Perform Detrending. 0 = Skip it.
P.ff.HighCut = 0.08;
P.ff.LowCut = 0.01;
P.ff.FilterOrder = 4;
P.ff.prefix = 'ffBPS_';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Finally we perform nuisance regression. The specification here entails providing a set of mask images representing regions whose signal you want to remove.

In this example we have masks for white matter, ventricles, and whole brain, and we also include the first derivatives of each of those signals.

% Nuisance Regressor Parameters
P.res.do = 1; % 1 = Perform the Nuisance variable regression. 0 = Skip it.
P.res.Masks = {...
               which('avg152T1_brain_MNI.nii') ... % Masks used for computing physiological signals to be regressed out. Add as many as you would like.
               which('avg152T1_ventricles_MNI.nii') ...
               which('avg152T1_WM_MNI.nii') ...
               };
P.res.MaskDeriv = 1; % 1 = Also regress out the first derivative of physio vars. 0 = Do NOT regress out the first derivative of physio vars.
P.res.prefix = 'res_';

Once the parameters are specified, the parameter script can be passed into the processing script.

fcPreproc('filename.nii','fcParams')

or

fcPreproc({'filename1.nii' 'filename2.nii'},'fcParams')

fcPreproc has a third input parameter that can be used to selectively run certain preprocessing steps.

For instance if I want to perform bandpass filtering as the final step I can do:

fcPreproc('filename.nii','fcParams', [1 2 3 4 5 6 8 7]);

Additionally I can choose to only run certain steps. For instance if I only want to perform motion regression I can do:

fcPreproc('filename.nii','fcParams', 6);

Importantly fcPreproc will start with whatever image you pass it, so the last command should probably be something like:

fcPreproc('ss_nn_st_dv_filename.nii','fcParams', 6);

That’s all there is to it. These functions can then be used to batch process multiple subjects and multiple sessions.

As another example check out RestVolPath.m this is a wrapper function that will only perform the processing steps that have not been completed. You pass in the original Nifti filename, and the script will use the specified prefixes in the param file to see which steps have been performed and which have not, and it will then run those steps that have not been performed yet. (Note that if a later processing step exists it assumes that earlier steps should not be re-run).

I use this type of script a lot to loop over large datasets and make sure that everything is processed.

Share this page: