Extracting ROIs for PPI Analysis Using SPM Batch
Posted on Wed 03 January 2018 in Neuroscience
A psycho-physiological interactions (PPI) analysis of fMRI data requires that some kind of region of interest be extracted from a subject-level contrast of interest. For small studies of only a few participants, manually extracting these ROIs using SPM's GUI is simple enough, but for larger studies, SPM's batch processing capabilities really pay off. Unfortunately, SPM's terse documentation of its batch ROI extraction utility often leaves users guessing how to use it effectively. In this post I am going to discuss a few of the important things to know about ROI extraction in SPM, and provide a few examples of how to automate ROI extraction.
Good Things to Know First
Extracting an ROI requires that you make a number of important decisions. You will need to decide what your ROI will look like, where it will be, and whether and by how much to threshold the image before extraction.
ROI Types
The location, size, and shape of your ROI will depend on your use-case. One particularly common and well-motivated choice is to use masks of pre-defined anatomical regions, in which case the ROI will typically be stationary across subjects, though the subset of voxels that are ultiimately included depends on whether and how they are thresholded. Another commonly used type of ROI is a spherical ROI centered around a certain voxel, such as a peak in a group-level contrast; for ROIs of this type, you may want to allow the exact location of the sphere to vary between subjects.
Functional Heterogeneity
ROIs may not be functionally homogenous, in which case a single summary measure such as the arithmetic mean may poorly represent activity in the ROI [1]. Furthermore, not all activity in the ROI may be of interest. Thresholding the SPM with respect to some contrast reduces functional heterogeneity by excluding the contributions of insignificant voxels to the extracted timecourse, relative to the effect of interest. Another way of addressing within-ROI functional heterogeneity is by taking the first eigenimage, a weighted mean of the ROI. SPM's VOI utility, spm_regions.m, permits a hybrid approach: it extracts the eigenimage, but also asks the user to specify a threshold as well, which can be uncorrected ('none'), or corrected for multiple comparisons ('FWE'). Setting the threshold to 1 means no voxels are excluded.
Adjustment
At the subject level, within ROI voxel data is extracted directly from your preprocessed data; it is then whitened and filtered using the high-pass filter specified in your SPM. From spm_regions.m:
%-Get raw data, whiten and filter
%--------------------------------------------------------------------------
y = spm_get_data(SPM.xY.VY,xSPM.XYZ(:,Q));
y = spm_filter(SPM.xX.K,SPM.xX.W*y);
Notably, this step does not exclude any effects of no-interest that were regressed out in your first-level analysis. That means, for example, that if you included realignment parameters as nuisance regressors in your first-level analysis, any variance explained by them is still in your ROI data. However, spm_regions.m gives the user the opportunity to adjust the data to remove unwanted variance.
Specifically, the user can decline to adjust the data by entering in 0, can adjust against a specific contrast, usually an F-contrast, by entering the index of that contrast, or they can adjust for everything by entering in NaN. The difference between these options can be obscure without looking at the code:
%-Remove null space of contrast
%--------------------------------------------------------------------------
if xY.Ic ~= 0
%-Parameter estimates: beta = xX.pKX*xX.K*y
%----------------------------------------------------------------------
beta = spm_get_data(SPM.Vbeta,xSPM.XYZ(:,Q));
%-subtract Y0 = XO*beta, Y = Yc + Y0 + e
%----------------------------------------------------------------------
if ~isnan(xY.Ic)
y = y - spm_FcUtil('Y0',SPM.xCon(xY.Ic),SPM.xX.xKXs,beta);
else
y = y - SPM.xX.xKXs.X * beta;
end
end
xY.Ic is the adjustment parameter. We see from the code that if it is 0, then the block is skipped entirely: the data is only whitened and filtered. If it is NaN, then the entire whitened and filtered design matrix, specified in SPM.xX.xKXs.X, is multiplied with your SPM's beta estimates and subtracted rom the whitened and filtered ROI Y values. Thus, when adjusting for everything, you are getting back the variance not accounted for by your model, i.e. the residuals. Otherwise, xY.Ic is taken as an index of a contrast. A call to spm_FcUtil.m returns the Y0 values explained by the design-matrix corresponding to the columns that are all zeros in the supplied contrast [2].
So, if, for example, you supplied an index to an F-contrast like:
1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
where the first three columns are effects of interest, and the last six are realignment parameters, then variance explained by those six realignment parameters will be subtracted out from your Ys.
Recipes to Extract ROIs
I've sketched out some code that can be used to extract ROIs in several typical scenarios. The heart of the code is a matlabbatch cell arry of structures, jobs, passed to SPM's (which in turn calls spm_regions.m for each job). However, to illustrate its use I will fill out the code somewhat, and encapsulate the routine within a function, which I will assume to be on the Matlab path.
Fixed Anatomical ROI
When extracting with a fixed anatomical ROI, all we really need to do in perform a conjunction on the brain map and the anatomical mask. The conjunction is specified by the boolean expression i1 & i2, where i1 and i2 refer to roi{1} and roi{2} respectively.
function extract_anatomical_rois(cfg)
% Extracts ROIs from one subject.
% Arguments:
% cfg: A struct with the following fields
% target_dir: path to subject's first level analysis
% results.
%
% adjust: 0 to not adjust, NaN to adjust
% for everything, else index of F-contrast to use.
%
% contrast: Index of contrast of interest.
%
% threshold: The p value at which to threshold the
% contrast image.
%
% threshdesc: Whether to correct for comparisons.
% 'none' for no correction, 'FWE' for correction.
%
% extent: Exclude clusters with number of significant
% voxels that is less than extent.
%
% name: Name of VOI output.
% 'VOI_' is prefixed to it by spm_regions.m
%
% mask: String paths to ROI mask.
%
spm('defaults', 'FMRI');
spm_jobman('initcfg');
% spm_regions needs to be in SPM.mat directory.
starting_dir = pwd;
cd(cfg(1).target_dir);
for ii = 1 : numel(cfg)
jobs{ii}.spm.util.voi.spmmat = {fullfile(...
cfg(1).target_dir, 'SPM.mat')};
jobs{ii}.spm.util.voi.adjust = cfg(ii).adjust;
jobs{ii}.spm.util.voi.session = cfg(ii).session;
jobs{ii}.spm.util.voi.name = cfg(ii).name;
% The image data
jobs{ii}.spm.util.voi.roi{1}.spm.spmmat = {''};
jobs{ii}.spm.util.voi.roi{1}.spm.contrast = cfg(ii).contrast;
jobs{ii}.spm.util.voi.roi{1}.spm.threshdesc = cfg(ii).threshdesc;
jobs{ii}.spm.util.voi.roi{1}.spm.thresh = cfg(ii).threshold;
jobs{ii}.spm.util.voi.roi{1}.spm.extent = cfg(ii).extent;
% The anatomical ROI mask
jobs{ii}.spm.util.voi.roi{2}.mask.image = {cfg(ii).mask)};
% The logical expression combining images.
jobs{ii}.spm.util.voi.expression = 'i1 & i2';
end
try
spm_jobman('run', jobs);
catch ME
cd(starting_dir);
rethrow(ME);
end
cd(starting_dir)
end
For a group of subjects, one might use this like so:
function run_extract_anatomical_rois()
data_root_dir = '/path/to/root/data/directory';
subject_folder_pattern = 'SUB*';
masks = {
'/path/to/dACC.nii'
'/path/to/left_insula.nii'
'/path/to/right_insula.nii'
};
names = {
'dACC'
'lInsula'
'rInsula'
};
threshold = 0.001;
threshdesc = 'none';
session = 1;
extent = 0;
F_contrast_index = 1;
contrast_of_interest_index = 3;
subject_folders = dir(fullfile(data_root_dir, ...
subject_folder_pattern));
for ii = 1 : numel(subject_folders)
% parfor is good alternative if available
cfg = struct();
for jj = 1 : numel(masks)
cfg(jj).target_dir = fullfile(data_root_dir, ...
subject_folders(ii).name); % should stay same
cfg(jj).mask = masks{jj}
cfg(jj).adjust = F_contrast_index;
cfg(jj).session = session;
cfg(jj).name = names{jj}
cfg(jj).contrast = contrast_of_interest_index;
cfg(jj).threshold = threshold;
cfg(jj).threshdesc = threshdesc;
cfg(jj).extent = extent;
end
extract_anatomical_rois(cfg);
end
end
Fixed Geometric ROI
Similarly to a fixed anatomical ROI, one can define an ROI using a geometric shape like a box or a sphere. One common use-case for this is when you wish to extrace ROIs around a peak in a group-level contrast. The code is very similar to the preceding, except that instead of supplying a mask, one supplies the relevant parameters of the geometric mask. For example, to create a spherical ROI such that the center of the sphere is identical for every subject, one could do the following:
function extract_spherical_rois(cfg)
% Extracts ROIs from one subject.
% Arguments:
% cfg: A struct with the following fields
% target_dir: path to subject's first level analysis
% results.
%
% adjust: 0 to not adjust, NaN to adjust
% for everything, else index of F-contrast to use.
%
% contrast: Index of contrast of interest.
%
% threshold: The p value at which to threshold the
% contrast image.
%
% threshdesc: Whether to correct for multiple comparisons.
% 'none' for no correction, 'FWE' for correction.
%
% extent: Exclude clusters with number of significant
% voxels that is less than extent.
%
% name: Name of VOI output.
% 'VOI_' is prefixed to it by spm_regions.m
%
% center: 1x3 matrices of coordinates
% of sphere center, e.g. [3 14 49]
%
% radius: nx1 matrix of sphere radius, in millimeters
spm('defaults', 'FMRI');
spm_jobman('initcfg');
% spm_regions needs to be in SPM.mat directory.
starting_dir = pwd;
cd(cfg(1).target_dir);
for ii = 1 : numel(cfg)
jobs{ii}.spm.util.voi.spmmat = {fullfile(...
cfg(ii).target_dir, 'SPM.mat')};
jobs{ii}.spm.util.voi.adjust = cfg(ii).adjust;
jobs{ii}.spm.util.voi.session = cfg(ii).session;
jobs{ii}.spm.util.voi.name = cfg(ii).name;
% The image data
jobs{ii}.spm.util.voi.roi{1}.spm.spmmat = {''};
jobs{ii}.spm.util.voi.roi{1}.spm.contrast = cfg(ii).contrast;
jobs{ii}.spm.util.voi.roi{1}.spm.threshdesc = cfg(ii).threshdesc;
jobs{ii}.spm.util.voi.roi{1}.spm.thresh = cfg(ii).threshold;
jobs{ii}.spm.util.voi.roi{1}.spm.extent = cfg(ii).extent;
% The ROI mask
jobs{ii}.spm.util.voi.roi{2}.sphere.centre = cfg(ii).center;
jobs{ii}.spm.util.voi.roi{2}.sphere.radius = cfg(ii).radius;
jobs{ii}.spm.util.voi.roi{2}.sphere.move.fixed = 1;
% The conjunctive expression combining images
jobs{ii}.spm.util.voi.expression = 'i1 & i2';
end
try
spm_jobman('run', jobs);
catch ME
cd(starting_dir);
rethrow(ME);
end
cd(starting_dir)
end
Flexible Geometric ROIs
One potential disadvantage in the approach take above is that it does not allow for inter-subject heterogeneity. For example, suppose we had a peak in a group-level T-test at MNI coordinates [3 14 49]. It is very probably the case that the peak activation for each subject is not exactly at [3 14 49], but is, instead, clustered around that. Arguably, it may therefore be a good idea to adjust each subject's ROI to be centered at their own individualized peak.
However, me must be careful to bound how far the ROI can move, otherwise some subjects ROIs may end up in completely different parts of the brain. To do this, we need to introduce a bounding region.
In this example, we will take the conjunction of the first and third image, and constrain the location of the third to within the second, a bounding sphere:
function extract_spherical_rois(cfg)
% Extracts ROIs from one subject.
% Arguments:
% cfg: A struct with the following fields
% target_dir: path to subject's first level analysis
% results.
%
% adjust: 0 to not adjust, NaN to adjust
% for everything, else index of F-contrast to use.
%
% contrast: Index of contrast of interest.
%
% threshold: The p value at which to threshold the
% contrast image.
%
% threshdesc: Whether to correct for multiple comparisons.
% 'none' for no correction, 'FWE' for correction.
%
% extent: Exclude clusters with extent less than extent.
%
% name: Name of VOI output.
% 'VOI_' is prefixed to it by spm_regions.m
%
% bounding_center: 1x3 matrix of coordinates
% of bounding sphere centers, e.g. [3 14 49]
%
% bounding_radius: nx1 matrix of bounding sphere radius, in
% millimeters
%
% roi_radius: nx1 matrix of ROi radius, in millimeters.
spm('defaults', 'FMRI');
spm_jobman('initcfg');
% spm_regions needs to be in SPM.mat directory.
starting_dir = pwd;
cd(cfg(1).target_dir);
for ii = 1 : numel(cfg)
jobs{ii}.spm.util.voi.spmmat = {fullfile(...
cfg(ii).target_dir, 'SPM.mat')};
jobs{ii}.spm.util.voi.adjust = cfg(ii).adjust;
jobs{ii}.spm.util.voi.session = cfg(ii).session;
jobs{ii}.spm.util.voi.name = cfg(ii).name;
% The image data
jobs{ii}.spm.util.voi.roi{1}.spm.spmmat = {''};
jobs{ii}.spm.util.voi.roi{1}.spm.contrast = cfg(ii).contrast;
jobs{ii}.spm.util.voi.roi{1}.spm.threshdesc = cfg(ii).threshdesc;
jobs{ii}.spm.util.voi.roi{1}.spm.thresh = cfg(ii).threshold;
jobs{ii}.spm.util.voi.roi{1}.spm.extent = cfg(ii).extent;
% The mask bounding the ROI's mask.
jobs{ii}.spm.util.voi.roi{2}.sphere.centre = cfg(ii).bounding_center;
jobs{ii}.spm.util.voi.roi{2}.sphere.radius = cfg(ii).bounding_radius;
jobs{ii}.spm.util.voi.roi{2}.sphere.move.fixed = 1;
% The extracted ROI's mask. Note how it references the second image as
% its mask. 'move.local' is used to move to nearest local maximum
% in mask. Use 'move.global' to instead move to largest peak in mask.
jobs{ii}.spm.util.voi.roi{3}.sphere.centre = [0 0 0]; % not used
jobs{ii}.spm.util.voi.roi{3}.sphere.radius = cfg(ii).roi_radius
jobs{ii}.spm.util.voi.roi{3}.sphere.move.local.spm = 1;
jobs{ii}.spm.util.voi.roi{3}.sphere.move.local.mask = 'i2'
% The conjunctive expression combining images
jobs{ii}.spm.util.voi.expression = 'i1 & i3';
end
try
spm_jobman('run', jobs);
catch ME
cd(starting_dir);
rethrow(ME);
end
cd(starting_dir)
end
Since the ROI mask is masked by the bounding mask, if the center of the ROI mask moves too close to the edge of the bounding mask, you may find that your ROIs are smaller than expected. When allowing ROIs to move, it is important to perform quality assurance on your extracted ROIs before using them in an analysis.
An attractive alternative to a bounding sphere or cube is to use an anatomical mask to constrain the possible locations of the extracted ROI. To do that, one simply replaces the second ROI in the code above with the mask, So, instead of:
jobs{ii}.spm.util.voi.roi{2}.sphere.centre = cfg(ii).bounding_center;
jobs{ii}.spm.util.voi.roi{2}.sphere.radius = cfg(ii).bounding_radius(ii);
jobs{ii}.spm.util.voi.roi{2}.sphere.move.fixed = 1;
we use:
matlabbatch{1}.spm.util.voi.roi{2}.mask.image = {cfg(ii).mask};
And that should do it!
| [1] | Tong, Y., Chen, Q., Nichols, T. E., Rasetti, R., Callicott, J. H., Berman, K. F., … Mattay, V. S. (2016). Seeking Optimal Region-Of-Interest (ROI) Single-Value Summary Measures for fMRI Studies in Imaging Genetics. PLoS ONE, 11(3), e0151391. http://doi.org/10.1371/journal.pone.0151391 |
| [2] | As of writing, spm_FcUtil.m indicates, in comments in the code, that it checks whether the supplied contrast is an F-contrast and raises an error otherwise. Indeed, its help documentation indicates that spm_FcUtil.m can be used for just that purpose by calling it with the argument IsFcon. However, inspection of the code shows that it only checks whether the supplied contrast structure has all the necessary fields for an F-contrast, which are also necessary fields for T-contrasts; testing confirms that it classifies a T-contrast as an F-contrast, which, I suppose, it is, sort of. |