Best Quality Alternatives to Bandpassing - Dataset Had NO Physiology .1D or .dat Info Collected

[sandra]

AFNI version info (afni -ver): Precompiled binary linux_ubuntu_16_64: May 8 2023 (Version AFNI_23.1.05 'Publius Helvius Pertinax')

To Paul Taylor – I really enjoyed your presentation on quality in Dallas at the Resting State Connectivity conference. I appreciate all your help!

AFNI,

Without cardiac and respiratory 1D or .dat files (older dataset) how can we improve our alignments and increase our degrees of freedom and overall quality on this resting fMRI dataset? We’re trying to remove (or greatly reduce) nuisance artifacts like motion, heart rate and respiration.

Our bandpass filters are set for
• HR = 60 – 90 bpm = 1 – 1.5 Hz
• BR = 10 – 14 = 0.166 – 0.234

Our dataset matches the specifics of TT_N27.
TR / sampling period =1.875
72 slices / acquisitions / temporal data points
TRs to remove = 7
04_BASELINE+orig = 4 min total resting data

Results: https://aisandra.com/research/6-baseline/index.html

afni_proc.py -subj_id resting-bandpass          \
             -dsets  04_BASELINE+orig.HEAD           \
             -blocks despike tshift align tlrc volreg blur mask regress \
             -copy_anat 14_SAG_T2_FLAIR_3D+orig                 \
             -tcat_remove_first_trs   7             \
             -tlrc_base TT_N27+tlrc                 \
             -tlrc_NL_warp                           \
             -volreg_align_to    MIN_OUTLIER                   \
             -volreg_align_e2a                               \
             -volreg_tlrc_warp                                \
             -blur_size 6.0                                \
             -regress_bandpass 0.01 0.1                       \
             -regress_bandpass 1.0 1.5                       \
             -regress_bandpass 0.166 0.234                       \
             -regress_motion_per_run                           \
             -regress_censor_motion 0.3                       \
             -regress_censor_outliers 0.1

[ptaylor]

Hi, Sandra-

It was nice to meet you at the conference and discuss the processing. Thanks for sending the full afni_proc.py (AP) command you are currently using. Some suggestions:

Firstly, as a reference have started to put more codes online for users when we/collaborators publish papers, to help comment on current recommendations. The homepage for these Code Examples (= Codex) is here.

From a recent paper on performing QC in FMRI (and part of a larger project of building a set of QC tutorials for the field across various software), we have this script, or more specifically the resting state AP command with options about the comments we use and why here.

From that, my primary comments would be to add these two options to your AP command, with the comments explaining what each does:

• Generally recommended starting point for EPI-anatomical alignment in human FMRI proc, where this lpc+ZZ cost function has been designed to align different contrast images (and for that last sub-option, left-right flipping can still occur...):

  -align_opts_aea           -cost lpc+ZZ -giant_move -check_flip 
  

• This option can help improve EPI-anatomical alignment, esp. if the EPI has brightness inhomogeneity (and it doesn't seem to hurt alignment even if that is not the case); generally recommended with human FMRI data processing nowadays:

  -align_unifize_epi        local 
  

• This adds useful APQC HTML items, radial correlation images of initial and volume-registered data (might see artifacts):

  -radial_correlate_blocks  tcat volreg regress
  

• Add a post-volreg TSNR plot to the APQC HTML:

  -volreg_compute_tsnr      yes
  

• Create useful mask from EPI-anatomical mask intersection (not applied to the EPI data here, but used to identify brain region):

  -mask_epi_anat            yes
  

• I might suggest a lower outlier censoring fraction, like 0.05 (=5%), just to be stricter on motion-related volumes:

  -regress_censor_outliers  0.05
  

• In resting state data, sometimes we include both motion regressors and their derivates as regressors, to help deal with some spikier noise aspects; to get the latter included in your processing, you could add this opt:

   -regress_apply_mot_types  demean deriv
  

• On a minor note, it is useful to state explicitly whether your input anatomical has a skull on still or not; I think the command is dealing with the input anatomical (that does have its skull currently; see a comment below about skullstripping before running AP...), but it would be clearer to add:

    -anat_has_skull           yes
  

Regarding bandpassing:

• First thing to note, the arguments supplied to -regress_bandpass .. are one or more pairs of frequency ranges to keep in your data (not to excluded, which it looks like your current syntax is doing). You have an FMRI sampling rate of TR=1.875 s.

• From a signal processing point of view, that means the Nyquist frequency of your data (the highest uniquely represented frequency in the signal) is: f_ny = 1/(2*1.875) = 0.267 Hz.

• Many people in resting state processing try to keep only low frequency fluctuations (LFFs), with the aim to exclude non-neuronal BOLD components like breathing and heartrate. While that makes some sense, and the range goes back to the original Biswal et al. 1995 paper, there are a lot of things lost in this argumentation. Namely, for the typical TRs used in MRI, most of those higher frequency heartrate will be aliased down into lower frequencies, and therefore be unremovable by bandpassing to LFFs (breathing miiiight be excludable in some cases, but breathing is also not constant); more importantly, perhaps, there is still useful information related to neuronal signals at higher frequencies, as Gohel and Biswal 2015 and others have shown. And super importantly, for standard TRs, that kind of bandpassing removes about 60% or more of the statistical degrees of freedom in the time series---greatly increasing the uncertainty bands of output results. This is somehow still not widely appreciated in the field, unfortunately, and more notes on it are here.

• Your estimated ranges of interest for heartrate (HR) and breathing rate (BR) were 1-1.5 Hz and 0.166-0.234 Hz, respectively.

• Since HR is totally above the Nyquist frequency, it cannot directly be removed by bandpassing: it will be aliased to other, lower frequencies below Nyquist, and in practice will be unable to be bandpassed out.

• Your estimated BR band, within which you would want to exclude frequencies, was in the interval 0.166-0.234 Hz. This range is already enirely above your selected LFF range, so you wouldn't need to include that as a separate band in the processing.

• With your data of f_ny = 0.26667 Hz, bandpassing to keep only 0.01-0.1Hz means you are removing 1-(0.1-0.01)/0.2667 = 66% of your degrees of freedom just from bandpassing---and more will be removed by censoring and the other regressors you have. That is a large cost.

• If you wanted to, perhaps you could keep a band in the interval 0.01-0.15 Hz, say, so you lose fewer degrees of freedom (only about 48% of them), while still removing the breathing-type frequencies. That might be strike a better balance, if you really want to try to exclude breathing-related frequencies that way.

Updating AP command:
So, if the above seem reasonable, then an updated AP command for your processing could look like:

afni_proc.py                                                                 \
    -subj_id                  resting-bandpass                               \
    -dsets                    04_BASELINE+orig.HEAD                          \
    -copy_anat                14_SAG_T2_FLAIR_3D+orig                        \
    -anat_has_skull           yes                                            \
    -blocks                   despike tshift align tlrc volreg blur mask     \
                              regress                                        \
    -radial_correlate_blocks  tcat volreg regress                            \
    -tcat_remove_first_trs    7                                              \
    -tlrc_base                TT_N27+tlrc                                    \
    -tlrc_NL_warp                                                            \
    -align_unifize_epi        local                                          \
    -align_opts_aea           -cost lpc+ZZ -giant_move -check_flip           \
    -volreg_compute_tsnr      yes                                            \
    -volreg_align_to          MIN_OUTLIER                                    \
    -volreg_align_e2a                                                        \
    -volreg_tlrc_warp                                                        \
    -blur_size                6.0                                            \
    -mask_epi_anat            yes                                            \
    -regress_bandpass         0.01 0.15                                      \
    -regress_motion_per_run                                                  \
    -regress_apply_mot_types  demean deriv                                   \
    -regress_censor_motion    0.3                                            \
    -regress_censor_outliers  0.05

Question:
Currently, you have 65 time points in input EPI/FMRI dataset (and you removed 7). However, for 4mins of acquired data with TR=1.875, you should have about 128 time points. What happened there?

Anatomical to template alignment:
The anatomical to template alignment looks like it should be better. I'm a bit surprised about that. Typically, we would recommend running @SSwarper prior to afni_proc.py, to perform skullstripping (SS) and nonlinear alignment (warping) at the same time. The results are then fed into afni_proc.py; this is demonstrated in the github repository I linked to above. We can look at that after some of the initial EPI-anatomical alignment and other options are implemented.

--pt

[sandra]

You’re right! I checked with my radiologist and this specific dataset is from a study of around TWO minutes not four! Good catch!

First, I ran the exact command you provided with “-anat_has_skull yes”. It had better alignment and brought the degrees of freedom from 17 to 23. It was still not high quality, so I ran this @SSwarper command:

@SSwarper -input 14_SAG_T2_FLAIR_3D+orig -base TT_N27_SSW.nii.gz -subid 6-baseline -odir sswarp-bp 2>&1 |tee output.SSwarper

Output pics, resulting quality page, and command output are all posted here: https://aisandra.com/research/6-baseline-SSwarpBP-noSkull/sswarper/index.html

Next, I updated all AFNI binaries and reran @SSwarper + afni_proc.py. The results were equivalent. Results: https://aisandra.com/research/6-baseline-SSwarpBP-newAfni/sswarp-newAfni/index.html

Questions:

• Is this the best atlas for our dataset? For afni_proc.py, we used TT_N27+tlrc and @SSwarper, we used TT_N27_SSW.nii.gz.
• If the atlas is an issue, is it worth creating custom atlases in Freesurfer and referencing them for each subject, or does using a custom atlas defeat the purpose of warping?
  • If so, is this automatable?
  • If not, is there anything else I can do here to improve quality?
• Would Freesurfer custom atlases be best in cases with tumors and other outlier-type populations?
• Would FreeSurfer recon-all (or any other Freesurfer command) help us here, as was done in the QC Practices in FMRI paper?

Our imaging:

• EPI = T2* (echoplanar)
• Anatomical = T2 weighted FLAIR

After that, I ran do_21_ap_rest_NL.tcsh in the directory containing data for subject 6:

tcsh -xef do_21_ap_rest_NL.tcsh |& tee output.do_21_ap_rest_NL

I received the following error:

source /etc/profile.d/modules.csh
/etc/profile.d/modules.csh: No such file or directory.

More questions:

• Should I download or install something before running this script?
• What should I modify within the script, or is it self-setting using shell or AFNI environment variables of some kind?
• Anything specific I should uncomment for our situation?
• Is there a template I should use to access my recommended atlas?
• Is this a script I run within a subject directory, or is it an entire modification to my AFNI install?

Even more questions:

How do I

• Get numerical values of
  • nodes associated with a specific seed correlation?
  • connectivity associated with multiple selected nodes?
  • nodes within a network and between networks?
  • of ALL seed correlation networks within the subject, for different individual networks?
• How do we display connectivity values between nodes?

Thanks again for all your work on quality. It’s a huge benefit to society as a whole.

[ptaylor]

Hi, Sandra-

That is awesome that you can share the various steps like that via the webpage---that helps a lot.

Below, there is a long reply, and actually it is divided into 2 pieces.

The following abbrevs are used to save needless wear and tear on me:

• AP = afni_proc.py
• SSW = @SSwarper
• FS = FreeSurfer

To answer your questions:

• That is a find template to use as a reference here; on a minor note, since you use TT_N27_SSW.nii.gz specifically within SSW, you should use that in the afni_proc.py command, too, instead of TT_N27+tlrc, but they should match, anyways.
• You can leave the final results of AP in subject anatomical space, yes; this is often most useful if you want to do an ROI-based analysis using output FS-estimated regions, either on the cortical surface or in the volume. You can pass the FS-estimated atlases/ROI maps as "follower datasets" to AP, so that they end up efficiently on the same final grid as the EPI data.
• In the case of input volumes with tumors, it is possible that staying in anatomical space might be useful, though I am not sure what that does to the quality of estimation of ROIs in the FS step.
• 3dQwarp, the AFNI program underlying SSW, actually takes an "exclusion mask" as input, to help when datasets being aligned have very different substructures. One example case for why it was added was to align a dset with a tumor to a template, as best as possible. SSW actually doesn't have this as part of its input, but I should be able to add that in relatively easily, if that would help.

---

Now, a separate point about the processing you have done so far and shown.

There are two pieces to this:

1. the SSW skullstripping+nonlinear alignment
2. The AP FMRI processing (making use of the SSW outputs)

There are a couple things to fix about the SSW warper part first, so let's focus on that.

The @SSwarper command by default assumes that the input is a T1w volume, which looks like one of the SSW template dsets (important because of how the default cost function works, to drive the alignment internally). That means having relative tissue brightness of (in increasing order): ventricles, GM, WM. For example, that is what the TT_N27_SSW reference dset looks like, as in the first column here:
https://afni.nimh.nih.gov/pub/dist/doc/htmldoc/_images/ALL_TT_N27_SSW.jpg
However, while your input anatomical FLAIR dset looks nice and sharp, it has different relative tissue brightness, which is (also in increasing order): ventricle, WM, GM.

This means that I would suspect that a different cost function than the default one be used, and the proof of this is shown in this QC image from your SSW output:

... where the overlaid edges of the template don't match well with the underlaid anatomical structures (and this also shows the relative tissue brightness of the input, that I had mentioned above).

On a sidenote, the skullstripping looks generally excellent here still:

https://aisandra.com/research/6-baseline-SSwarpBP-newAfni/sswarp-newAfni/QC_anatSS.6-baseline.jpg(image larger than 4 MB)

Anyways, to address the anatomical alignment for this kind of input dataset contrast, I think this should be re-run with a different cost function specified. We also note that alignment happens in two stages here: an initial, global affine fit, and then the local patch refinement. We will specify the cost functions for each. I don't have experience with this kind of input anatomical contrast, so I am not sure which of the following might work best. I suggest running each of these combinations (sorry for the extra work---but that is the prices of trailblazing!):

• Trial run A, use both:

-cost_nl_init   lpa+ZZ
-cost_nl_final  lpa

• Trial run B, use both:

-cost_nl_init   nmi
-cost_nl_final  nmi

When we look more at the EPI-anatomical alignment, as well, the fact that this is a FLAIR and not a T1w volume might lead to cost function adjustment, too.

And again, let me know if the anatomical-with-tumor-to-template alignment is a practical use case for your analyses.

--pt