Anat-EPI Alignment Troubleshooting

[brainboi]

Hi AFNI team,

I’m having Anat_EPI alignment issues with most of our subjects. For troubleshooting purposes, I’ve focused on subjects that have the best initial alignment between anat and EPI. For these subjects, when running align_epi_anat.py independently, the output alignment looks excellent regardless of using lpc+ZZ or lpa+ZZ. But, the overlap between anat_final and final_epi_base_min_outlier is still noticeably off. Which leads me to believe something may be going awry during the 3dVolreg command? I thought the non-linear warp output of SSwarper looked good for the most part, but in case that could be contributing to this I’ll include one of those QC images too. AFNI ver is 23.0.01, Commodus, Jan 18 2023. Happy to provide anything else upon request! Attached is the output (APQC, etc) for one of our subjects for you to take a look at (I was going to include the whole thing but for upload size I took a screenshot of just the first half). Thank you for your continued help.

Sincerely,
Ian

# ================================= align ==================================
# for e2a: compute anat alignment transformation to EPI registration base
# (new anat will be current anatSS.${subj}+orig)
# run (localized) uniformity correction on EPI base
3dLocalUnifize -input vr_base_min_outlier+orig -prefix \
    vr_base_min_outlier_unif

align_epi_anat.py -anat2epi -anat anatSS.${subj}+orig   \
       -suffix _al_junk                                \
       -epi vr_base_min_outlier_unif+orig -epi_base 0  \
       -epi_strip 3dAutomask                           \
       -anat_has_skull no                              \
       -cost lpc+ZZ -check_flip                        \
       -volreg off -tshift off

# ================================== tlrc ==================================

# nothing to do: have external -tlrc_NL_warped_dsets

# warped anat     : anatQQ.${subj}+tlrc
# affine xform    : anatQQ.${subj}.aff12.1D
# non-linear warp : anatQQ.${subj}_WARP.nii

# ================================= volreg =================================
# align each dset to base volume, to anat, warp to tlrc space

# verify that we have a +tlrc warp dataset
if ( ! -f anatQQ.${subj}+tlrc.HEAD ) then
    echo "** missing +tlrc warp dataset: anatQQ.${subj}+tlrc.HEAD" 
    exit
endif

# register and warp
foreach run ( $runs )
    # register each volume to the base image
    3dvolreg -verbose -zpad 1 -base vr_base_min_outlier+orig             \
             -1Dfile dfile.r$run.1D -prefix rm.epi.volreg.r$run          \
             -cubic                                                      \
             -1Dmatrix_save mat.r$run.vr.aff12.1D                        \
             pb01.$subj.r$run.tshift+orig

    # create an all-1 dataset to mask the extents of the warp
    3dcalc -overwrite -a pb01.$subj.r$run.tshift+orig -expr 1            \
           -prefix rm.epi.all1

    # catenate volreg/epi2anat/tlrc xforms
    cat_matvec -ONELINE                                                  \
               anatQQ.${subj}.aff12.1D                                    \
               anatSS.${subj}_al_junk_mat.aff12.1D -I                     \
               mat.r$run.vr.aff12.1D > mat.r$run.warp.aff12.1D

    # apply catenated xform: volreg/epi2anat/tlrc/NLtlrc
    # then apply non-linear standard-space warp
    3dNwarpApply -master anatQQ.${subj}+tlrc -dxyz 3                      \
                 -source pb01.$subj.r$run.tshift+orig                    \
                 -nwarp "anatQQ.${subj}_WARP.nii mat.r$run.warp.aff12.1D" \
                 -prefix rm.epi.nomask.r$run

    # warp the all-1 dataset for extents masking 
    3dNwarpApply -master anatQQ.${subj}+tlrc -dxyz 3                      \
                 -source rm.epi.all1+orig                                \
                 -nwarp "anatQQ.${subj}_WARP.nii mat.r$run.warp.aff12.1D" \
                 -interp cubic                                           \
                 -ainterp NN -quiet                                      \
                 -prefix rm.epi.1.r$run

    # make an extents intersection mask of this run
    3dTstat -min -prefix rm.epi.min.r$run rm.epi.1.r$run+tlrc
end

# make a single file of registration params
cat dfile.r*.1D > dfile_rall.1D

# ----------------------------------------
# create the extents mask: mask_epi_extents+tlrc
# (this is a mask of voxels that have valid data at every TR)
3dMean -datum short -prefix rm.epi.mean rm.epi.min.r*.HEAD 
3dcalc -a rm.epi.mean+tlrc -expr 'step(a-0.999)' -prefix mask_epi_extents

# and apply the extents mask to the EPI data 
# (delete any time series with missing data)
foreach run ( $runs )
    3dcalc -a rm.epi.nomask.r$run+tlrc -b mask_epi_extents+tlrc          \
           -expr 'a*b' -prefix pb02.$subj.r$run.volreg
end

# warp the volreg base EPI dataset to make a final version
cat_matvec -ONELINE                                                      \
           anatQQ.${subj}.aff12.1D                                        \
           anatSS.${subj}_al_junk_mat.aff12.1D -I  > mat.basewarp.aff12.1D

3dNwarpApply -master anatQQ.${subj}+tlrc -dxyz 3                          \
             -source vr_base_min_outlier+orig                            \
             -nwarp "anatQQ.${subj}_WARP.nii mat.basewarp.aff12.1D"       \
             -prefix final_epi_vr_base_min_outlier

# create an anat_final dataset, aligned with stats
3dcopy anatQQ.${subj}+tlrc anat_final.$subj

# record final registration costs
3dAllineate -base final_epi_vr_base_min_outlier+tlrc -allcostX           \
            -input anat_final.$subj+tlrc |& tee out.allcostX.txt

# --------------------------------------
# create a TSNR dataset, just from run 1
3dTstat -mean -prefix rm.signal.vreg.r01 pb02.$subj.r01.volreg+tlrc
3dDetrend -polort 2 -prefix rm.noise.det -overwrite pb02.$subj.r01.volreg+tlrc
3dTstat -stdev -prefix rm.noise.vreg.r01 rm.noise.det+tlrc
3dcalc -a rm.signal.vreg.r01+tlrc                                        \
       -b rm.noise.vreg.r01+tlrc                                         \
       -c mask_epi_extents+tlrc                                          \
       -expr 'c*a/b' -prefix TSNR.vreg.r01.$subj

# -----------------------------------------
# warp anat follower datasets (non-linear)
# warp follower dataset copy_af_anat_orig+orig
3dNwarpApply -source copy_af_anat_orig+orig                              \
             -master anat_final.$subj+tlrc                               \
             -ainterp wsinc5 -nwarp anatQQ.${subj}_WARP.nii               \
             anatQQ.${subj}.aff12.1D                                      \
             -prefix follow_anat_anat_orig

[rickr]

Hi Ian,

Try adding the afni_proc.py option “-align_epi_strip_method 3dSkullStrip”.

It looks like 3dAutomask is not doing a good job of removing the skull area of the EPI (which is pretty apparent in your data), and 3dSkullStrip will probably do a better job.

• rick

[brainboi]

Hi Rick,

Thanks for the response! I tried 3dSkullStrip instead but the results look the same. In this case though, I do get this in the terminal output:

#++ removing skull or area outside brain
#Script is running (command trimmed):
  3dSkullStrip -orig_vol -input ./__tt_vr_base_min_outlier_unif_ts_rs+orig -prefix ./vr_base_min_outlier_unif_ts_rs_ns
oo     Warning 3dSkullStrip (SUMA_3dSkullStrip.c:1649):
Input dataset has a very low value range.
If stripping fails, repeat with option -fac 1000

At this suggestion, I added -skullstrip_opts -fac 1000 to my afni_proc.py command so the proc script looks like this:

align_epi_anat.py -anat2epi -anat anatSS.AB1056+orig        \
       -suffix _al_junk                                     \
       -epi vr_base_min_outlier_unif+orig -epi_base 0       \
       -epi_strip 3dSkullStrip                              \
       -anat_has_skull no                                   \
       -cost lpc+ZZ -check_flip -skullstrip_opts -fac 1000 \
       -volreg off -tshift off

However, I am still getting the same result and terminal output as before, so I’m unsure if 3dSkullstrip is applying the -fac 1000 option. Do you happen to have any other suggestions?

Sincerely,
Ian

[ptaylor]

Hi, Ian-

So, I assume your afni_proc.py (AP) command uses:

-align_unifize_epi        local                                  \

is that right? We recommend that now in most cases (at least for scans of humans) to deal with brightness inhomogeneity in the EPIs. Actually, the reason we don’t necessarily recommend it in animal scans at present is because the skulls/nonbrain material tend to be so much bigger and more difficult to remove, and the inhomogeneity doesn’t interact so well with that. Given the difficulty of auto-skullstripping here, I wonder if that kind of interaction is actually the problem (given discussion about trying 3dSkullStrip). But, for the moment, I still recommend keeping that option on, because of the inhomogeneity of the EPI.

What does your “align_opts_aea …” option look like? In this case, your initial EPI-anatomical alignment is quite close to start (the ‘vorig’ QC block image where obliquity is applied in the APQC HTML is the appropriate one to look at for judging that)—so, you probably shouldn’t need “-giant_move” or anything like that there. I might just start with:

-align_opts_aea           -cost lpc+ZZ                           \
                                      -check_flip                            \

If lpc+ZZ is not working well (which surprises me, but MRI datasets caaaan be surprising), then I would next try “-cost nmi”.

So, please try adjusting nmi, and then it would be helpful to know your other relevant EPI-anatomical AP options (or, perhaps just plopping the whole AP command for simplicity).

–pt

[brainboi]

Hi Paul,

Thanks for the reply, and sorry for the delay! I included the relevant portion of the proc script in my initial message if you’d like to have a look. I ended up removing the 3dLocalUnifize command and that has fixed the issue for the majority of our subjects! I’m not sure why that was causing the problem, but I’m glad the QC looks good now. And we’re using LPC+ZZ for the cost function which also looks great.

While I have this thread open, could you also take a look at this motion/regressor QC comparison between a low and high motion subject? We study a population that has higher motion, so we have set our motion and outlier frac limit to 1.0 and 0.15, respectively, and I want to double-check that this looks reasonable (I’m also not entirely sure how to interpret the global and radial correlate outputs). I’ve attached partial QC snapshots to this message.

Sincerely,
Ian

[ptaylor]

That is interesting about the EPI-anatomical alignment, and surprising. But if the alignment is working for you now, great.

Your censor thresholds, Enorm>1.0 and outlier_frac > 0.15, are higher than I have typically seen. From looking at the 1D profiles in the motion QC block, I would say there is some leftover motion artifact in the data, esp. in the first attachment. Some of that will manifest in higher corr_brain and radcor coverage, and this definitely appears in the 2nd attachment.

I think it would be worth checking out this very recent article on QCing FMRI data for some descriptions and examples:
Reynolds RC, Taylor PA, Glen DR (2023). Quality control practices in FMRI analysis: Philosophy, methods and examples using AFNI. Front. Neurosci. 16:1073800. doi: 10.3389/fnins.2022.1073800
https://www.frontiersin.org/articles/10.3389/fnins.2022.1073800/full/

–pt

[ptaylor]

As a follow-up, I might note that having Enorm censoring at 0.2 or 0.3 (which are in units of ~mm) would be more typical, and outlier frac of 0.05 (= finding spots where a brainmask has 5% or more outliers).

I think you still have plenty of degrees of freedom to give, and this shouldn’t result in too much time point loss, just eyeballing these motion plots. Hopefully that will reduce some of the widespread coverage in the various corr_brain and radcor images.

–pt

[brainboi]

Hi Paul,

First, thanks so much for the link to that incredible paper! I just finished reading it all through and it is so helpful.

Second, for our study we are looking at an often high-motion group (ADHD) and a (typically) lower-motion control group. We’ve deliberated extensively regarding the censor thresholds we should use, given the ADHD group has disproportionately higher levels of motion. Using the default of 0.3mm and 0.1 outlier fraction often brings the DOF below 50-60% for these subjects, hence why we adjusted them to 1mm and 0.15 (based on what we could find from other studies and looking at our own data). This does have the downside of missing some motion in subjects that don’t reach the higher threshold, but at least the outlier fraction tends to capture some of the motion still.

As far as GCOR, I tested some subjects by lowering the censor thresholds but the visual GCOR image tends to not change very much. And the radcor images look exactly the same regardless of the censor threshold. Even in some subjects with low motion and otherwise clean data there can be fairly high visual GCOR activity. Should this be a concern for our data? Following the lead of the paper you shared, would you suggest using a specific GCOR cutoff percentage for exclusionary purposes?

Sincerely,
Ian