Adding a field map for rs-fMRI distortion correction

Hi there!

I’ve been using the same lovely AFNI pipeline to pre-process the various datasets that I’m working on.

We’ve recently started a new study where we’ll be collecting field maps in the hope of doing distortion correction and improving the quality of our data.
I was wondering how exactly I could integrate field maps into the pipeline without breaking something.

The command(s) I’m using to generate the script is as follows:
module load afni/21.1.07
module load anaconda/3.8
module load tedana/0.0.11

afni_proc.py -blocks despike tshift align tlrc volreg mask combine blur scale regress -copy_anat t1.nii.gz -anat_has_skull yes -dsets_me_echo fmri_e1.nii.gz -dsets_me_echo fmri_e2.nii.gz -dsets_me_echo fmri_e3.nii.gz -dsets_me_echo fmri_e4.nii.gz -echo_times 12.9 28.9 44.9 60.9 -reg_echo 2 -tcat_remove_first_trs 10 -align_opts_aea -cost nmi -check_flip -skullstrip_opts -push_to_edge -tlrc_base MNI152_T1_2009c+tlrc -tlrc_NL_warp -volreg_align_to MIN_OUTLIER -volreg_align_e2a -volreg_tlrc_warp -mask_epi_anat yes -mask_segment_anat yes -mask_segment_erode yes -combine_method tedana -blur_size 6 -blur_in_mask yes -regress_bandpass 0.01 0.1 -regress_ROI WMe CSF -regress_motion_per_run -regress_censor_motion 3.0 -regress_censor_outliers 0.05 -regress_apply_mot_types demean deriv -regress_est_blur_epits -html_review_style pythonic

In the combine block code, I remove the tedana wrapper bit and add this instead:
tedana -d pb03.SUBJ.r01.e01.volreg+tlrc.BRIK
pb03.SUBJ.r01.e02.volreg+tlrc.BRIK
pb03.SUBJ.r01.e03.volreg+tlrc.BRIK
pb03.SUBJ.r01.e04.volreg+tlrc.BRIK
-e 12.9 28.9 44.9 60.9 --mask mask_epi_anat.SUBJ+tlrc.BRIK --debug

Finally, I change the 3dcopy bit to:
3dcopy desc-optcomDenoised_bold.nii.gz pb04.$subj.r$run.combine+orig

I’ve got a field map (in radians), which has been unwrapped.
https://afni.nimh.nih.gov/pub/dist/doc/htmldoc/programs/epi_b0_correct.py_sphx.html seems like the thing I should be using but https://afni.nimh.nih.gov/pub/dist/doc/htmldoc/programs/3dQwarp_sphx.html is the thing that is used in the blip block.
I’m not exactly sure how to go about things, so any tips would be highly appreciated!

Looking forward to your reply!

Best Wishes,
Misho

Hi Misho,

Well, I have good news and bad news.

The good news is that you can use tedana directly in afni_proc.py (without the wrapper) by applying it as the “m_tedana” combine method. Give that a try, and compare it with a previous result.

The bad news is that I have not gotten around to applying epi_b0_correct.py properly in afni_proc.py (it almost happened a few years back). Paul and I were pondering some of the details, but then other things came up.

From what I recall, since there are many options to consider, we thought it might be easiest (at least initially) to have the user run epi_b0_correct.py to generate a distortion correction matrix, and then pass it to afni_proc.py, which could then be applied akin to using the blip results. But still, I have not done it.

We could revisit that if you are in a situation to apply it…

  • rick

Hi Rick,

Thanks for getting back to me so promptly!

I didn’t know about the m_tedana option, it would be quite useful. :slight_smile:

Re: the field map:
When is the ‘best time’ to apply the distortion correction? After temporal (despike and tshift) and before structural (align etc.) steps or could it be done right at the start?
If the latter is ok to do, I could apply the correction and then simply run the pre-processing script as usual?
If the former way is the correct one, however, how could I use the matrix (which I assume is the ‘+ text file of params :’ that is mentioned in the epi_b0_correct.py output section) in a blip block?

By the way, here is my pre-processing script: https://www.dropbox.com/t/wHbjFpEO8VWK3uBp

Cheers,
Misho

Hi again Rick,

Just wanted to check if you’ve managed to come up with something. :smiley:

Cheers,
Misho

Hi there!

Is there any update on implementing epi_b0_correct.py in afni_proc.py? I am facing the same issue as Misho and am not sure if I should use it before running afni_proc.py or if it needs to be applied after a specific processing block within afni_proc.py.

Thanks in advanced for your time on this topic!

Best,
Sophie

hi Sophie,

pretty sure it hasn't been implemented in afni_proc yet, so you should do it beforehand. i did the AFNI parts for this recent paper and used the following script before my afni_proc (feel free to calculate "-scale_freq" using anything, and feel free to get rid of this whole loop structure:

foreach dir ( `ls -1d sub-*` )

	set subj = `echo $dir | awk -F "-" '{ print $2 }'`
	cd sub-${subj}/ses-01/run-01/

	set phasediffmap = sub-${subj}_phasediff.nii.gz
	set phasename = `basename -s .nii.gz $phasediffmap`
	set magnitudemap = sub-${subj}_magnitude1.nii.gz

	#assume all the BOLD epis have the same parameters
	set inepijson = /Users/sam/Desktop/vanderbilt/ds001796/jsons/sub-202_func_sub-202_task-rest_bold.json
	set inepi = BOLD_d.nii.gz

	#apparently freq maps need to be shorts not floats?
	3dcalc -datum short -a $phasediffmap -expr 'a' -prefix s${subj}_${phasename}_shrt.nii

	#make a mask with first echo magnitude map
	3dUnifize ${magnitudemap}'[0]'
	echo "skull stripping..."
	3dSkullStrip -input Unifized+orig. -prefix almostbin.nii
	3dcalc -a almostbin.nii -expr 'step(a)' -prefix ${subj}_mask_bin.nii
	\rm Unifized+orig* almostbin.nii

	#-scale_freq was calculated thusly:
	#mag1's TE - mag2's TE (or visa versa?) = 0.00738 - 0.00492 = 0.00246
	#ipython
	#import numpy as np
	#(2*np.pi)/(2*4096*0.00246)
	#0.31178471298

	epi_b0_correct.py -in_freq s${subj}_${phasename}_shrt.nii -in_epi_json $inepijson \
		-scale_freq 0.31178471298 -in_mask ${subj}_mask_bin.nii -in_epi $inepi 	\
		-blur_sigma 8 -prefix b0_corr_s$subj

	cd ../../..
end

note: this was Siemens data, and -blur_sigma is tweakable and makes a difference.

1 Like

Hi!

This is incredibly helpful and is exactly what I was looking for! Thank you so much for sharing!

Best,
Sophie

1 Like

Sorry for being behind on this, but indeed, it was never quite implemented in afni_proc.py (though it was really high on the todo list at a couple of points in time). It should not be too hard to add, so if there is interest, I could look at it again.

Note that it might be the same command being run before afni_proc.py, but the plan was to pass the resulting non-linear transformation matrix to AP, for it to combine with the other alignment transformations, basically replacing the AP block for reverse blip distortion correction.

Adding that to AP would reduce the effect of blurring with the additional transformation.

  • rick
1 Like