spm12r
functionsRequires MATLAB (installs SPM to R library)
spm12_slice_timing
- slices are not taken instantaneouslyspm12_realign
- people movespm12_coregister
- structural imaging is higher resolutionspm12_segment
- where’s the gray matter?spm12_normalize
- brains look better in MNIspm12_smooth
- turn that noise downWe will use the subject 113 from the Kirby21 data set to illustrate some basic operations of functional magnetic resonance imaging (fMRI). We will load in the T1 anatomical image and the fMRI from the respective packages.
library(kirby21.t1)
library(kirby21.fmri)
stopifnot(download_fmri_data())
stopifnot(download_t1_data())
functional = get_fmri_filenames(ids = 113, visit = 1)
anatomical = get_t1_filenames(ids = 113, visit = 1)
files = c(anatomical = anatomical,
functional = functional)
files
anatomical
"/Library/Frameworks/R.framework/Versions/3.3/Resources/library/kirby21.t1/visit_1/113/113-01-T1.nii.gz"
functional
"/Library/Frameworks/R.framework/Versions/3.3/Resources/library/kirby21.fmri/visit_1/113/113-01-fMRI.nii.gz"
We know the repetition time (TR) is 2 seconds for this data. It may be encoded in the NIfTI file itself or come from a parameter file from the scanner. We will drop the first 20 seconds to allow for signal stabilization.
library(neurobase)
tr = 2 # seconds
DROP = 10 # 20 seconds for stabilization
fmri = readnii(files["functional"])
times = (DROP + 1):ntim(fmri)
run_fmri = copyNIfTIHeader(fmri, fmri[,,,times], drop = TRUE)
Now run_fmri
contains a nifti
object with the first 10 volumes dropped. We will pass this into the subsequent functions.
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 731659 39.1 1168576 62.5 1168576 62.5
Vcells 48832126 372.6 494071668 3769.5 590072426 4501.9
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 731661 39.1 1168576 62.5 1168576 62.5
Vcells 48832182 372.6 395257334 3015.6 590072426 4501.9
As SPM requires MATLAB and calls all the functions through the matlabr
package, we will have checks in this vignette/workflow that the user has MATLAB. The have_matlab()
function returns a logical that is TRUE
when matlabr
can find MATLAB to run the subsequent commands.
library(matlabr)
have_matlab()
[1] TRUE
If this is not TRUE
, almost none of the functionality below will run because it would simply result in errors.
We will show how to do spatial realignment, slice-timing correction, spatial normalization to the MNI template (2 different ways), and spatial smoothing. Overall, there are many different ways to order these operations, with different options, so this represents just one way to organize a preprocessing pipeline.
Realignment is referring to in this case as within-subject registration of the 4D fMRI data.
library(spm12r)
####################################
# Realignment
####################################
if (have_matlab()) {
realigned = spm12_realign(
filename = run_fmri,
register_to = "mean",
reslice = "mean",
clean = FALSE
)
print(realigned)
}
# Getting Number of Time Points
# Reslice is mean
# Adding SPMDIR: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/spm12r/spm12
# Running script /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpko4uTU/Executable.m
which calls /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpko4uTU/Realign_job.m
# Result is 0
$outfiles
[1] "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpko4uTU/file24dc6012c109.nii"
$rp
[1] "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpko4uTU/rp_file24dc6012c109.txt"
$mean
[1] "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpko4uTU/meanfile24dc6012c109.nii"
$mat
[1] "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpko4uTU/file24dc6012c109.mat"
Overall the spm12_realign
does the realignment. There is some discussion of performing realignment before slice-timing correction because estimation of motion parameters may be skewed after slice-timing correction. We see that the output realigned
has the output 4D fMRI data (outfiles
), the realignment parameters (rp
), voxel-wise mean after realignment (mean
), and the matrix of transformations for the 4D series (mat
).
Here we can read in the rp
file to show the estiamted parameters. These can be used as regressors in motion correction for further analyses.
####################################
# Read in Motion data
####################################
if (have_matlab()) {
rpfile = realigned[['rp']]
rp = read.table(file = rpfile, header = FALSE)
colnames(rp) = c("x", "y", "z",
"roll", "pitch", "yaw")
rp = as.matrix(rp)
print(head(rp))
print(dim(rp))
}
x y z roll pitch
[1,] 0.00000000 0.000000000 -4.410886e-16 -2.507721e-37 2.507721e-37
[2,] -0.03250006 0.022276041 6.105989e-02 -3.616417e-04 -1.389442e-04
[3,] 0.02884634 -0.005311455 5.649025e-02 -3.774759e-04 -1.798493e-04
[4,] 0.02379671 -0.056644719 2.574345e-02 5.180747e-05 -1.756995e-04
[5,] 0.04293237 -0.056942376 2.414558e-02 -2.210120e-05 -8.367122e-05
[6,] 0.02858868 -0.038461745 8.289177e-02 -2.785027e-04 -3.532936e-04
yaw
[1,] 6.288666e-74
[2,] -2.732142e-04
[3,] 1.048090e-04
[4,] -6.489734e-06
[5,] 2.291027e-04
[6,] -8.279382e-06
[1] 200 6
A slice-timing correction does interpolation since each slice was not actually taken at the same time point, but a shifted time point over the course of an entire TR. The correction requires you to input the reference slice (in this case the median, ref_slice
), the repetition time (tr
), time between the first and the last slice within one scan (ta
), and the order the slices were acquired. In our case, it was done in an ascending, contiguous order, so we created the slice order as such. If you used descending or interleaved acquisition, then this must be changed accordingly.
####################################
# Slice Timing Correction
####################################
nslices = oro.nifti::nsli(run_fmri)
slice_order = 1:nslices
ref_slice = slice_order[median(seq(nslices))]
ta = tr - tr/nslices
n_time_points = ntim(run_fmri)
if (have_matlab()) {
aimg = spm12_slice_timing(
filename = realigned[['outfiles']],
nslices = nslices,
tr = tr,
slice_order = slice_order,
ta = ta,
ref_slice = ref_slice,
prefix = "a",
clean = FALSE,
retimg = FALSE)
print(aimg)
mean_img = realigned[["mean"]]
mean_nifti = readnii(mean_img)
}
# Getting Number of Time Points
# Adding SPMDIR: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/spm12r/spm12
# Running script /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpko4uTU/Executable.m
which calls /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpko4uTU/Slice_Timing_job.m
# Result is 0
[1] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpko4uTU/file24dc2f9e99f4/afile24dc6012c109.nii"
We see the output aimg
has the filename of the slice-timing corrected 4D image. We also read in the mean image into a nifti
object (mean_nifti
).
For the subsequent image normalization steps, SPM assumes the data is aligned along the anterior commissure (AC) posterior commissure (PC) line (AC-PC). The acpc_reorient
function (based on nii_setOrigin
from Dr. Chris Rorden) will do this. The syntax is that the first file (mean_img
) is used to estimate the line/plane and the subsequent files are reoriented using this estimation (aimg
). These are changes to the header of the image and the image with the new header is written to the same file as the input file.
if (have_matlab()) {
acpc_reorient(
infiles = c(mean_img, aimg),
modality = "T1")
}
#Reorientation /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpko4uTU/meanfile24dc6012c109.nii
/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpko4uTU/file24dc3b6bbd09.m
/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpko4uTU/meanfile24dc6012c109.nii
"/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpko4uTU/meanfile24dc6012c109.nii"
/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpko4uTU/file24dc2f9e99f4/afile24dc6012c109.nii
"/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpko4uTU/file24dc2f9e99f4/afile24dc6012c109.nii"
Here, we will spatially normalize the fMRI to the MNI template using spm12_normalize
. Here we specify a standard bounding box for a 2mm MNI template. We are taking the mean image and directly registering it to the MNI template (T1-weighted), and applying that transform to the other.files
, in this case the mean image and the 4D fMRI image.
if (have_matlab()) {
bbox = matrix(
c(-90, -126, -72,
90, 90, 108),
nrow = 2, byrow = TRUE)
print(bbox)
direct_norm = spm12_normalize(
filename = mean_img,
other.files = c(mean_img, aimg),
bounding_box = bbox,
clean = FALSE
)
print(direct_norm)
dnorm_files = direct_norm$outfiles
dnorm_mean_img = readnii(dnorm_files[1])
}
[,1] [,2] [,3]
[1,] -90 -126 -72
[2,] 90 90 108
# Adding SPMDIR: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/spm12r/spm12
# Running script /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpko4uTU/Executable.m
which calls /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpko4uTU/Normalize_Estimate_and_Write_job.m
# Result is 0
$outfiles
[1] "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpko4uTU/wmeanfile24dc6012c109.nii"
[2] "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpko4uTU/file24dc2f9e99f4/wafile24dc6012c109.nii"
We see tbe output direct_norm
is a list with the output files (outfiles
). The order of these files is the same as the order of those put in. In this case, the first file of outfiles
is the normalized mean image and second is the normalized 4D image. Here we read in the spatially normalized mean image to compare to the template later.
Indirect normalization refers to spatially normalizing the co-registered anatomical T1-weighted image to the MNI template. This transformation is applied to the mean image and 4D fMRI image. This is also referred to as Unified Segmentation.
Here we will perform the registration of the T1-weigted anatomical image into the space of the mean fMRI image after realignment. This is referred to as “co-registration” as it is within-subject registration, but across modalities (where we referred to within-subject, within-modality as realignment).
Here, we also reorient the anatomical image the AC-PC line. We then perform the coregistration using spm12_coregister
, where the fixed image is the mean image and the moving image is the anatomical.
if (have_matlab()) {
anatomical = files["anatomical"]
anat_img = checknii(anatomical)
print(anat_img)
acpc_reorient(
infiles = anat_img,
modality = "T1")
coreg = spm12_coregister(
fixed = mean_img,
moving = anat_img,
prefix = "r")
coreg_anat = coreg$outfile
coreg_img = readnii(coreg_anat)
double_ortho(coreg_img, mean_nifti)
}
[1] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpko4uTU/113-01-T1.nii"
attr(,"nbrOfBytes")
[1] 44564832
#Reorientation /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpko4uTU/113-01-T1.nii
/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpko4uTU/file24dc225e18ff.m
# Adding SPMDIR: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/spm12r/spm12
# Running script /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpko4uTU/Executable.m
which calls /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpko4uTU/Coregister_job.m
# Result is 0
# Removing scripts