Functional Preprocessing Nuisance Regression- derivatives not regressing out

mmahale · October 16, 2022, 2:46am

Hi all,

We are trying to run functional preprocessing nuisance regression with the 3Ddeconvolve command and it doesnt look like the CSF and ventricle derivatives are being regressed out

Here is the script we’re using- any idea why thats happening?

All and any advice is appreciated- thank you in advance!!

echo

echo “Creating WM and ventricular masks…”

cd tmp

# Convert FAST output to MNI space

/usr/local/fsl/bin/applywarp --interp=nn --in=T1_fast_pve_CSF --ref=${subj}_anat_biascorr_MNI -w ${subj}_T1_to_MNI_nonlin_field -o T1_fast_pve_CSF_MNI

/usr/local/fsl/bin/applywarp --interp=nn --in=T1_fast_pve_WM --ref=${subj}_anat_biascorr_MNI -w ${subj}_T1_to_MNI_nonlin_field -o T1_fast_pve_WM_MNI

# Make two masks of the MNI brain: one w/ cerebrum and one w/ ventricles

3dcalc -a /usr/local/fsl/data/standard/MNI152_t1_2mm_strucseg.nii.gz -expr “amongst(a,1,5)” -prefix ${subj}_MNI152_t1_2mm_MNI_cerebrum.nii.gz

3dcalc -a /usr/local/fsl/data/standard/MNI152_t1_2mm_strucseg.nii.gz -expr “equals(a,5)” -prefix ${subj}_MNI152_t1_2mm_MNI_ventricles.nii.gz

# Make and trim the WM and ventricle masks

3dcalc -a ${subj}_MNI152_t1_2mm_MNI_ventricles.nii.gz -b T1_fast_pve_CSF_MNI.nii.gz -expr “abequals(b,1)” -prefix ${subj}_ventricles_all.nii.gz

3dcalc -a ${subj}_MNI152_t1_2mm_MNI_cerebrum.nii.gz -b T1_fast_pve_WM_MNI.nii.gz -expr “abequals(b,1)” -prefix ${subj}_WM_all.nii.gz

##STEP1B##

cd tmp

3dcalc -a ${subj}_WM_all.nii.gz -b a+i -c a-i -d a+j -e a-j -f a+k -g a-k \

-expr ‘a*(1-amongst(0,b,c,d,e,f,g))’ -prefix ${subj}_WM_trim1.nii.gz

3dcalc -a ${subj}_WM_trim1.nii.gz -b a+i -c a-i -d a+j -e a-j -f a+k -g a-k \

-expr ‘a*(1-amongst(0,b,c,d,e,f,g))’ -prefix ${subj}_wm_mask.nii.gz

3dcalc -a ${subj}_ventricles_all.nii.gz -b a+i -c a-i -d a+j -e a-j -f a+k -g a-k \

-expr ‘a*(1-amongst(0,b,c,d,e,f,g))’ -prefix ${subj}_vent_trim1.nii.gz

3dcalc -a ${subj}_vent_trim1.nii.gz -b a+i -c a-i -d a+j -e a-j -f a+k -g a-k \

-expr ‘a*(1-amongst(0,b,c,d,e,f,g))’ -prefix ${subj}_vent_mask.nii.gz

##STEP3##

echo

echo “Extracting regressor signals…”

rm ${subj}_*_regressor.1D

3dmaskave -q -mask ${subj}_wm_mask.nii.gz rest/${subj}_rest_mcf_MNI_b6.nii.gz > ${subj}_wm_regressor.1D

3dmaskave -q -mask ${subj}_vent_mask.nii.gz rest/${subj}_rest_mcf_MNI_b6.nii.gz > ${subj}_vent_regressor.1D

3dmaskave -q -mask rest/${subj}_automask_MNI.nii.gz rest/${subj}_rest_mcf_MNI_b6.nii.gz > ${subj}_gs_regressor.1D

##STEP4##

echo

echo “Computing derivatives…”

cd ${subjectsDirectory}/${subj}/rest/qualitycontrol/motion

1d_tool.py -infile ${subj}_fMRI_vr_motion.1D -derivative -overwrite -write ${subj}_fMRI_vr_motion_deriv.1D

cd …/…/…

1d_tool.py -infile ${subj}_wm_regressor.1D -derivative -overwrite -write ${subj}_wm_regressor_deriv.1D

1d_tool.py -infile ${subj}_vent_regressor.1D -derivative -overwrite -write ${subj}_vent_regressor_deriv.1D

1d_tool.py -infile ${subj}_gs_regressor.1D -derivative -overwrite -write ${subj}_gs_regressor_deriv.1D

echo

cd ${subjectsDirectory}/${subj}

echo “Removing nuisance regressors…”

3dDeconvolve -float -input rest/${subj}_rest_mcf_MNI_b6.nii.gz \

-nfirst 0 -polort 4 \

-mask rest/${subj}_automask_MNI.nii.gz \

-GOFORIT 5 \

-num_stimts 16 \

-stim_file 1 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[0]’ -stim_label 1 roll -stim_base 1 \

-stim_file 2 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[1]’ -stim_label 2 pitch -stim_base 2 \

-stim_file 3 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[2]’ -stim_label 3 yaw -stim_base 3 \

-stim_file 4 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[3]’ -stim_label 4 x -stim_base 4 \

-stim_file 5 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[4]’ -stim_label 5 y -stim_base 5 \

-stim_file 6 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[5]’ -stim_label 6 z -stim_base 6 \

-stim_file 7 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[0]’ -stim_label 7 droll -stim_base 7 \

-stim_file 8 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[1]’ -stim_label 8 dpitch -stim_base 8 \

-stim_file 9 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[2]’ -stim_label 9 dyaw -stim_base 9 \

-stim_file 10 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[3]’ -stim_label 10 dx -stim_base 10 \

-stim_file 11 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[4]’ -stim_label 11 dy -stim_base 11 \

-stim_file 12 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[5]’ -stim_label 12 dz -stim_base 12 \

-stim_file 13 ${subj}_wm_regressor.1D’[0]’ -stim_label 13 wm -stim_base 13 \

-stim_file 14 ${subj}_vent_regressor.1D’[0]’ -stim_label 14 vent -stim_base 14 \

-stim_file 15 ${subj}_wm_regressor_deriv.1D’[0]’ -stim_label 15 dwm -stim_base 15 \

-stim_file 16 ${subj}_vent_regressor_deriv.1D’[0]’ -stim_label 16 dvent -stim_base 16 \

-jobs 2 \

-errts rest/${subj}_rest_mcf_MNI_b6_nuisreg.nii.gz

echo

cd ${subjectsDirectory}/${subj}

rm Decon*

3dDeconvolve -float -input rest/${subj}_rest_mcf_orig_b6.nii.gz \

-nfirst 0 -polort 4 \

-mask rest/${subj}_automask_orig.nii.gz \

-GOFORIT 5 \

-num_stimts 16 \

-stim_file 1 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[0]’ -stim_label 1 roll -stim_base 1 \

-stim_file 2 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[1]’ -stim_label 2 pitch -stim_base 2 \

-stim_file 3 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[2]’ -stim_label 3 yaw -stim_base 3 \

-stim_file 4 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[3]’ -stim_label 4 x -stim_base 4 \

-stim_file 5 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[4]’ -stim_label 5 y -stim_base 5 \

-stim_file 6 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[5]’ -stim_label 6 z -stim_base 6 \

-stim_file 7 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[0]’ -stim_label 7 droll -stim_base 7 \

-stim_file 8 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[1]’ -stim_label 8 dpitch -stim_base 8 \

-stim_file 9 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[2]’ -stim_label 9 dyaw -stim_base 9 \

-stim_file 10 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[3]’ -stim_label 10 dx -stim_base 10 \

-stim_file 11 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[4]’ -stim_label 11 dy -stim_base 11 \

-stim_file 12 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[5]’ -stim_label 12 dz -stim_base 12 \

-stim_file 13 ${subj}_wm_regressor.1D’[0]’ -stim_label 13 wm -stim_base 13 \

-stim_file 14 ${subj}_vent_regressor.1D’[0]’ -stim_label 14 vent -stim_base 14 \

-stim_file 15 ${subj}_wm_regressor_deriv.1D’[0]’ -stim_label 15 dwm -stim_base 15 \

-stim_file 16 ${subj}_vent_regressor_deriv.1D’[0]’ -stim_label 16 dvent -stim_base 16 \

-jobs 2 \

-errts rest/${subj}_rest_mcf_orig_b6_nuisreg.nii.gz

rm Decon*

cd ${subjectsDirectory}/${subj}/rest
echo
echo "Filtering..."
3dTproject -input ${subj}_rest_mcf_MNI_b6_nuisreg.nii.gz -passband 0.008 0.08 -quiet -prefix ${subj}_rest_mcf_MNI_b6_nuisreg_bp.nii.gz 
3dTproject -input ${subj}_rest_mcf_orig_b6_nuisreg.nii.gz -passband 0.008 0.08 -quiet -prefix ${subj}_rest_mcf_orig_b6_nuisreg_bp.nii.gz

end

rickr · October 17, 2022, 4:45pm

Hello,

Let me try to clarify some things.

What is your TR, and how many time points are there?
Are you trying to erode those masks twice? If so, maybe 3dmask_tool would be easier and more direct to use.
Do the regressor and regressor_deriv files look like valid time series?
Why do you think that the time series were not projected out?
Are you attempting to bandpass separate from nuisance regression? This might be part of the problem, in that the output from this will no longer necessarily be orthogonal to the nuisance regressors.

The correct way to do this would be to run 3dDeconvolve using -x1D_stop (i.e. do not actually make an output dataset), just generating the projection matrix X.xmat.1D. Then pass that along with your bandpass options to 3dTproject, so that it can all be part of the same model.

Doing projections in multiple steps has 2 problems:
a. Unless the projection is performed on the second (or later) set of regressors, the second projection will project back in some of what was removed in the first one.
b. Degrees of freedom are not properly accounted for. Bandpassing removes a large fraction of your degrees of freedom, essentially using MANY regressors of no interest. Sometimes just projections would essentially take one to few or even negative remaining DoF.

rick

RJ1 · October 18, 2022, 11:47pm

Hi Rick,

Thanks so much for your reply and suggestions!

TR is .72 with 667 volumes
Initially we were completing it twice, but are no longer trimming the masks
It does look like valid time series
We think the regressors did not project out based on the results of the one-sample and paired-samples t-test, we are still finding signal in the ventricles, WM, cerebellum etc.,
Yes, we first ran nuisance regression and then band pass filtered that processed functional file

In that case if we include -x1D_stop, is the script below correct?
##STEP4##

echo

cd ${subjectsDirectory}/

echo “Removing nuisance regressors…”

3dDeconvolve -float -input rest/rest_mcf_MNI_b6.nii.gz \

-nfirst 0 -polort 4 \

-mask rest/automask_MNI.nii.gz \

-GOFORIT 5 \

-num_stimts 16 \

-stim_file 1 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[0]’ -stim_label 1 roll -stim_base 1 \

-stim_file 2 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[1]’ -stim_label 2 pitch -stim_base 2 \

-stim_file 3 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[2]’ -stim_label 3 yaw -stim_base 3 \

-stim_file 4 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[3]’ -stim_label 4 x -stim_base 4 \

-stim_file 5 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[4]’ -stim_label 5 y -stim_base 5 \

-stim_file 6 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion.1D’[5]’ -stim_label 6 z -stim_base 6 \

-stim_file 7 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[0]’ -stim_label 7 droll -stim_base 7 \

-stim_file 8 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[1]’ -stim_label 8 dpitch -stim_base 8 \

-stim_file 9 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[2]’ -stim_label 9 dyaw -stim_base 9 \

-stim_file 10 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[3]’ -stim_label 10 dx -stim_base 10 \

-stim_file 11 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[4]’ -stim_label 11 dy -stim_base 11 \

-stim_file 12 rest/qualitycontrol/motion/${subj}_fMRI_vr_motion_deriv.1D’[5]’ -stim_label 12 dz -stim_base 12 \

-stim_file 13 ${subj}_wm_regressor.1D’[0]’ -stim_label 13 wm -stim_base 13 \

-stim_file 14 ${subj}_vent_regressor.1D’[0]’ -stim_label 14 vent -stim_base 14 \

-stim_file 15 ${subj}_wm_regressor_deriv.1D’[0]’ -stim_label 15 dwm -stim_base 15 \

-stim_file 16 ${subj}_vent_regressor_deriv.1D’[0]’ -stim_label 16 dvent -stim_base 16 \

-jobs 2 \

-errts X.xmat.1D

-x1D_stop

##STEP6##

rm Decon*

cd ${subjectsDirectory}
echo

echo “Filtering…”
3dTproject -input Decon.xmat.1D -passband 0.008 0.08 -quiet -prefix rest/${subj}_rest_mcf_MNI_b6_nuisreg_bp.nii.gz

Able to complete 3dDeconvolve and generate a .1D file, but 3dTproject doesn’t run.

If you can please assist, thank you so much!

rickr · October 19, 2022, 4:13am

Regarding still finding correlations in the ventricles, that is not a surprise. If you computed principle components within a ventricle mask, there would probably be 3-5 major signals, possibly including some BOLD. So regressing out a mean ventricle signal will not prevent the possibility of correlations within the ventricles (unless all of the ventricles had just a single signal).

Regarding bandpassing in your case, note that you will likely wipe out most of the data with that. To be specific, with a short TR=0.72, the Nyquist frequency is 1/1.44 = 0.694444 (the highest frequency detectable in your data). Bandpassing below 0.08 means keeping approximately 0.08/0.69444 = 11.5% of the frequencies. With the slow components and motion, that means your projection matrix will have about 600 regressors of no interest to project from a 667 time point time series. It might be good to ponder that a bit.

Yes, the updated script looks better. Though heavy bandpassing with such a short TR might have its own consequences.

rick

RJ1 · October 20, 2022, 4:21am

Noted, thank you! We’ll look into the bandpassing with our acquired short TRs. Appreciate the suggestions and assistance!