Nuisance Regression with AFNI 3dDeconvolve

Author

Michael Pascale

Published

July 21, 2025

#| echo: false
#| eval: false
# module load fsl/6.0.7.8

# Check that the NIFTI file is correct.
src/bridle.sh data/tmp/s4101_bold_rest.nii.gz NIFTI-1+ FLOAT32 400

Motion correction was done using mcflirt. The -plots flag provides a parameters timeseries output

(X, Y, Z, r o t X, r o t Y, r o t Z)

with one row per volume.

echo "Line count: $(wc -l < data/tmp/mc/rest_run-1_motion_correction.par)"
head data/tmp/mc/rest_run-1_motion_correction.par | column -t

The -rmsabs flag is undocumented but may be found in the source code. The output associated with this flag is root mean square deviation relative to the reference file with one row per volume.

Note the first line is not zero. The reference image is specified on the mcflirt command line (see CAVEAT/scripts/prep_fx/mcflirt.sh).

echo "Line count: $(wc -l < data/tmp/mc/rest_run-1_motion_correction_abs.rms)"
head data/tmp/mc/rest_run-1_motion_correction_abs.rms

The -rmsrel flag similarly output the relative motion between volumes. As this is relative to the first volume in the series, the first row ( $0.000 m m$ ) is ommitted.

wc -l data/tmp/mc/rest_run-1_motion_correction_rel.rms

Following Toro-Serey, Tobyne, and McGuire (2020), we will censor those volumes where framewise-displacement is $>= 0.5 m m$ . AFNI’s 3dDeconvolve will take a file listing of those volume numbers to the -censor option.

import numpy as np

fd = np.loadtxt('data/tmp/mc/rest_run-1_motion_correction_rel.rms') 
assert fd.shape == (399,), "Expected a file 399x1."

# Censor high motion volumes.
censor = (fd > 0.5).astype(int)

# Censor volumes immediately after high motion.
censor[1:] = censor[1:] | censor[:-1]

# The first volume has no relative motion.
censor = 1 - np.insert(censor, 0, 0)

np.savetxt('data/tmp/s4101_bold_rest_desc-motion_censor.1D', censor, fmt='%d')

array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1])

# module load afni

# Demonstrate the design matrix.
3dDeconvolve \
    -nodata 400 1.475 \
    -censor data/tmp/s4101_bold_rest_desc-motion_censor.1D \
    -polort 5 \
    -num_stimts 1 \
    -stim_file 1 data/tmp/s4101_bold_ventricles.1D \
    -stim_label 1 ventricular \
    -stim_base 1 \
    -ortvec data/tmp/mc/rest_run-1_motion_correction.par motion

x <- read.delim('3ddec_out.1D')
x['t'] <- as.numeric(rownames(x)) * 1.475
plot(x)

mkdir -p data/tmp/afni

# Actually fit the regression.
3dDeconvolve \
    -input data/tmp/s4101_bold_rest.nii.gz \
    -censor data/tmp/s4101_bold_rest_desc-motion_censor.1D \
    -polort A \
    -num_stimts 1 \
    -stim_file 1 data/tmp/s4101_bold_ventricles.1D \
    -stim_label 1 ventricular \
    -stim_base 1 \
    -ortvec data/tmp/mc/rest_run-1_motion_correction.par motion \
    -fitts data/tmp/afni/rest_fit_ts.nii.gz \
    -errts data/tmp/afni/rest_fit_resid.nii.gz \
    -rout -tout -fout \
    -bucket data/tmp/afni/decon \
    -noFDR \
    -xjpeg data/tmp/afni/des_matrix.png \
    -nox1D \

# TODO: must mask brain, voxels of interest?
# TODO: forgot demean baseline option

import nibabel as nib
from nilearn.plotting import view_img
from nilearn import image

res = nib.load('data/tmp/afni/rest_fit_resid.nii.gz')

# Sample an image from teh timeseries.
view_img(image.index_img(res, 10), None, black_bg=True)

References

Toro-Serey, Claudio, Sean M. Tobyne, and Joseph T. McGuire. 2020. “Spectral Partitioning Identifies Individual Heterogeneity in the Functional Network Topography of Ventral and Anterior Medial Prefrontal Cortex.” NeuroImage 205 (January): 116305. https://doi.org/10.1016/j.neuroimage.2019.116305.